Actual source code: bddc.c
1: /* TODOLIST
3: Solvers
4: - Add support for cholesky for coarse solver (similar to local solvers)
5: - Propagate ksp prefixes for solvers to mat objects?
7: User interface
8: - ** DM attached to pc?
10: Debugging output
11: - * Better management of verbosity levels of debugging output
13: Extra
14: - *** Is it possible to work with PCBDDCGraph on boundary indices only (less memory consumed)?
15: - BDDC with MG framework?
17: MATIS related operations contained in BDDC code
18: - Provide general case for subassembling
20: */
22: #include <petsc/private/pcbddcimpl.h>
23: #include <petsc/private/pcbddcprivateimpl.h>
24: #include <petscblaslapack.h>
26: static PetscBool PCBDDCPackageInitialized = PETSC_FALSE;
28: static PetscBool cited = PETSC_FALSE;
29: static const char citation[] = "@article{ZampiniPCBDDC,\n"
30: "author = {Stefano Zampini},\n"
31: "title = {{PCBDDC}: A Class of Robust Dual-Primal Methods in {PETS}c},\n"
32: "journal = {SIAM Journal on Scientific Computing},\n"
33: "volume = {38},\n"
34: "number = {5},\n"
35: "pages = {S282-S306},\n"
36: "year = {2016},\n"
37: "doi = {10.1137/15M1025785},\n"
38: "URL = {http://dx.doi.org/10.1137/15M1025785},\n"
39: "eprint = {http://dx.doi.org/10.1137/15M1025785}\n"
40: "}\n";
42: PetscLogEvent PC_BDDC_Topology[PETSC_PCBDDC_MAXLEVELS];
43: PetscLogEvent PC_BDDC_LocalSolvers[PETSC_PCBDDC_MAXLEVELS];
44: PetscLogEvent PC_BDDC_LocalWork[PETSC_PCBDDC_MAXLEVELS];
45: PetscLogEvent PC_BDDC_CorrectionSetUp[PETSC_PCBDDC_MAXLEVELS];
46: PetscLogEvent PC_BDDC_ApproxSetUp[PETSC_PCBDDC_MAXLEVELS];
47: PetscLogEvent PC_BDDC_ApproxApply[PETSC_PCBDDC_MAXLEVELS];
48: PetscLogEvent PC_BDDC_CoarseSetUp[PETSC_PCBDDC_MAXLEVELS];
49: PetscLogEvent PC_BDDC_CoarseSolver[PETSC_PCBDDC_MAXLEVELS];
50: PetscLogEvent PC_BDDC_AdaptiveSetUp[PETSC_PCBDDC_MAXLEVELS];
51: PetscLogEvent PC_BDDC_Scaling[PETSC_PCBDDC_MAXLEVELS];
52: PetscLogEvent PC_BDDC_Schurs[PETSC_PCBDDC_MAXLEVELS];
53: PetscLogEvent PC_BDDC_Solves[PETSC_PCBDDC_MAXLEVELS][3];
55: const char *const PCBDDCInterfaceExtTypes[] = {"DIRICHLET", "LUMP", "PCBDDCInterfaceExtType", "PC_BDDC_INTERFACE_EXT_", NULL};
57: static PetscErrorCode PCApply_BDDC(PC, Vec, Vec);
59: static PetscErrorCode PCSetFromOptions_BDDC(PC pc, PetscOptionItems *PetscOptionsObject)
60: {
61: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
62: PetscInt nt, i;
64: PetscFunctionBegin;
65: PetscOptionsHeadBegin(PetscOptionsObject, "BDDC options");
66: /* Verbose debugging */
67: PetscCall(PetscOptionsInt("-pc_bddc_check_level", "Verbose output for PCBDDC (intended for debug)", "none", pcbddc->dbg_flag, &pcbddc->dbg_flag, NULL));
68: /* Approximate solvers */
69: PetscCall(PetscOptionsEnum("-pc_bddc_interface_ext_type", "Use DIRICHLET or LUMP to extend interface corrections to interior", "PCBDDCSetInterfaceExtType", PCBDDCInterfaceExtTypes, (PetscEnum)pcbddc->interface_extension, (PetscEnum *)&pcbddc->interface_extension, NULL));
70: if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_DIRICHLET) {
71: PetscCall(PetscOptionsBool("-pc_bddc_dirichlet_approximate", "Inform PCBDDC that we are using approximate Dirichlet solvers", "none", pcbddc->NullSpace_corr[0], &pcbddc->NullSpace_corr[0], NULL));
72: PetscCall(PetscOptionsBool("-pc_bddc_dirichlet_approximate_scale", "Inform PCBDDC that we need to scale the Dirichlet solve", "none", pcbddc->NullSpace_corr[1], &pcbddc->NullSpace_corr[1], NULL));
73: } else {
74: /* This flag is needed/implied by lumping */
75: pcbddc->switch_static = PETSC_TRUE;
76: }
77: PetscCall(PetscOptionsBool("-pc_bddc_neumann_approximate", "Inform PCBDDC that we are using approximate Neumann solvers", "none", pcbddc->NullSpace_corr[2], &pcbddc->NullSpace_corr[2], NULL));
78: PetscCall(PetscOptionsBool("-pc_bddc_neumann_approximate_scale", "Inform PCBDDC that we need to scale the Neumann solve", "none", pcbddc->NullSpace_corr[3], &pcbddc->NullSpace_corr[3], NULL));
79: /* Primal space customization */
80: PetscCall(PetscOptionsBool("-pc_bddc_use_local_mat_graph", "Use or not adjacency graph of local mat for interface analysis", "none", pcbddc->use_local_adj, &pcbddc->use_local_adj, NULL));
81: PetscCall(PetscOptionsInt("-pc_bddc_graph_maxcount", "Maximum number of shared subdomains for a connected component", "none", pcbddc->graphmaxcount, &pcbddc->graphmaxcount, NULL));
82: PetscCall(PetscOptionsBool("-pc_bddc_corner_selection", "Activates face-based corner selection", "none", pcbddc->corner_selection, &pcbddc->corner_selection, NULL));
83: PetscCall(PetscOptionsBool("-pc_bddc_use_vertices", "Use or not corner dofs in coarse space", "none", pcbddc->use_vertices, &pcbddc->use_vertices, NULL));
84: PetscCall(PetscOptionsBool("-pc_bddc_use_edges", "Use or not edge constraints in coarse space", "none", pcbddc->use_edges, &pcbddc->use_edges, NULL));
85: PetscCall(PetscOptionsBool("-pc_bddc_use_faces", "Use or not face constraints in coarse space", "none", pcbddc->use_faces, &pcbddc->use_faces, NULL));
86: PetscCall(PetscOptionsInt("-pc_bddc_vertex_size", "Connected components smaller or equal to vertex size will be considered as primal vertices", "none", pcbddc->vertex_size, &pcbddc->vertex_size, NULL));
87: PetscCall(PetscOptionsBool("-pc_bddc_use_nnsp", "Use near null space attached to the matrix to compute constraints", "none", pcbddc->use_nnsp, &pcbddc->use_nnsp, NULL));
88: PetscCall(PetscOptionsBool("-pc_bddc_use_nnsp_true", "Use near null space attached to the matrix to compute constraints as is", "none", pcbddc->use_nnsp_true, &pcbddc->use_nnsp_true, NULL));
89: PetscCall(PetscOptionsBool("-pc_bddc_use_qr_single", "Use QR factorization for single constraints on cc (QR is always used when multiple constraints are present)", "none", pcbddc->use_qr_single, &pcbddc->use_qr_single, NULL));
90: /* Change of basis */
91: PetscCall(PetscOptionsBool("-pc_bddc_use_change_of_basis", "Use or not internal change of basis on local edge nodes", "none", pcbddc->use_change_of_basis, &pcbddc->use_change_of_basis, NULL));
92: PetscCall(PetscOptionsBool("-pc_bddc_use_change_on_faces", "Use or not internal change of basis on local face nodes", "none", pcbddc->use_change_on_faces, &pcbddc->use_change_on_faces, NULL));
93: if (!pcbddc->use_change_of_basis) pcbddc->use_change_on_faces = PETSC_FALSE;
94: /* Switch between M_2 (default) and M_3 preconditioners (as defined by C. Dohrmann in the ref. article) */
95: PetscCall(PetscOptionsBool("-pc_bddc_switch_static", "Switch on static condensation ops around the interface preconditioner", "none", pcbddc->switch_static, &pcbddc->switch_static, NULL));
96: PetscCall(PetscOptionsInt("-pc_bddc_coarse_eqs_per_proc", "Target number of equations per process for coarse problem redistribution (significant only at the coarsest level)", "none", pcbddc->coarse_eqs_per_proc, &pcbddc->coarse_eqs_per_proc, NULL));
97: i = pcbddc->coarsening_ratio;
98: PetscCall(PetscOptionsInt("-pc_bddc_coarsening_ratio", "Set coarsening ratio used in multilevel coarsening", "PCBDDCSetCoarseningRatio", i, &i, NULL));
99: PetscCall(PCBDDCSetCoarseningRatio(pc, i));
100: i = pcbddc->max_levels;
101: PetscCall(PetscOptionsInt("-pc_bddc_levels", "Set maximum number of levels for multilevel", "PCBDDCSetLevels", i, &i, NULL));
102: PetscCall(PCBDDCSetLevels(pc, i));
103: PetscCall(PetscOptionsInt("-pc_bddc_coarse_eqs_limit", "Set maximum number of equations on coarsest grid to aim for", "none", pcbddc->coarse_eqs_limit, &pcbddc->coarse_eqs_limit, NULL));
104: PetscCall(PetscOptionsBool("-pc_bddc_use_coarse_estimates", "Use estimated eigenvalues for coarse problem", "none", pcbddc->use_coarse_estimates, &pcbddc->use_coarse_estimates, NULL));
105: PetscCall(PetscOptionsBool("-pc_bddc_use_deluxe_scaling", "Use deluxe scaling for BDDC", "none", pcbddc->use_deluxe_scaling, &pcbddc->use_deluxe_scaling, NULL));
106: PetscCall(PetscOptionsBool("-pc_bddc_schur_rebuild", "Whether or not the interface graph for Schur principal minors has to be rebuilt (i.e. define the interface without any adjacency)", "none", pcbddc->sub_schurs_rebuild, &pcbddc->sub_schurs_rebuild, NULL));
107: PetscCall(PetscOptionsInt("-pc_bddc_schur_layers", "Number of dofs' layers for the computation of principal minors (i.e. -1 uses all dofs)", "none", pcbddc->sub_schurs_layers, &pcbddc->sub_schurs_layers, NULL));
108: PetscCall(PetscOptionsBool("-pc_bddc_schur_use_useradj", "Whether or not the CSR graph specified by the user should be used for computing successive layers (default is to use adj of local mat)", "none", pcbddc->sub_schurs_use_useradj, &pcbddc->sub_schurs_use_useradj, NULL));
109: PetscCall(PetscOptionsBool("-pc_bddc_schur_exact", "Whether or not to use the exact Schur complement instead of the reduced one (which excludes size 1 cc)", "none", pcbddc->sub_schurs_exact_schur, &pcbddc->sub_schurs_exact_schur, NULL));
110: PetscCall(PetscOptionsBool("-pc_bddc_deluxe_zerorows", "Zero rows and columns of deluxe operators associated with primal dofs", "none", pcbddc->deluxe_zerorows, &pcbddc->deluxe_zerorows, NULL));
111: PetscCall(PetscOptionsBool("-pc_bddc_deluxe_singlemat", "Collapse deluxe operators", "none", pcbddc->deluxe_singlemat, &pcbddc->deluxe_singlemat, NULL));
112: PetscCall(PetscOptionsBool("-pc_bddc_adaptive_userdefined", "Use user-defined constraints (should be attached via MatSetNearNullSpace to pmat) in addition to those adaptively generated", "none", pcbddc->adaptive_userdefined, &pcbddc->adaptive_userdefined, NULL));
113: nt = 2;
114: PetscCall(PetscOptionsRealArray("-pc_bddc_adaptive_threshold", "Thresholds to be used for adaptive selection of constraints", "none", pcbddc->adaptive_threshold, &nt, NULL));
115: if (nt == 1) pcbddc->adaptive_threshold[1] = pcbddc->adaptive_threshold[0];
116: PetscCall(PetscOptionsInt("-pc_bddc_adaptive_nmin", "Minimum number of constraints per connected components", "none", pcbddc->adaptive_nmin, &pcbddc->adaptive_nmin, NULL));
117: PetscCall(PetscOptionsInt("-pc_bddc_adaptive_nmax", "Maximum number of constraints per connected components", "none", pcbddc->adaptive_nmax, &pcbddc->adaptive_nmax, NULL));
118: PetscCall(PetscOptionsBool("-pc_bddc_symmetric", "Symmetric computation of primal basis functions", "none", pcbddc->symmetric_primal, &pcbddc->symmetric_primal, NULL));
119: PetscCall(PetscOptionsInt("-pc_bddc_coarse_adj", "Number of processors where to map the coarse adjacency list", "none", pcbddc->coarse_adj_red, &pcbddc->coarse_adj_red, NULL));
120: PetscCall(PetscOptionsBool("-pc_bddc_benign_trick", "Apply the benign subspace trick to saddle point problems with discontinuous pressures", "none", pcbddc->benign_saddle_point, &pcbddc->benign_saddle_point, NULL));
121: PetscCall(PetscOptionsBool("-pc_bddc_benign_change", "Compute the pressure change of basis explicitly", "none", pcbddc->benign_change_explicit, &pcbddc->benign_change_explicit, NULL));
122: PetscCall(PetscOptionsBool("-pc_bddc_benign_compute_correction", "Compute the benign correction during PreSolve", "none", pcbddc->benign_compute_correction, &pcbddc->benign_compute_correction, NULL));
123: PetscCall(PetscOptionsBool("-pc_bddc_nonetflux", "Automatic computation of no-net-flux quadrature weights", "none", pcbddc->compute_nonetflux, &pcbddc->compute_nonetflux, NULL));
124: PetscCall(PetscOptionsBool("-pc_bddc_detect_disconnected", "Detects disconnected subdomains", "none", pcbddc->detect_disconnected, &pcbddc->detect_disconnected, NULL));
125: PetscCall(PetscOptionsBool("-pc_bddc_detect_disconnected_filter", "Filters out small entries in the local matrix when detecting disconnected subdomains", "none", pcbddc->detect_disconnected_filter, &pcbddc->detect_disconnected_filter, NULL));
126: PetscCall(PetscOptionsBool("-pc_bddc_eliminate_dirichlet", "Whether or not we want to eliminate dirichlet dofs during presolve", "none", pcbddc->eliminate_dirdofs, &pcbddc->eliminate_dirdofs, NULL));
127: PetscOptionsHeadEnd();
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: static PetscErrorCode PCView_BDDC(PC pc, PetscViewer viewer)
132: {
133: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
134: PC_IS *pcis = (PC_IS *)pc->data;
135: PetscBool isascii;
136: PetscSubcomm subcomm;
137: PetscViewer subviewer;
139: PetscFunctionBegin;
140: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
141: /* ASCII viewer */
142: if (isascii) {
143: PetscMPIInt color, rank, size;
144: PetscInt64 loc[7], gsum[6], gmax[6], gmin[6], totbenign;
145: PetscScalar interface_size;
146: PetscReal ratio1 = 0., ratio2 = 0.;
147: Vec counter;
149: if (!pc->setupcalled) PetscCall(PetscViewerASCIIPrintf(viewer, " Partial information available: preconditioner has not been setup yet\n"));
150: PetscCall(PetscViewerASCIIPrintf(viewer, " Use verbose output: %" PetscInt_FMT "\n", pcbddc->dbg_flag));
151: PetscCall(PetscViewerASCIIPrintf(viewer, " Use user-defined CSR: %d\n", !!pcbddc->mat_graph->nvtxs_csr));
152: PetscCall(PetscViewerASCIIPrintf(viewer, " Use local mat graph: %d\n", pcbddc->use_local_adj && !pcbddc->mat_graph->nvtxs_csr));
153: if (pcbddc->mat_graph->twodim) {
154: PetscCall(PetscViewerASCIIPrintf(viewer, " Connectivity graph topological dimension: 2\n"));
155: } else {
156: PetscCall(PetscViewerASCIIPrintf(viewer, " Connectivity graph topological dimension: 3\n"));
157: }
158: if (pcbddc->graphmaxcount != PETSC_MAX_INT) PetscCall(PetscViewerASCIIPrintf(viewer, " Graph max count: %" PetscInt_FMT "\n", pcbddc->graphmaxcount));
159: PetscCall(PetscViewerASCIIPrintf(viewer, " Corner selection: %d (selected %d)\n", pcbddc->corner_selection, pcbddc->corner_selected));
160: PetscCall(PetscViewerASCIIPrintf(viewer, " Use vertices: %d (vertex size %" PetscInt_FMT ")\n", pcbddc->use_vertices, pcbddc->vertex_size));
161: PetscCall(PetscViewerASCIIPrintf(viewer, " Use edges: %d\n", pcbddc->use_edges));
162: PetscCall(PetscViewerASCIIPrintf(viewer, " Use faces: %d\n", pcbddc->use_faces));
163: PetscCall(PetscViewerASCIIPrintf(viewer, " Use true near null space: %d\n", pcbddc->use_nnsp_true));
164: PetscCall(PetscViewerASCIIPrintf(viewer, " Use QR for single constraints on cc: %d\n", pcbddc->use_qr_single));
165: PetscCall(PetscViewerASCIIPrintf(viewer, " Use change of basis on local edge nodes: %d\n", pcbddc->use_change_of_basis));
166: PetscCall(PetscViewerASCIIPrintf(viewer, " Use change of basis on local face nodes: %d\n", pcbddc->use_change_on_faces));
167: PetscCall(PetscViewerASCIIPrintf(viewer, " User defined change of basis matrix: %d\n", !!pcbddc->user_ChangeOfBasisMatrix));
168: PetscCall(PetscViewerASCIIPrintf(viewer, " Has change of basis matrix: %d\n", !!pcbddc->ChangeOfBasisMatrix));
169: PetscCall(PetscViewerASCIIPrintf(viewer, " Eliminate dirichlet boundary dofs: %d\n", pcbddc->eliminate_dirdofs));
170: PetscCall(PetscViewerASCIIPrintf(viewer, " Switch on static condensation ops around the interface preconditioner: %d\n", pcbddc->switch_static));
171: PetscCall(PetscViewerASCIIPrintf(viewer, " Use exact dirichlet trick: %d\n", pcbddc->use_exact_dirichlet_trick));
172: PetscCall(PetscViewerASCIIPrintf(viewer, " Interface extension: %s\n", PCBDDCInterfaceExtTypes[pcbddc->interface_extension]));
173: PetscCall(PetscViewerASCIIPrintf(viewer, " Multilevel max levels: %" PetscInt_FMT "\n", pcbddc->max_levels));
174: PetscCall(PetscViewerASCIIPrintf(viewer, " Multilevel coarsening ratio: %" PetscInt_FMT "\n", pcbddc->coarsening_ratio));
175: PetscCall(PetscViewerASCIIPrintf(viewer, " Use estimated eigs for coarse problem: %d\n", pcbddc->use_coarse_estimates));
176: PetscCall(PetscViewerASCIIPrintf(viewer, " Use deluxe scaling: %d\n", pcbddc->use_deluxe_scaling));
177: PetscCall(PetscViewerASCIIPrintf(viewer, " Use deluxe zerorows: %d\n", pcbddc->deluxe_zerorows));
178: PetscCall(PetscViewerASCIIPrintf(viewer, " Use deluxe singlemat: %d\n", pcbddc->deluxe_singlemat));
179: PetscCall(PetscViewerASCIIPrintf(viewer, " Rebuild interface graph for Schur principal minors: %d\n", pcbddc->sub_schurs_rebuild));
180: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of dofs' layers for the computation of principal minors: %" PetscInt_FMT "\n", pcbddc->sub_schurs_layers));
181: PetscCall(PetscViewerASCIIPrintf(viewer, " Use user CSR graph to compute successive layers: %d\n", pcbddc->sub_schurs_use_useradj));
182: if (pcbddc->adaptive_threshold[1] != pcbddc->adaptive_threshold[0]) {
183: PetscCall(PetscViewerASCIIPrintf(viewer, " Adaptive constraint selection thresholds (active %d, userdefined %d): %g,%g\n", pcbddc->adaptive_selection, pcbddc->adaptive_userdefined, (double)pcbddc->adaptive_threshold[0], (double)pcbddc->adaptive_threshold[1]));
184: } else {
185: PetscCall(PetscViewerASCIIPrintf(viewer, " Adaptive constraint selection threshold (active %d, userdefined %d): %g\n", pcbddc->adaptive_selection, pcbddc->adaptive_userdefined, (double)pcbddc->adaptive_threshold[0]));
186: }
187: PetscCall(PetscViewerASCIIPrintf(viewer, " Min constraints / connected component: %" PetscInt_FMT "\n", pcbddc->adaptive_nmin));
188: PetscCall(PetscViewerASCIIPrintf(viewer, " Max constraints / connected component: %" PetscInt_FMT "\n", pcbddc->adaptive_nmax));
189: PetscCall(PetscViewerASCIIPrintf(viewer, " Invert exact Schur complement for adaptive selection: %d\n", pcbddc->sub_schurs_exact_schur));
190: PetscCall(PetscViewerASCIIPrintf(viewer, " Symmetric computation of primal basis functions: %d\n", pcbddc->symmetric_primal));
191: PetscCall(PetscViewerASCIIPrintf(viewer, " Num. Procs. to map coarse adjacency list: %" PetscInt_FMT "\n", pcbddc->coarse_adj_red));
192: PetscCall(PetscViewerASCIIPrintf(viewer, " Coarse eqs per proc (significant at the coarsest level): %" PetscInt_FMT "\n", pcbddc->coarse_eqs_per_proc));
193: PetscCall(PetscViewerASCIIPrintf(viewer, " Detect disconnected: %d (filter %d)\n", pcbddc->detect_disconnected, pcbddc->detect_disconnected_filter));
194: PetscCall(PetscViewerASCIIPrintf(viewer, " Benign subspace trick: %d (change explicit %d)\n", pcbddc->benign_saddle_point, pcbddc->benign_change_explicit));
195: PetscCall(PetscViewerASCIIPrintf(viewer, " Benign subspace trick is active: %d\n", pcbddc->benign_have_null));
196: PetscCall(PetscViewerASCIIPrintf(viewer, " Algebraic computation of no-net-flux: %d\n", pcbddc->compute_nonetflux));
197: if (!pc->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
199: /* compute interface size */
200: PetscCall(VecSet(pcis->vec1_B, 1.0));
201: PetscCall(MatCreateVecs(pc->pmat, &counter, NULL));
202: PetscCall(VecSet(counter, 0.0));
203: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, counter, INSERT_VALUES, SCATTER_REVERSE));
204: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, counter, INSERT_VALUES, SCATTER_REVERSE));
205: PetscCall(VecSum(counter, &interface_size));
206: PetscCall(VecDestroy(&counter));
208: /* compute some statistics on the domain decomposition */
209: gsum[0] = 1;
210: gsum[1] = gsum[2] = gsum[3] = gsum[4] = gsum[5] = 0;
211: loc[0] = !!pcis->n;
212: loc[1] = pcis->n - pcis->n_B;
213: loc[2] = pcis->n_B;
214: loc[3] = pcbddc->local_primal_size;
215: loc[4] = pcis->n;
216: loc[5] = pcbddc->n_local_subs > 0 ? pcbddc->n_local_subs : (pcis->n ? 1 : 0);
217: loc[6] = pcbddc->benign_n;
218: PetscCallMPI(MPI_Reduce(loc, gsum, 6, MPIU_INT64, MPI_SUM, 0, PetscObjectComm((PetscObject)pc)));
219: if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = loc[5] = -1;
220: PetscCallMPI(MPI_Reduce(loc, gmax, 6, MPIU_INT64, MPI_MAX, 0, PetscObjectComm((PetscObject)pc)));
221: if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = loc[5] = PETSC_MAX_INT;
222: PetscCallMPI(MPI_Reduce(loc, gmin, 6, MPIU_INT64, MPI_MIN, 0, PetscObjectComm((PetscObject)pc)));
223: PetscCallMPI(MPI_Reduce(&loc[6], &totbenign, 1, MPIU_INT64, MPI_SUM, 0, PetscObjectComm((PetscObject)pc)));
224: if (pcbddc->coarse_size) {
225: ratio1 = pc->pmat->rmap->N / (1. * pcbddc->coarse_size);
226: ratio2 = PetscRealPart(interface_size) / pcbddc->coarse_size;
227: }
228: PetscCall(PetscViewerASCIIPrintf(viewer, "********************************** STATISTICS AT LEVEL %" PetscInt_FMT " **********************************\n", pcbddc->current_level));
229: PetscCall(PetscViewerASCIIPrintf(viewer, " Global dofs sizes: all %" PetscInt_FMT " interface %" PetscInt_FMT " coarse %" PetscInt_FMT "\n", pc->pmat->rmap->N, (PetscInt)PetscRealPart(interface_size), pcbddc->coarse_size));
230: PetscCall(PetscViewerASCIIPrintf(viewer, " Coarsening ratios: all/coarse %" PetscInt_FMT " interface/coarse %" PetscInt_FMT "\n", (PetscInt)ratio1, (PetscInt)ratio2));
231: PetscCall(PetscViewerASCIIPrintf(viewer, " Active processes : %" PetscInt_FMT "\n", (PetscInt)gsum[0]));
232: PetscCall(PetscViewerASCIIPrintf(viewer, " Total subdomains : %" PetscInt_FMT "\n", (PetscInt)gsum[5]));
233: if (pcbddc->benign_have_null) PetscCall(PetscViewerASCIIPrintf(viewer, " Benign subs : %" PetscInt_FMT "\n", (PetscInt)totbenign));
234: PetscCall(PetscViewerASCIIPrintf(viewer, " Dofs type :\tMIN\tMAX\tMEAN\n"));
235: PetscCall(PetscViewerASCIIPrintf(viewer, " Interior dofs :\t%" PetscInt_FMT "\t%" PetscInt_FMT "\t%" PetscInt_FMT "\n", (PetscInt)gmin[1], (PetscInt)gmax[1], (PetscInt)(gsum[1] / gsum[0])));
236: PetscCall(PetscViewerASCIIPrintf(viewer, " Interface dofs :\t%" PetscInt_FMT "\t%" PetscInt_FMT "\t%" PetscInt_FMT "\n", (PetscInt)gmin[2], (PetscInt)gmax[2], (PetscInt)(gsum[2] / gsum[0])));
237: PetscCall(PetscViewerASCIIPrintf(viewer, " Primal dofs :\t%" PetscInt_FMT "\t%" PetscInt_FMT "\t%" PetscInt_FMT "\n", (PetscInt)gmin[3], (PetscInt)gmax[3], (PetscInt)(gsum[3] / gsum[0])));
238: PetscCall(PetscViewerASCIIPrintf(viewer, " Local dofs :\t%" PetscInt_FMT "\t%" PetscInt_FMT "\t%" PetscInt_FMT "\n", (PetscInt)gmin[4], (PetscInt)gmax[4], (PetscInt)(gsum[4] / gsum[0])));
239: PetscCall(PetscViewerASCIIPrintf(viewer, " Local subs :\t%" PetscInt_FMT "\t%" PetscInt_FMT "\n", (PetscInt)gmin[5], (PetscInt)gmax[5]));
240: PetscCall(PetscViewerFlush(viewer));
242: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
244: /* local solvers */
245: PetscCall(PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)pcbddc->ksp_D), &subviewer));
246: if (rank == 0) {
247: PetscCall(PetscViewerASCIIPrintf(subviewer, "--- Interior solver (rank 0)\n"));
248: PetscCall(PetscViewerASCIIPushTab(subviewer));
249: PetscCall(KSPView(pcbddc->ksp_D, subviewer));
250: PetscCall(PetscViewerASCIIPopTab(subviewer));
251: PetscCall(PetscViewerASCIIPrintf(subviewer, "--- Correction solver (rank 0)\n"));
252: PetscCall(PetscViewerASCIIPushTab(subviewer));
253: PetscCall(KSPView(pcbddc->ksp_R, subviewer));
254: PetscCall(PetscViewerASCIIPopTab(subviewer));
255: PetscCall(PetscViewerFlush(subviewer));
256: }
257: PetscCall(PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)pcbddc->ksp_D), &subviewer));
258: PetscCall(PetscViewerFlush(viewer));
260: /* the coarse problem can be handled by a different communicator */
261: if (pcbddc->coarse_ksp) color = 1;
262: else color = 0;
263: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
264: PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)pc), &subcomm));
265: PetscCall(PetscSubcommSetNumber(subcomm, PetscMin(size, 2)));
266: PetscCall(PetscSubcommSetTypeGeneral(subcomm, color, rank));
267: PetscCall(PetscViewerGetSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer));
268: if (color == 1) {
269: PetscCall(PetscViewerASCIIPrintf(subviewer, "--- Coarse solver\n"));
270: PetscCall(PetscViewerASCIIPushTab(subviewer));
271: PetscCall(KSPView(pcbddc->coarse_ksp, subviewer));
272: PetscCall(PetscViewerASCIIPopTab(subviewer));
273: PetscCall(PetscViewerFlush(subviewer));
274: }
275: PetscCall(PetscViewerRestoreSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer));
276: PetscCall(PetscSubcommDestroy(&subcomm));
277: PetscCall(PetscViewerFlush(viewer));
278: }
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: static PetscErrorCode PCBDDCSetDiscreteGradient_BDDC(PC pc, Mat G, PetscInt order, PetscInt field, PetscBool global, PetscBool conforming)
283: {
284: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
286: PetscFunctionBegin;
287: PetscCall(PetscObjectReference((PetscObject)G));
288: PetscCall(MatDestroy(&pcbddc->discretegradient));
289: pcbddc->discretegradient = G;
290: pcbddc->nedorder = order > 0 ? order : -order;
291: pcbddc->nedfield = field;
292: pcbddc->nedglobal = global;
293: pcbddc->conforming = conforming;
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: /*@
298: PCBDDCSetDiscreteGradient - Sets the discrete gradient to be used by the `PCBDDC` preconditioner
300: Collective
302: Input Parameters:
303: + pc - the preconditioning context
304: . G - the discrete gradient matrix (in `MATAIJ` format)
305: . order - the order of the Nedelec space (1 for the lowest order)
306: . field - the field id of the Nedelec dofs (not used if the fields have not been specified)
307: . global - the type of global ordering for the rows of `G`
308: - conforming - whether the mesh is conforming or not
310: Level: advanced
312: Note:
313: The discrete gradient matrix `G` is used to analyze the subdomain edges, and it should not contain any zero entry.
314: For variable order spaces, the order should be set to zero.
315: If `global` is `PETSC_TRUE`, the rows of `G` should be given in global ordering for the whole dofs;
316: if `PETSC_FALSE`, the ordering should be global for the Nedelec field.
317: In the latter case, it should hold gid[i] < gid[j] iff geid[i] < geid[j], with gid the global orderding for all the dofs
318: and geid the one for the Nedelec field.
320: .seealso: `PCBDDC`, `PCBDDCSetDofsSplitting()`, `PCBDDCSetDofsSplittingLocal()`, `MATAIJ`, `PCBDDCSetDivergenceMat()`
321: @*/
322: PetscErrorCode PCBDDCSetDiscreteGradient(PC pc, Mat G, PetscInt order, PetscInt field, PetscBool global, PetscBool conforming)
323: {
324: PetscFunctionBegin;
331: PetscCheckSameComm(pc, 1, G, 2);
332: PetscTryMethod(pc, "PCBDDCSetDiscreteGradient_C", (PC, Mat, PetscInt, PetscInt, PetscBool, PetscBool), (pc, G, order, field, global, conforming));
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: static PetscErrorCode PCBDDCSetDivergenceMat_BDDC(PC pc, Mat divudotp, PetscBool trans, IS vl2l)
337: {
338: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
340: PetscFunctionBegin;
341: PetscCall(PetscObjectReference((PetscObject)divudotp));
342: PetscCall(MatDestroy(&pcbddc->divudotp));
343: pcbddc->divudotp = divudotp;
344: pcbddc->divudotp_trans = trans;
345: pcbddc->compute_nonetflux = PETSC_TRUE;
346: if (vl2l) {
347: PetscCall(PetscObjectReference((PetscObject)vl2l));
348: PetscCall(ISDestroy(&pcbddc->divudotp_vl2l));
349: pcbddc->divudotp_vl2l = vl2l;
350: }
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: /*@
355: PCBDDCSetDivergenceMat - Sets the linear operator representing \int_\Omega \div {\bf u} \cdot p dx for the `PCBDDC` preconditioner
357: Collective
359: Input Parameters:
360: + pc - the preconditioning context
361: . divudotp - the matrix (must be of type `MATIS`)
362: . trans - if `PETSC_FALSE` (resp. `PETSC_TRUE`), then pressures are in the test (trial) space and velocities are in the trial (test) space.
363: - vl2l - optional index set describing the local (wrt the local matrix in `divudotp`) to local (wrt the local matrix
364: in the preconditioning matrix) map for the velocities
366: Level: advanced
368: Notes:
369: This auxiliary matrix is used to compute quadrature weights representing the net-flux across subdomain boundaries
371: If `vl2l` is `NULL`, the local ordering for velocities in `divudotp` should match that of the preconditioning matrix
373: .seealso: `PCBDDC`, `PCBDDCSetDiscreteGradient()`
374: @*/
375: PetscErrorCode PCBDDCSetDivergenceMat(PC pc, Mat divudotp, PetscBool trans, IS vl2l)
376: {
377: PetscBool ismatis;
379: PetscFunctionBegin;
382: PetscCheckSameComm(pc, 1, divudotp, 2);
385: PetscCall(PetscObjectTypeCompare((PetscObject)divudotp, MATIS, &ismatis));
386: PetscCheck(ismatis, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Divergence matrix needs to be of type MATIS");
387: PetscTryMethod(pc, "PCBDDCSetDivergenceMat_C", (PC, Mat, PetscBool, IS), (pc, divudotp, trans, vl2l));
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: static PetscErrorCode PCBDDCSetChangeOfBasisMat_BDDC(PC pc, Mat change, PetscBool interior)
392: {
393: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
395: PetscFunctionBegin;
396: PetscCall(PetscObjectReference((PetscObject)change));
397: PetscCall(MatDestroy(&pcbddc->user_ChangeOfBasisMatrix));
398: pcbddc->user_ChangeOfBasisMatrix = change;
399: pcbddc->change_interior = interior;
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: /*@
404: PCBDDCSetChangeOfBasisMat - Set user defined change of basis for dofs
406: Collective
408: Input Parameters:
409: + pc - the preconditioning context
410: . change - the change of basis matrix
411: - interior - whether or not the change of basis modifies interior dofs
413: Level: intermediate
415: .seealso: `PCBDDC`
416: @*/
417: PetscErrorCode PCBDDCSetChangeOfBasisMat(PC pc, Mat change, PetscBool interior)
418: {
419: PetscFunctionBegin;
422: PetscCheckSameComm(pc, 1, change, 2);
423: if (pc->mat) {
424: PetscInt rows_c, cols_c, rows, cols;
425: PetscCall(MatGetSize(pc->mat, &rows, &cols));
426: PetscCall(MatGetSize(change, &rows_c, &cols_c));
427: PetscCheck(rows_c == rows, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Invalid number of rows for change of basis matrix! %" PetscInt_FMT " != %" PetscInt_FMT, rows_c, rows);
428: PetscCheck(cols_c == cols, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Invalid number of columns for change of basis matrix! %" PetscInt_FMT " != %" PetscInt_FMT, cols_c, cols);
429: PetscCall(MatGetLocalSize(pc->mat, &rows, &cols));
430: PetscCall(MatGetLocalSize(change, &rows_c, &cols_c));
431: PetscCheck(rows_c == rows, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Invalid number of local rows for change of basis matrix! %" PetscInt_FMT " != %" PetscInt_FMT, rows_c, rows);
432: PetscCheck(cols_c == cols, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Invalid number of local columns for change of basis matrix! %" PetscInt_FMT " != %" PetscInt_FMT, cols_c, cols);
433: }
434: PetscTryMethod(pc, "PCBDDCSetChangeOfBasisMat_C", (PC, Mat, PetscBool), (pc, change, interior));
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: static PetscErrorCode PCBDDCSetPrimalVerticesIS_BDDC(PC pc, IS PrimalVertices)
439: {
440: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
441: PetscBool isequal = PETSC_FALSE;
443: PetscFunctionBegin;
444: PetscCall(PetscObjectReference((PetscObject)PrimalVertices));
445: if (pcbddc->user_primal_vertices) PetscCall(ISEqual(PrimalVertices, pcbddc->user_primal_vertices, &isequal));
446: PetscCall(ISDestroy(&pcbddc->user_primal_vertices));
447: PetscCall(ISDestroy(&pcbddc->user_primal_vertices_local));
448: pcbddc->user_primal_vertices = PrimalVertices;
449: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
450: PetscFunctionReturn(PETSC_SUCCESS);
451: }
453: /*@
454: PCBDDCSetPrimalVerticesIS - Set additional user defined primal vertices in `PCBDDC`
456: Collective
458: Input Parameters:
459: + pc - the preconditioning context
460: - PrimalVertices - index set of primal vertices in global numbering (can be empty)
462: Level: intermediate
464: Note:
465: Any process can list any global node
467: .seealso: `PCBDDC`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
468: @*/
469: PetscErrorCode PCBDDCSetPrimalVerticesIS(PC pc, IS PrimalVertices)
470: {
471: PetscFunctionBegin;
474: PetscCheckSameComm(pc, 1, PrimalVertices, 2);
475: PetscTryMethod(pc, "PCBDDCSetPrimalVerticesIS_C", (PC, IS), (pc, PrimalVertices));
476: PetscFunctionReturn(PETSC_SUCCESS);
477: }
479: static PetscErrorCode PCBDDCGetPrimalVerticesIS_BDDC(PC pc, IS *is)
480: {
481: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
483: PetscFunctionBegin;
484: *is = pcbddc->user_primal_vertices;
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: /*@
489: PCBDDCGetPrimalVerticesIS - Get user defined primal vertices set with `PCBDDCSetPrimalVerticesIS()`
491: Collective
493: Input Parameter:
494: . pc - the preconditioning context
496: Output Parameter:
497: . is - index set of primal vertices in global numbering (`NULL` if not set)
499: Level: intermediate
501: .seealso: `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
502: @*/
503: PetscErrorCode PCBDDCGetPrimalVerticesIS(PC pc, IS *is)
504: {
505: PetscFunctionBegin;
507: PetscAssertPointer(is, 2);
508: PetscUseMethod(pc, "PCBDDCGetPrimalVerticesIS_C", (PC, IS *), (pc, is));
509: PetscFunctionReturn(PETSC_SUCCESS);
510: }
512: static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices)
513: {
514: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
515: PetscBool isequal = PETSC_FALSE;
517: PetscFunctionBegin;
518: PetscCall(PetscObjectReference((PetscObject)PrimalVertices));
519: if (pcbddc->user_primal_vertices_local) PetscCall(ISEqual(PrimalVertices, pcbddc->user_primal_vertices_local, &isequal));
520: PetscCall(ISDestroy(&pcbddc->user_primal_vertices));
521: PetscCall(ISDestroy(&pcbddc->user_primal_vertices_local));
522: pcbddc->user_primal_vertices_local = PrimalVertices;
523: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: /*@
528: PCBDDCSetPrimalVerticesLocalIS - Set additional user defined primal vertices in `PCBDDC`
530: Collective
532: Input Parameters:
533: + pc - the preconditioning context
534: - PrimalVertices - index set of primal vertices in local numbering (can be empty)
536: Level: intermediate
538: .seealso: `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
539: @*/
540: PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices)
541: {
542: PetscFunctionBegin;
545: PetscCheckSameComm(pc, 1, PrimalVertices, 2);
546: PetscTryMethod(pc, "PCBDDCSetPrimalVerticesLocalIS_C", (PC, IS), (pc, PrimalVertices));
547: PetscFunctionReturn(PETSC_SUCCESS);
548: }
550: static PetscErrorCode PCBDDCGetPrimalVerticesLocalIS_BDDC(PC pc, IS *is)
551: {
552: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
554: PetscFunctionBegin;
555: *is = pcbddc->user_primal_vertices_local;
556: PetscFunctionReturn(PETSC_SUCCESS);
557: }
559: /*@
560: PCBDDCGetPrimalVerticesLocalIS - Get user defined primal vertices set with `PCBDDCSetPrimalVerticesLocalIS()`
562: Collective
564: Input Parameter:
565: . pc - the preconditioning context
567: Output Parameter:
568: . is - index set of primal vertices in local numbering (`NULL` if not set)
570: Level: intermediate
572: .seealso: `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`
573: @*/
574: PetscErrorCode PCBDDCGetPrimalVerticesLocalIS(PC pc, IS *is)
575: {
576: PetscFunctionBegin;
578: PetscAssertPointer(is, 2);
579: PetscUseMethod(pc, "PCBDDCGetPrimalVerticesLocalIS_C", (PC, IS *), (pc, is));
580: PetscFunctionReturn(PETSC_SUCCESS);
581: }
583: static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc, PetscInt k)
584: {
585: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
587: PetscFunctionBegin;
588: pcbddc->coarsening_ratio = k;
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: /*@
593: PCBDDCSetCoarseningRatio - Set coarsening ratio used in the multi-level version of `PCBDDC`
595: Logically Collective
597: Input Parameters:
598: + pc - the preconditioning context
599: - k - coarsening ratio (H/h at the coarser level)
601: Options Database Key:
602: . -pc_bddc_coarsening_ratio <int> - Set the coarsening ratio used in multi-level coarsening
604: Level: intermediate
606: Note:
607: Approximately `k` subdomains at the finer level will be aggregated into a single subdomain at the coarser level
609: .seealso: `PCBDDC`, `PCBDDCSetLevels()`
610: @*/
611: PetscErrorCode PCBDDCSetCoarseningRatio(PC pc, PetscInt k)
612: {
613: PetscFunctionBegin;
616: PetscTryMethod(pc, "PCBDDCSetCoarseningRatio_C", (PC, PetscInt), (pc, k));
617: PetscFunctionReturn(PETSC_SUCCESS);
618: }
620: /* The following functions (PCBDDCSetUseExactDirichlet PCBDDCSetLevel) are not public */
621: static PetscErrorCode PCBDDCSetUseExactDirichlet_BDDC(PC pc, PetscBool flg)
622: {
623: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
625: PetscFunctionBegin;
626: pcbddc->use_exact_dirichlet_trick = flg;
627: PetscFunctionReturn(PETSC_SUCCESS);
628: }
630: PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc, PetscBool flg)
631: {
632: PetscFunctionBegin;
635: PetscTryMethod(pc, "PCBDDCSetUseExactDirichlet_C", (PC, PetscBool), (pc, flg));
636: PetscFunctionReturn(PETSC_SUCCESS);
637: }
639: static PetscErrorCode PCBDDCSetLevel_BDDC(PC pc, PetscInt level)
640: {
641: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
643: PetscFunctionBegin;
644: pcbddc->current_level = level;
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
648: PetscErrorCode PCBDDCSetLevel(PC pc, PetscInt level)
649: {
650: PetscFunctionBegin;
653: PetscTryMethod(pc, "PCBDDCSetLevel_C", (PC, PetscInt), (pc, level));
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }
657: static PetscErrorCode PCBDDCSetLevels_BDDC(PC pc, PetscInt levels)
658: {
659: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
661: PetscFunctionBegin;
662: PetscCheck(levels < PETSC_PCBDDC_MAXLEVELS, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Maximum number of additional levels for BDDC is %d", PETSC_PCBDDC_MAXLEVELS - 1);
663: pcbddc->max_levels = levels;
664: PetscFunctionReturn(PETSC_SUCCESS);
665: }
667: /*@
668: PCBDDCSetLevels - Sets the maximum number of additional levels allowed for multilevel `PCBDDC`
670: Logically Collective
672: Input Parameters:
673: + pc - the preconditioning context
674: - levels - the maximum number of levels
676: Options Database Key:
677: . -pc_bddc_levels <int> - Set maximum number of levels for multilevel
679: Level: intermediate
681: Note:
682: The default value is 0, that gives the classical two-levels BDDC algorithm
684: .seealso: `PCBDDC`, `PCBDDCSetCoarseningRatio()`
685: @*/
686: PetscErrorCode PCBDDCSetLevels(PC pc, PetscInt levels)
687: {
688: PetscFunctionBegin;
691: PetscTryMethod(pc, "PCBDDCSetLevels_C", (PC, PetscInt), (pc, levels));
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc, IS DirichletBoundaries)
696: {
697: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
698: PetscBool isequal = PETSC_FALSE;
700: PetscFunctionBegin;
701: PetscCall(PetscObjectReference((PetscObject)DirichletBoundaries));
702: if (pcbddc->DirichletBoundaries) PetscCall(ISEqual(DirichletBoundaries, pcbddc->DirichletBoundaries, &isequal));
703: /* last user setting takes precedence -> destroy any other customization */
704: PetscCall(ISDestroy(&pcbddc->DirichletBoundariesLocal));
705: PetscCall(ISDestroy(&pcbddc->DirichletBoundaries));
706: pcbddc->DirichletBoundaries = DirichletBoundaries;
707: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
708: PetscFunctionReturn(PETSC_SUCCESS);
709: }
711: /*@
712: PCBDDCSetDirichletBoundaries - Set the `IS` defining Dirichlet boundaries for the global problem.
714: Collective
716: Input Parameters:
717: + pc - the preconditioning context
718: - DirichletBoundaries - parallel `IS` defining the Dirichlet boundaries
720: Level: intermediate
722: Note:
723: Provide the information if you used `MatZeroRows()` or `MatZeroRowsColumns()`. Any process can list any global node
725: .seealso: `PCBDDC`, `PCBDDCSetDirichletBoundariesLocal()`, `MatZeroRows()`, `MatZeroRowsColumns()`
726: @*/
727: PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc, IS DirichletBoundaries)
728: {
729: PetscFunctionBegin;
732: PetscCheckSameComm(pc, 1, DirichletBoundaries, 2);
733: PetscTryMethod(pc, "PCBDDCSetDirichletBoundaries_C", (PC, IS), (pc, DirichletBoundaries));
734: PetscFunctionReturn(PETSC_SUCCESS);
735: }
737: static PetscErrorCode PCBDDCSetDirichletBoundariesLocal_BDDC(PC pc, IS DirichletBoundaries)
738: {
739: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
740: PetscBool isequal = PETSC_FALSE;
742: PetscFunctionBegin;
743: PetscCall(PetscObjectReference((PetscObject)DirichletBoundaries));
744: if (pcbddc->DirichletBoundariesLocal) PetscCall(ISEqual(DirichletBoundaries, pcbddc->DirichletBoundariesLocal, &isequal));
745: /* last user setting takes precedence -> destroy any other customization */
746: PetscCall(ISDestroy(&pcbddc->DirichletBoundariesLocal));
747: PetscCall(ISDestroy(&pcbddc->DirichletBoundaries));
748: pcbddc->DirichletBoundariesLocal = DirichletBoundaries;
749: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
750: PetscFunctionReturn(PETSC_SUCCESS);
751: }
753: /*@
754: PCBDDCSetDirichletBoundariesLocal - Set the `IS` defining Dirichlet boundaries for the global problem in local ordering.
756: Collective
758: Input Parameters:
759: + pc - the preconditioning context
760: - DirichletBoundaries - parallel `IS` defining the Dirichlet boundaries (in local ordering)
762: Level: intermediate
764: .seealso: `PCBDDC`, `PCBDDCSetDirichletBoundaries()`, `MatZeroRows()`, `MatZeroRowsColumns()`
765: @*/
766: PetscErrorCode PCBDDCSetDirichletBoundariesLocal(PC pc, IS DirichletBoundaries)
767: {
768: PetscFunctionBegin;
771: PetscCheckSameComm(pc, 1, DirichletBoundaries, 2);
772: PetscTryMethod(pc, "PCBDDCSetDirichletBoundariesLocal_C", (PC, IS), (pc, DirichletBoundaries));
773: PetscFunctionReturn(PETSC_SUCCESS);
774: }
776: static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc, IS NeumannBoundaries)
777: {
778: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
779: PetscBool isequal = PETSC_FALSE;
781: PetscFunctionBegin;
782: PetscCall(PetscObjectReference((PetscObject)NeumannBoundaries));
783: if (pcbddc->NeumannBoundaries) PetscCall(ISEqual(NeumannBoundaries, pcbddc->NeumannBoundaries, &isequal));
784: /* last user setting takes precedence -> destroy any other customization */
785: PetscCall(ISDestroy(&pcbddc->NeumannBoundariesLocal));
786: PetscCall(ISDestroy(&pcbddc->NeumannBoundaries));
787: pcbddc->NeumannBoundaries = NeumannBoundaries;
788: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
789: PetscFunctionReturn(PETSC_SUCCESS);
790: }
792: /*@
793: PCBDDCSetNeumannBoundaries - Set the `IS` defining Neumann boundaries for the global problem.
795: Collective
797: Input Parameters:
798: + pc - the preconditioning context
799: - NeumannBoundaries - parallel `IS` defining the Neumann boundaries
801: Level: intermediate
803: Note:
804: Any process can list any global node
806: .seealso: `PCBDDC`, `PCBDDCSetNeumannBoundariesLocal()`
807: @*/
808: PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc, IS NeumannBoundaries)
809: {
810: PetscFunctionBegin;
813: PetscCheckSameComm(pc, 1, NeumannBoundaries, 2);
814: PetscTryMethod(pc, "PCBDDCSetNeumannBoundaries_C", (PC, IS), (pc, NeumannBoundaries));
815: PetscFunctionReturn(PETSC_SUCCESS);
816: }
818: static PetscErrorCode PCBDDCSetNeumannBoundariesLocal_BDDC(PC pc, IS NeumannBoundaries)
819: {
820: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
821: PetscBool isequal = PETSC_FALSE;
823: PetscFunctionBegin;
824: PetscCall(PetscObjectReference((PetscObject)NeumannBoundaries));
825: if (pcbddc->NeumannBoundariesLocal) PetscCall(ISEqual(NeumannBoundaries, pcbddc->NeumannBoundariesLocal, &isequal));
826: /* last user setting takes precedence -> destroy any other customization */
827: PetscCall(ISDestroy(&pcbddc->NeumannBoundariesLocal));
828: PetscCall(ISDestroy(&pcbddc->NeumannBoundaries));
829: pcbddc->NeumannBoundariesLocal = NeumannBoundaries;
830: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
831: PetscFunctionReturn(PETSC_SUCCESS);
832: }
834: /*@
835: PCBDDCSetNeumannBoundariesLocal - Set the `IS` defining Neumann boundaries for the global problem in local ordering.
837: Collective
839: Input Parameters:
840: + pc - the preconditioning context
841: - NeumannBoundaries - parallel `IS` defining the subdomain part of Neumann boundaries (in local ordering)
843: Level: intermediate
845: .seealso: `PCBDDC`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCGetDirichletBoundaries()`
846: @*/
847: PetscErrorCode PCBDDCSetNeumannBoundariesLocal(PC pc, IS NeumannBoundaries)
848: {
849: PetscFunctionBegin;
852: PetscCheckSameComm(pc, 1, NeumannBoundaries, 2);
853: PetscTryMethod(pc, "PCBDDCSetNeumannBoundariesLocal_C", (PC, IS), (pc, NeumannBoundaries));
854: PetscFunctionReturn(PETSC_SUCCESS);
855: }
857: static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc, IS *DirichletBoundaries)
858: {
859: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
861: PetscFunctionBegin;
862: *DirichletBoundaries = pcbddc->DirichletBoundaries;
863: PetscFunctionReturn(PETSC_SUCCESS);
864: }
866: /*@
867: PCBDDCGetDirichletBoundaries - Get parallel `IS` for Dirichlet boundaries
869: Collective
871: Input Parameter:
872: . pc - the preconditioning context
874: Output Parameter:
875: . DirichletBoundaries - index set defining the Dirichlet boundaries
877: Level: intermediate
879: Note:
880: The `IS` returned (if any) is the same passed in earlier by the user with `PCBDDCSetDirichletBoundaries()`
882: .seealso: `PCBDDC`, `PCBDDCSetDirichletBoundaries()`
883: @*/
884: PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc, IS *DirichletBoundaries)
885: {
886: PetscFunctionBegin;
888: PetscUseMethod(pc, "PCBDDCGetDirichletBoundaries_C", (PC, IS *), (pc, DirichletBoundaries));
889: PetscFunctionReturn(PETSC_SUCCESS);
890: }
892: static PetscErrorCode PCBDDCGetDirichletBoundariesLocal_BDDC(PC pc, IS *DirichletBoundaries)
893: {
894: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
896: PetscFunctionBegin;
897: *DirichletBoundaries = pcbddc->DirichletBoundariesLocal;
898: PetscFunctionReturn(PETSC_SUCCESS);
899: }
901: /*@
902: PCBDDCGetDirichletBoundariesLocal - Get parallel `IS` for Dirichlet boundaries (in local ordering)
904: Collective
906: Input Parameter:
907: . pc - the preconditioning context
909: Output Parameter:
910: . DirichletBoundaries - index set defining the subdomain part of Dirichlet boundaries
912: Level: intermediate
914: Note:
915: The `IS` returned could be the same passed in earlier by the user (if provided with `PCBDDCSetDirichletBoundariesLocal()`)
916: or a global-to-local map of the global `IS` (if provided with `PCBDDCSetDirichletBoundaries()`).
917: In the latter case, the `IS` will be available only after `PCSetUp()`.
919: .seealso: `PCBDDC`, `PCBDDCGetDirichletBoundaries()`, `PCBDDCSetDirichletBoundaries()`
920: @*/
921: PetscErrorCode PCBDDCGetDirichletBoundariesLocal(PC pc, IS *DirichletBoundaries)
922: {
923: PetscFunctionBegin;
925: PetscUseMethod(pc, "PCBDDCGetDirichletBoundariesLocal_C", (PC, IS *), (pc, DirichletBoundaries));
926: PetscFunctionReturn(PETSC_SUCCESS);
927: }
929: static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc, IS *NeumannBoundaries)
930: {
931: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
933: PetscFunctionBegin;
934: *NeumannBoundaries = pcbddc->NeumannBoundaries;
935: PetscFunctionReturn(PETSC_SUCCESS);
936: }
938: /*@
939: PCBDDCGetNeumannBoundaries - Get parallel `IS` for Neumann boundaries
941: Not Collective
943: Input Parameter:
944: . pc - the preconditioning context
946: Output Parameter:
947: . NeumannBoundaries - index set defining the Neumann boundaries
949: Level: intermediate
951: Note:
952: The `IS` returned (if any) is the same passed in earlier by the user with `PCBDDCSetNeumannBoundaries()`
954: .seealso: `PCBDDC`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCGetDirichletBoundaries()`, `PCBDDCSetDirichletBoundaries()`
955: @*/
956: PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc, IS *NeumannBoundaries)
957: {
958: PetscFunctionBegin;
960: PetscUseMethod(pc, "PCBDDCGetNeumannBoundaries_C", (PC, IS *), (pc, NeumannBoundaries));
961: PetscFunctionReturn(PETSC_SUCCESS);
962: }
964: static PetscErrorCode PCBDDCGetNeumannBoundariesLocal_BDDC(PC pc, IS *NeumannBoundaries)
965: {
966: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
968: PetscFunctionBegin;
969: *NeumannBoundaries = pcbddc->NeumannBoundariesLocal;
970: PetscFunctionReturn(PETSC_SUCCESS);
971: }
973: /*@
974: PCBDDCGetNeumannBoundariesLocal - Get parallel `IS` for Neumann boundaries (in local ordering)
976: Not Collective
978: Input Parameter:
979: . pc - the preconditioning context
981: Output Parameter:
982: . NeumannBoundaries - index set defining the subdomain part of Neumann boundaries
984: Level: intermediate
986: Note:
987: The `IS` returned could be the same passed in earlier by the user (if provided with `PCBDDCSetNeumannBoundariesLocal()`
988: or a global-to-local map of the global `IS` (if provided with `PCBDDCSetNeumannBoundaries()`).
989: In the latter case, the `IS` will be available after `PCSetUp()`.
991: .seealso: `PCBDDC`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCSetNeumannBoundariesLocal)`, `PCBDDCGetNeumannBoundaries()`
992: @*/
993: PetscErrorCode PCBDDCGetNeumannBoundariesLocal(PC pc, IS *NeumannBoundaries)
994: {
995: PetscFunctionBegin;
997: PetscUseMethod(pc, "PCBDDCGetNeumannBoundariesLocal_C", (PC, IS *), (pc, NeumannBoundaries));
998: PetscFunctionReturn(PETSC_SUCCESS);
999: }
1001: static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs, const PetscInt xadj[], const PetscInt adjncy[], PetscCopyMode copymode)
1002: {
1003: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1004: PCBDDCGraph mat_graph = pcbddc->mat_graph;
1005: PetscBool same_data = PETSC_FALSE;
1007: PetscFunctionBegin;
1008: if (!nvtxs) {
1009: if (copymode == PETSC_OWN_POINTER) {
1010: PetscCall(PetscFree(xadj));
1011: PetscCall(PetscFree(adjncy));
1012: }
1013: PetscCall(PCBDDCGraphResetCSR(mat_graph));
1014: PetscFunctionReturn(PETSC_SUCCESS);
1015: }
1016: if (mat_graph->nvtxs == nvtxs && mat_graph->freecsr) { /* we own the data */
1017: if (mat_graph->xadj == xadj && mat_graph->adjncy == adjncy) same_data = PETSC_TRUE;
1018: if (!same_data && mat_graph->xadj[nvtxs] == xadj[nvtxs]) {
1019: PetscCall(PetscArraycmp(xadj, mat_graph->xadj, nvtxs + 1, &same_data));
1020: if (same_data) PetscCall(PetscArraycmp(adjncy, mat_graph->adjncy, xadj[nvtxs], &same_data));
1021: }
1022: }
1023: if (!same_data) {
1024: /* free old CSR */
1025: PetscCall(PCBDDCGraphResetCSR(mat_graph));
1026: /* get CSR into graph structure */
1027: if (copymode == PETSC_COPY_VALUES) {
1028: PetscCall(PetscMalloc1(nvtxs + 1, &mat_graph->xadj));
1029: PetscCall(PetscMalloc1(xadj[nvtxs], &mat_graph->adjncy));
1030: PetscCall(PetscArraycpy(mat_graph->xadj, xadj, nvtxs + 1));
1031: PetscCall(PetscArraycpy(mat_graph->adjncy, adjncy, xadj[nvtxs]));
1032: mat_graph->freecsr = PETSC_TRUE;
1033: } else if (copymode == PETSC_OWN_POINTER) {
1034: mat_graph->xadj = (PetscInt *)xadj;
1035: mat_graph->adjncy = (PetscInt *)adjncy;
1036: mat_graph->freecsr = PETSC_TRUE;
1037: } else if (copymode == PETSC_USE_POINTER) {
1038: mat_graph->xadj = (PetscInt *)xadj;
1039: mat_graph->adjncy = (PetscInt *)adjncy;
1040: mat_graph->freecsr = PETSC_FALSE;
1041: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported copy mode %d", copymode);
1042: mat_graph->nvtxs_csr = nvtxs;
1043: pcbddc->recompute_topography = PETSC_TRUE;
1044: }
1045: PetscFunctionReturn(PETSC_SUCCESS);
1046: }
1048: /*@
1049: PCBDDCSetLocalAdjacencyGraph - Set adjacency structure (CSR graph) of the local degrees of freedom.
1051: Not collective
1053: Input Parameters:
1054: + pc - the preconditioning context.
1055: . nvtxs - number of local vertices of the graph (i.e., the number of local dofs).
1056: . xadj - CSR format row pointers for the connectivity of the dofs
1057: . adjncy - CSR format column pointers for the connectivity of the dofs
1058: - copymode - supported modes are `PETSC_COPY_VALUES`, `PETSC_USE_POINTER` or `PETSC_OWN_POINTER`.
1060: Level: intermediate
1062: Note:
1063: A dof is considered connected with all local dofs if xadj[dof+1]-xadj[dof] == 1 and adjncy[xadj[dof]] is negative.
1065: .seealso: `PCBDDC`, `PetscCopyMode`
1066: @*/
1067: PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc, PetscInt nvtxs, const PetscInt xadj[], const PetscInt adjncy[], PetscCopyMode copymode)
1068: {
1069: void (*f)(void) = NULL;
1071: PetscFunctionBegin;
1073: if (nvtxs) {
1074: PetscAssertPointer(xadj, 3);
1075: if (xadj[nvtxs]) PetscAssertPointer(adjncy, 4);
1076: }
1077: PetscTryMethod(pc, "PCBDDCSetLocalAdjacencyGraph_C", (PC, PetscInt, const PetscInt[], const PetscInt[], PetscCopyMode), (pc, nvtxs, xadj, adjncy, copymode));
1078: /* free arrays if PCBDDC is not the PC type */
1079: PetscCall(PetscObjectQueryFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", &f));
1080: if (!f && copymode == PETSC_OWN_POINTER) {
1081: PetscCall(PetscFree(xadj));
1082: PetscCall(PetscFree(adjncy));
1083: }
1084: PetscFunctionReturn(PETSC_SUCCESS);
1085: }
1087: static PetscErrorCode PCBDDCSetDofsSplittingLocal_BDDC(PC pc, PetscInt n_is, IS ISForDofs[])
1088: {
1089: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1090: PetscInt i;
1091: PetscBool isequal = PETSC_FALSE;
1093: PetscFunctionBegin;
1094: if (pcbddc->n_ISForDofsLocal == n_is) {
1095: for (i = 0; i < n_is; i++) {
1096: PetscBool isequalt;
1097: PetscCall(ISEqual(ISForDofs[i], pcbddc->ISForDofsLocal[i], &isequalt));
1098: if (!isequalt) break;
1099: }
1100: if (i == n_is) isequal = PETSC_TRUE;
1101: }
1102: for (i = 0; i < n_is; i++) PetscCall(PetscObjectReference((PetscObject)ISForDofs[i]));
1103: /* Destroy ISes if they were already set */
1104: for (i = 0; i < pcbddc->n_ISForDofsLocal; i++) PetscCall(ISDestroy(&pcbddc->ISForDofsLocal[i]));
1105: PetscCall(PetscFree(pcbddc->ISForDofsLocal));
1106: /* last user setting takes precedence -> destroy any other customization */
1107: for (i = 0; i < pcbddc->n_ISForDofs; i++) PetscCall(ISDestroy(&pcbddc->ISForDofs[i]));
1108: PetscCall(PetscFree(pcbddc->ISForDofs));
1109: pcbddc->n_ISForDofs = 0;
1110: /* allocate space then set */
1111: if (n_is) PetscCall(PetscMalloc1(n_is, &pcbddc->ISForDofsLocal));
1112: for (i = 0; i < n_is; i++) pcbddc->ISForDofsLocal[i] = ISForDofs[i];
1113: pcbddc->n_ISForDofsLocal = n_is;
1114: if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
1115: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1116: PetscFunctionReturn(PETSC_SUCCESS);
1117: }
1119: /*@
1120: PCBDDCSetDofsSplittingLocal - Set the `IS` defining fields of the local subdomain matrix
1122: Collective
1124: Input Parameters:
1125: + pc - the preconditioning context
1126: . n_is - number of index sets defining the fields, must be the same on all MPI processes
1127: - ISForDofs - array of `IS` describing the fields in local ordering
1129: Level: intermediate
1131: Note:
1132: Not all nodes need to be listed, unlisted nodes will belong to the complement field.
1134: .seealso: `PCBDDC`, `PCBDDCSetDofsSplitting()`
1135: @*/
1136: PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc, PetscInt n_is, IS ISForDofs[])
1137: {
1138: PetscInt i;
1140: PetscFunctionBegin;
1143: for (i = 0; i < n_is; i++) {
1144: PetscCheckSameComm(pc, 1, ISForDofs[i], 3);
1146: }
1147: PetscTryMethod(pc, "PCBDDCSetDofsSplittingLocal_C", (PC, PetscInt, IS[]), (pc, n_is, ISForDofs));
1148: PetscFunctionReturn(PETSC_SUCCESS);
1149: }
1151: static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc, PetscInt n_is, IS ISForDofs[])
1152: {
1153: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1154: PetscInt i;
1155: PetscBool isequal = PETSC_FALSE;
1157: PetscFunctionBegin;
1158: if (pcbddc->n_ISForDofs == n_is) {
1159: for (i = 0; i < n_is; i++) {
1160: PetscBool isequalt;
1161: PetscCall(ISEqual(ISForDofs[i], pcbddc->ISForDofs[i], &isequalt));
1162: if (!isequalt) break;
1163: }
1164: if (i == n_is) isequal = PETSC_TRUE;
1165: }
1166: for (i = 0; i < n_is; i++) PetscCall(PetscObjectReference((PetscObject)ISForDofs[i]));
1167: /* Destroy ISes if they were already set */
1168: for (i = 0; i < pcbddc->n_ISForDofs; i++) PetscCall(ISDestroy(&pcbddc->ISForDofs[i]));
1169: PetscCall(PetscFree(pcbddc->ISForDofs));
1170: /* last user setting takes precedence -> destroy any other customization */
1171: for (i = 0; i < pcbddc->n_ISForDofsLocal; i++) PetscCall(ISDestroy(&pcbddc->ISForDofsLocal[i]));
1172: PetscCall(PetscFree(pcbddc->ISForDofsLocal));
1173: pcbddc->n_ISForDofsLocal = 0;
1174: /* allocate space then set */
1175: if (n_is) PetscCall(PetscMalloc1(n_is, &pcbddc->ISForDofs));
1176: for (i = 0; i < n_is; i++) pcbddc->ISForDofs[i] = ISForDofs[i];
1177: pcbddc->n_ISForDofs = n_is;
1178: if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
1179: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1180: PetscFunctionReturn(PETSC_SUCCESS);
1181: }
1183: /*@
1184: PCBDDCSetDofsSplitting - Set the `IS` defining fields of the global matrix
1186: Collective
1188: Input Parameters:
1189: + pc - the preconditioning context
1190: . n_is - number of index sets defining the fields
1191: - ISForDofs - array of `IS` describing the fields in global ordering
1193: Level: intermediate
1195: Note:
1196: Any process can list any global node. Not all nodes need to be listed: unlisted nodes will belong to the complement field.
1198: .seealso: `PCBDDC`, `PCBDDCSetDofsSplittingLocal()`
1199: @*/
1200: PetscErrorCode PCBDDCSetDofsSplitting(PC pc, PetscInt n_is, IS ISForDofs[])
1201: {
1202: PetscInt i;
1204: PetscFunctionBegin;
1207: for (i = 0; i < n_is; i++) {
1209: PetscCheckSameComm(pc, 1, ISForDofs[i], 3);
1210: }
1211: PetscTryMethod(pc, "PCBDDCSetDofsSplitting_C", (PC, PetscInt, IS[]), (pc, n_is, ISForDofs));
1212: PetscFunctionReturn(PETSC_SUCCESS);
1213: }
1215: /*
1216: PCPreSolve_BDDC - Changes the right hand side and (if necessary) the initial
1217: guess if a transformation of basis approach has been selected.
1219: Input Parameter:
1220: + pc - the preconditioner context
1222: Note:
1223: The interface routine PCPreSolve() is not usually called directly by
1224: the user, but instead is called by KSPSolve().
1225: */
1226: static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
1227: {
1228: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1229: PC_IS *pcis = (PC_IS *)(pc->data);
1230: Vec used_vec;
1231: PetscBool iscg, save_rhs = PETSC_TRUE, benign_correction_computed;
1233: PetscFunctionBegin;
1234: /* if we are working with CG, one dirichlet solve can be avoided during Krylov iterations */
1235: if (ksp) {
1236: PetscCall(PetscObjectTypeCompareAny((PetscObject)ksp, &iscg, KSPCG, KSPGROPPCG, KSPPIPECG, KSPPIPELCG, KSPPIPECGRR, ""));
1237: if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || !iscg || pc->mat != pc->pmat) PetscCall(PCBDDCSetUseExactDirichlet(pc, PETSC_FALSE));
1238: }
1239: if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || pc->mat != pc->pmat) PetscCall(PCBDDCSetUseExactDirichlet(pc, PETSC_FALSE));
1241: /* Creates parallel work vectors used in presolve */
1242: if (!pcbddc->original_rhs) PetscCall(VecDuplicate(pcis->vec1_global, &pcbddc->original_rhs));
1243: if (!pcbddc->temp_solution) PetscCall(VecDuplicate(pcis->vec1_global, &pcbddc->temp_solution));
1245: pcbddc->temp_solution_used = PETSC_FALSE;
1246: if (x) {
1247: PetscCall(PetscObjectReference((PetscObject)x));
1248: used_vec = x;
1249: } else { /* it can only happen when calling PCBDDCMatFETIDPGetRHS */
1250: PetscCall(PetscObjectReference((PetscObject)pcbddc->temp_solution));
1251: used_vec = pcbddc->temp_solution;
1252: PetscCall(VecSet(used_vec, 0.0));
1253: pcbddc->temp_solution_used = PETSC_TRUE;
1254: PetscCall(VecCopy(rhs, pcbddc->original_rhs));
1255: save_rhs = PETSC_FALSE;
1256: pcbddc->eliminate_dirdofs = PETSC_TRUE;
1257: }
1259: /* hack into ksp data structure since PCPreSolve comes earlier than setting to zero the guess in src/ksp/ksp/interface/itfunc.c */
1260: if (ksp) {
1261: /* store the flag for the initial guess since it will be restored back during PCPostSolve_BDDC */
1262: PetscCall(KSPGetInitialGuessNonzero(ksp, &pcbddc->ksp_guess_nonzero));
1263: if (!pcbddc->ksp_guess_nonzero) PetscCall(VecSet(used_vec, 0.0));
1264: }
1266: pcbddc->rhs_change = PETSC_FALSE;
1267: /* Take into account zeroed rows -> change rhs and store solution removed */
1268: if (rhs && pcbddc->eliminate_dirdofs) {
1269: IS dirIS = NULL;
1271: /* DirichletBoundariesLocal may not be consistent among neighbours; gets a dirichlet dofs IS from graph (may be cached) */
1272: PetscCall(PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph, &dirIS));
1273: if (dirIS) {
1274: Mat_IS *matis = (Mat_IS *)pc->pmat->data;
1275: PetscInt dirsize, i, *is_indices;
1276: PetscScalar *array_x;
1277: const PetscScalar *array_diagonal;
1279: PetscCall(MatGetDiagonal(pc->pmat, pcis->vec1_global));
1280: PetscCall(VecPointwiseDivide(pcis->vec1_global, rhs, pcis->vec1_global));
1281: PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_global, pcis->vec2_N, INSERT_VALUES, SCATTER_FORWARD));
1282: PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_global, pcis->vec2_N, INSERT_VALUES, SCATTER_FORWARD));
1283: PetscCall(VecScatterBegin(matis->rctx, used_vec, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD));
1284: PetscCall(VecScatterEnd(matis->rctx, used_vec, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD));
1285: PetscCall(ISGetLocalSize(dirIS, &dirsize));
1286: PetscCall(VecGetArray(pcis->vec1_N, &array_x));
1287: PetscCall(VecGetArrayRead(pcis->vec2_N, &array_diagonal));
1288: PetscCall(ISGetIndices(dirIS, (const PetscInt **)&is_indices));
1289: for (i = 0; i < dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]];
1290: PetscCall(ISRestoreIndices(dirIS, (const PetscInt **)&is_indices));
1291: PetscCall(VecRestoreArrayRead(pcis->vec2_N, &array_diagonal));
1292: PetscCall(VecRestoreArray(pcis->vec1_N, &array_x));
1293: PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_N, used_vec, INSERT_VALUES, SCATTER_REVERSE));
1294: PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_N, used_vec, INSERT_VALUES, SCATTER_REVERSE));
1295: pcbddc->rhs_change = PETSC_TRUE;
1296: PetscCall(ISDestroy(&dirIS));
1297: }
1298: }
1300: /* remove the computed solution or the initial guess from the rhs */
1301: if (pcbddc->rhs_change || (ksp && pcbddc->ksp_guess_nonzero)) {
1302: /* save the original rhs */
1303: if (save_rhs) {
1304: PetscCall(VecSwap(rhs, pcbddc->original_rhs));
1305: save_rhs = PETSC_FALSE;
1306: }
1307: pcbddc->rhs_change = PETSC_TRUE;
1308: PetscCall(VecScale(used_vec, -1.0));
1309: PetscCall(MatMultAdd(pc->mat, used_vec, pcbddc->original_rhs, rhs));
1310: PetscCall(VecScale(used_vec, -1.0));
1311: PetscCall(VecCopy(used_vec, pcbddc->temp_solution));
1312: pcbddc->temp_solution_used = PETSC_TRUE;
1313: if (ksp) PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_FALSE));
1314: }
1315: PetscCall(VecDestroy(&used_vec));
1317: /* compute initial vector in benign space if needed
1318: and remove non-benign solution from the rhs */
1319: benign_correction_computed = PETSC_FALSE;
1320: if (rhs && pcbddc->benign_compute_correction && (pcbddc->benign_have_null || pcbddc->benign_apply_coarse_only)) {
1321: /* compute u^*_h using ideas similar to those in Xuemin Tu's PhD thesis (see Section 4.8.1)
1322: Recursively apply BDDC in the multilevel case */
1323: if (!pcbddc->benign_vec) PetscCall(VecDuplicate(rhs, &pcbddc->benign_vec));
1324: /* keep applying coarse solver unless we no longer have benign subdomains */
1325: pcbddc->benign_apply_coarse_only = pcbddc->benign_have_null ? PETSC_TRUE : PETSC_FALSE;
1326: if (!pcbddc->benign_skip_correction) {
1327: PetscCall(PCApply_BDDC(pc, rhs, pcbddc->benign_vec));
1328: benign_correction_computed = PETSC_TRUE;
1329: if (pcbddc->temp_solution_used) PetscCall(VecAXPY(pcbddc->temp_solution, 1.0, pcbddc->benign_vec));
1330: PetscCall(VecScale(pcbddc->benign_vec, -1.0));
1331: /* store the original rhs if not done earlier */
1332: if (save_rhs) PetscCall(VecSwap(rhs, pcbddc->original_rhs));
1333: if (pcbddc->rhs_change) {
1334: PetscCall(MatMultAdd(pc->mat, pcbddc->benign_vec, rhs, rhs));
1335: } else {
1336: PetscCall(MatMultAdd(pc->mat, pcbddc->benign_vec, pcbddc->original_rhs, rhs));
1337: }
1338: pcbddc->rhs_change = PETSC_TRUE;
1339: }
1340: pcbddc->benign_apply_coarse_only = PETSC_FALSE;
1341: } else {
1342: PetscCall(VecDestroy(&pcbddc->benign_vec));
1343: }
1345: /* dbg output */
1346: if (pcbddc->dbg_flag && benign_correction_computed) {
1347: Vec v;
1349: PetscCall(VecDuplicate(pcis->vec1_global, &v));
1350: if (pcbddc->ChangeOfBasisMatrix) {
1351: PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, rhs, v));
1352: } else {
1353: PetscCall(VecCopy(rhs, v));
1354: }
1355: PetscCall(PCBDDCBenignGetOrSetP0(pc, v, PETSC_TRUE));
1356: PetscCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "LEVEL %" PetscInt_FMT ": is the correction benign?\n", pcbddc->current_level));
1357: PetscCall(PetscScalarView(pcbddc->benign_n, pcbddc->benign_p0, pcbddc->dbg_viewer));
1358: PetscCall(PetscViewerFlush(pcbddc->dbg_viewer));
1359: PetscCall(VecDestroy(&v));
1360: }
1362: /* set initial guess if using PCG */
1363: pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
1364: if (x && pcbddc->use_exact_dirichlet_trick) {
1365: PetscCall(VecSet(x, 0.0));
1366: if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) {
1367: if (benign_correction_computed) { /* we have already saved the changed rhs */
1368: PetscCall(VecLockReadPop(pcis->vec1_global));
1369: } else {
1370: PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, rhs, pcis->vec1_global));
1371: }
1372: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec1_global, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1373: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec1_global, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1374: } else {
1375: PetscCall(VecScatterBegin(pcis->global_to_D, rhs, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1376: PetscCall(VecScatterEnd(pcis->global_to_D, rhs, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1377: }
1378: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1379: PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D));
1380: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1381: PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D));
1382: if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) {
1383: PetscCall(VecSet(pcis->vec1_global, 0.));
1384: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, pcis->vec1_global, INSERT_VALUES, SCATTER_REVERSE));
1385: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, pcis->vec1_global, INSERT_VALUES, SCATTER_REVERSE));
1386: PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcis->vec1_global, x));
1387: } else {
1388: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, x, INSERT_VALUES, SCATTER_REVERSE));
1389: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, x, INSERT_VALUES, SCATTER_REVERSE));
1390: }
1391: if (ksp) PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
1392: pcbddc->exact_dirichlet_trick_app = PETSC_TRUE;
1393: } else if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior && benign_correction_computed && pcbddc->use_exact_dirichlet_trick) {
1394: PetscCall(VecLockReadPop(pcis->vec1_global));
1395: }
1396: PetscFunctionReturn(PETSC_SUCCESS);
1397: }
1399: /*
1400: PCPostSolve_BDDC - Changes the computed solution if a transformation of basis
1401: approach has been selected. Also, restores rhs to its original state.
1403: Input Parameter:
1404: + pc - the preconditioner context
1406: Application Interface Routine: PCPostSolve()
1408: Note:
1409: The interface routine PCPostSolve() is not usually called directly by
1410: the user, but instead is called by KSPSolve().
1411: */
1412: static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
1413: {
1414: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1416: PetscFunctionBegin;
1417: /* add solution removed in presolve */
1418: if (x && pcbddc->rhs_change) {
1419: if (pcbddc->temp_solution_used) {
1420: PetscCall(VecAXPY(x, 1.0, pcbddc->temp_solution));
1421: } else if (pcbddc->benign_compute_correction && pcbddc->benign_vec) {
1422: PetscCall(VecAXPY(x, -1.0, pcbddc->benign_vec));
1423: }
1424: /* restore to original state (not for FETI-DP) */
1425: if (ksp) pcbddc->temp_solution_used = PETSC_FALSE;
1426: }
1428: /* restore rhs to its original state (not needed for FETI-DP) */
1429: if (rhs && pcbddc->rhs_change) {
1430: PetscCall(VecSwap(rhs, pcbddc->original_rhs));
1431: pcbddc->rhs_change = PETSC_FALSE;
1432: }
1433: /* restore ksp guess state */
1434: if (ksp) {
1435: PetscCall(KSPSetInitialGuessNonzero(ksp, pcbddc->ksp_guess_nonzero));
1436: /* reset flag for exact dirichlet trick */
1437: pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
1438: }
1439: PetscFunctionReturn(PETSC_SUCCESS);
1440: }
1442: // PetscClangLinter pragma disable: -fdoc-sowing-chars
1443: /*
1444: PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner
1445: by setting data structures and options.
1447: Input Parameter:
1448: . pc - the preconditioner context
1450: Application Interface Routine: PCSetUp()
1452: */
1453: static PetscErrorCode PCSetUp_BDDC(PC pc)
1454: {
1455: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1456: PCBDDCSubSchurs sub_schurs;
1457: Mat_IS *matis;
1458: MatNullSpace nearnullspace;
1459: Mat lA;
1460: IS lP, zerodiag = NULL;
1461: PetscInt nrows, ncols;
1462: PetscMPIInt size;
1463: PetscBool computesubschurs;
1464: PetscBool computeconstraintsmatrix;
1465: PetscBool new_nearnullspace_provided, ismatis, rl;
1466: PetscBool isset, issym, isspd;
1468: PetscFunctionBegin;
1469: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATIS, &ismatis));
1470: PetscCheck(ismatis, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PCBDDC preconditioner requires matrix of type MATIS");
1471: PetscCall(MatGetSize(pc->pmat, &nrows, &ncols));
1472: PetscCheck(nrows == ncols, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCBDDC preconditioner requires a square preconditioning matrix");
1473: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
1475: matis = (Mat_IS *)pc->pmat->data;
1476: /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */
1477: /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup
1478: Also, BDDC builds its own KSP for the Dirichlet problem */
1479: rl = pcbddc->recompute_topography;
1480: if (!pc->setupcalled || pc->flag == DIFFERENT_NONZERO_PATTERN) rl = PETSC_TRUE;
1481: PetscCall(MPIU_Allreduce(&rl, &pcbddc->recompute_topography, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)pc)));
1482: if (pcbddc->recompute_topography) {
1483: pcbddc->graphanalyzed = PETSC_FALSE;
1484: computeconstraintsmatrix = PETSC_TRUE;
1485: } else {
1486: computeconstraintsmatrix = PETSC_FALSE;
1487: }
1489: /* check parameters' compatibility */
1490: if (!pcbddc->use_deluxe_scaling) pcbddc->deluxe_zerorows = PETSC_FALSE;
1491: pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_threshold[0] != 0.0 || pcbddc->adaptive_threshold[1] != 0.0);
1492: pcbddc->use_deluxe_scaling = (PetscBool)(pcbddc->use_deluxe_scaling && size > 1);
1493: pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_selection && size > 1);
1494: pcbddc->adaptive_userdefined = (PetscBool)(pcbddc->adaptive_selection && pcbddc->adaptive_userdefined);
1495: if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE;
1497: computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling);
1499: /* activate all connected components if the netflux has been requested */
1500: if (pcbddc->compute_nonetflux) {
1501: pcbddc->use_vertices = PETSC_TRUE;
1502: pcbddc->use_edges = PETSC_TRUE;
1503: pcbddc->use_faces = PETSC_TRUE;
1504: }
1506: /* Get stdout for dbg */
1507: if (pcbddc->dbg_flag) {
1508: if (!pcbddc->dbg_viewer) pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc));
1509: PetscCall(PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer));
1510: PetscCall(PetscViewerASCIIAddTab(pcbddc->dbg_viewer, 2 * pcbddc->current_level));
1511: }
1513: /* process topology information */
1514: PetscCall(PetscLogEventBegin(PC_BDDC_Topology[pcbddc->current_level], pc, 0, 0, 0));
1515: if (pcbddc->recompute_topography) {
1516: PetscCall(PCBDDCComputeLocalTopologyInfo(pc));
1517: if (pcbddc->discretegradient) PetscCall(PCBDDCNedelecSupport(pc));
1518: }
1519: if (pcbddc->corner_selected) pcbddc->use_vertices = PETSC_TRUE;
1521: /* change basis if requested by the user */
1522: if (pcbddc->user_ChangeOfBasisMatrix) {
1523: /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */
1524: pcbddc->use_change_of_basis = PETSC_FALSE;
1525: PetscCall(PCBDDCComputeLocalMatrix(pc, pcbddc->user_ChangeOfBasisMatrix));
1526: } else {
1527: PetscCall(MatDestroy(&pcbddc->local_mat));
1528: PetscCall(PetscObjectReference((PetscObject)matis->A));
1529: pcbddc->local_mat = matis->A;
1530: }
1532: /*
1533: Compute change of basis on local pressures (aka zerodiag dofs) with the benign trick
1534: This should come earlier then PCISSetUp for extracting the correct subdomain matrices
1535: */
1536: PetscCall(PCBDDCBenignShellMat(pc, PETSC_TRUE));
1537: if (pcbddc->benign_saddle_point) {
1538: PC_IS *pcis = (PC_IS *)pc->data;
1540: if (pcbddc->user_ChangeOfBasisMatrix || pcbddc->use_change_of_basis || !computesubschurs) pcbddc->benign_change_explicit = PETSC_TRUE;
1541: /* detect local saddle point and change the basis in pcbddc->local_mat */
1542: PetscCall(PCBDDCBenignDetectSaddlePoint(pc, (PetscBool)(!pcbddc->recompute_topography), &zerodiag));
1543: /* pop B0 mat from local mat */
1544: PetscCall(PCBDDCBenignPopOrPushB0(pc, PETSC_TRUE));
1545: /* give pcis a hint to not reuse submatrices during PCISCreate */
1546: if (pc->flag == SAME_NONZERO_PATTERN && pcis->reusesubmatrices == PETSC_TRUE) {
1547: if (pcbddc->benign_n && (pcbddc->benign_change_explicit || pcbddc->dbg_flag)) {
1548: pcis->reusesubmatrices = PETSC_FALSE;
1549: } else {
1550: pcis->reusesubmatrices = PETSC_TRUE;
1551: }
1552: } else {
1553: pcis->reusesubmatrices = PETSC_FALSE;
1554: }
1555: }
1557: /* propagate relevant information */
1558: PetscCall(MatIsSymmetricKnown(matis->A, &isset, &issym));
1559: if (isset) PetscCall(MatSetOption(pcbddc->local_mat, MAT_SYMMETRIC, issym));
1560: PetscCall(MatIsSPDKnown(matis->A, &isset, &isspd));
1561: if (isset) PetscCall(MatSetOption(pcbddc->local_mat, MAT_SPD, isspd));
1563: /* Set up all the "iterative substructuring" common block without computing solvers */
1564: {
1565: Mat temp_mat;
1567: temp_mat = matis->A;
1568: matis->A = pcbddc->local_mat;
1569: PetscCall(PCISSetUp(pc, PETSC_TRUE, PETSC_FALSE));
1570: pcbddc->local_mat = matis->A;
1571: matis->A = temp_mat;
1572: }
1574: /* Analyze interface */
1575: if (!pcbddc->graphanalyzed) {
1576: PetscCall(PCBDDCAnalyzeInterface(pc));
1577: computeconstraintsmatrix = PETSC_TRUE;
1578: PetscCheck(!(pcbddc->adaptive_selection && !pcbddc->use_deluxe_scaling && !pcbddc->mat_graph->twodim), PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot compute the adaptive primal space for a problem with 3D edges without deluxe scaling");
1579: if (pcbddc->compute_nonetflux) {
1580: MatNullSpace nnfnnsp;
1582: PetscCheck(pcbddc->divudotp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Missing divudotp operator");
1583: PetscCall(PCBDDCComputeNoNetFlux(pc->pmat, pcbddc->divudotp, pcbddc->divudotp_trans, pcbddc->divudotp_vl2l, pcbddc->mat_graph, &nnfnnsp));
1584: /* TODO what if a nearnullspace is already attached? */
1585: if (nnfnnsp) {
1586: PetscCall(MatSetNearNullSpace(pc->pmat, nnfnnsp));
1587: PetscCall(MatNullSpaceDestroy(&nnfnnsp));
1588: }
1589: }
1590: }
1591: PetscCall(PetscLogEventEnd(PC_BDDC_Topology[pcbddc->current_level], pc, 0, 0, 0));
1593: /* check existence of a divergence free extension, i.e.
1594: b(v_I,p_0) = 0 for all v_I (raise error if not).
1595: Also, check that PCBDDCBenignGetOrSetP0 works */
1596: if (pcbddc->benign_saddle_point && pcbddc->dbg_flag > 1) PetscCall(PCBDDCBenignCheck(pc, zerodiag));
1597: PetscCall(ISDestroy(&zerodiag));
1599: /* Setup local dirichlet solver ksp_D and sub_schurs solvers */
1600: if (computesubschurs && pcbddc->recompute_topography) PetscCall(PCBDDCInitSubSchurs(pc));
1601: /* SetUp Scaling operator (scaling matrices could be needed in SubSchursSetUp)*/
1602: if (!pcbddc->use_deluxe_scaling) PetscCall(PCBDDCScalingSetUp(pc));
1604: /* finish setup solvers and do adaptive selection of constraints */
1605: sub_schurs = pcbddc->sub_schurs;
1606: if (sub_schurs && sub_schurs->schur_explicit) {
1607: if (computesubschurs) PetscCall(PCBDDCSetUpSubSchurs(pc));
1608: PetscCall(PCBDDCSetUpLocalSolvers(pc, PETSC_TRUE, PETSC_FALSE));
1609: } else {
1610: PetscCall(PCBDDCSetUpLocalSolvers(pc, PETSC_TRUE, PETSC_FALSE));
1611: if (computesubschurs) PetscCall(PCBDDCSetUpSubSchurs(pc));
1612: }
1613: if (pcbddc->adaptive_selection) {
1614: PetscCall(PCBDDCAdaptiveSelection(pc));
1615: computeconstraintsmatrix = PETSC_TRUE;
1616: }
1618: /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */
1619: new_nearnullspace_provided = PETSC_FALSE;
1620: PetscCall(MatGetNearNullSpace(pc->pmat, &nearnullspace));
1621: if (pcbddc->onearnullspace) { /* already used nearnullspace */
1622: if (!nearnullspace) { /* near null space attached to mat has been destroyed */
1623: new_nearnullspace_provided = PETSC_TRUE;
1624: } else {
1625: /* determine if the two nullspaces are different (should be lightweight) */
1626: if (nearnullspace != pcbddc->onearnullspace) {
1627: new_nearnullspace_provided = PETSC_TRUE;
1628: } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */
1629: PetscInt i;
1630: const Vec *nearnullvecs;
1631: PetscObjectState state;
1632: PetscInt nnsp_size;
1633: PetscCall(MatNullSpaceGetVecs(nearnullspace, NULL, &nnsp_size, &nearnullvecs));
1634: for (i = 0; i < nnsp_size; i++) {
1635: PetscCall(PetscObjectStateGet((PetscObject)nearnullvecs[i], &state));
1636: if (pcbddc->onearnullvecs_state[i] != state) {
1637: new_nearnullspace_provided = PETSC_TRUE;
1638: break;
1639: }
1640: }
1641: }
1642: }
1643: } else {
1644: if (!nearnullspace) { /* both nearnullspaces are null */
1645: new_nearnullspace_provided = PETSC_FALSE;
1646: } else { /* nearnullspace attached later */
1647: new_nearnullspace_provided = PETSC_TRUE;
1648: }
1649: }
1651: /* Setup constraints and related work vectors */
1652: /* reset primal space flags */
1653: PetscCall(PetscLogEventBegin(PC_BDDC_LocalWork[pcbddc->current_level], pc, 0, 0, 0));
1654: pcbddc->new_primal_space = PETSC_FALSE;
1655: pcbddc->new_primal_space_local = PETSC_FALSE;
1656: if (computeconstraintsmatrix || new_nearnullspace_provided) {
1657: /* It also sets the primal space flags */
1658: PetscCall(PCBDDCConstraintsSetUp(pc));
1659: }
1660: /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */
1661: PetscCall(PCBDDCSetUpLocalWorkVectors(pc));
1663: if (pcbddc->use_change_of_basis) {
1664: PC_IS *pcis = (PC_IS *)(pc->data);
1666: PetscCall(PCBDDCComputeLocalMatrix(pc, pcbddc->ChangeOfBasisMatrix));
1667: if (pcbddc->benign_change) {
1668: PetscCall(MatDestroy(&pcbddc->benign_B0));
1669: /* pop B0 from pcbddc->local_mat */
1670: PetscCall(PCBDDCBenignPopOrPushB0(pc, PETSC_TRUE));
1671: }
1672: /* get submatrices */
1673: PetscCall(MatDestroy(&pcis->A_IB));
1674: PetscCall(MatDestroy(&pcis->A_BI));
1675: PetscCall(MatDestroy(&pcis->A_BB));
1676: PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcis->is_B_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &pcis->A_BB));
1677: PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcis->is_I_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &pcis->A_IB));
1678: PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcis->is_B_local, pcis->is_I_local, MAT_INITIAL_MATRIX, &pcis->A_BI));
1679: /* set flag in pcis to not reuse submatrices during PCISCreate */
1680: pcis->reusesubmatrices = PETSC_FALSE;
1681: } else if (!pcbddc->user_ChangeOfBasisMatrix && !pcbddc->benign_change) {
1682: PetscCall(MatDestroy(&pcbddc->local_mat));
1683: PetscCall(PetscObjectReference((PetscObject)matis->A));
1684: pcbddc->local_mat = matis->A;
1685: }
1687: /* interface pressure block row for B_C */
1688: PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_lP", (PetscObject *)&lP));
1689: PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_lA", (PetscObject *)&lA));
1690: if (lA && lP) {
1691: PC_IS *pcis = (PC_IS *)pc->data;
1692: Mat B_BI, B_BB, Bt_BI, Bt_BB;
1693: PetscBool issym;
1695: PetscCall(MatIsSymmetric(lA, PETSC_SMALL, &issym));
1696: if (issym) {
1697: PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_I_local, MAT_INITIAL_MATRIX, &B_BI));
1698: PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_B_local, MAT_INITIAL_MATRIX, &B_BB));
1699: PetscCall(MatCreateTranspose(B_BI, &Bt_BI));
1700: PetscCall(MatCreateTranspose(B_BB, &Bt_BB));
1701: } else {
1702: PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_I_local, MAT_INITIAL_MATRIX, &B_BI));
1703: PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_B_local, MAT_INITIAL_MATRIX, &B_BB));
1704: PetscCall(MatCreateSubMatrix(lA, pcis->is_I_local, lP, MAT_INITIAL_MATRIX, &Bt_BI));
1705: PetscCall(MatCreateSubMatrix(lA, pcis->is_B_local, lP, MAT_INITIAL_MATRIX, &Bt_BB));
1706: }
1707: PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_B_BI", (PetscObject)B_BI));
1708: PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_B_BB", (PetscObject)B_BB));
1709: PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_Bt_BI", (PetscObject)Bt_BI));
1710: PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_Bt_BB", (PetscObject)Bt_BB));
1711: PetscCall(MatDestroy(&B_BI));
1712: PetscCall(MatDestroy(&B_BB));
1713: PetscCall(MatDestroy(&Bt_BI));
1714: PetscCall(MatDestroy(&Bt_BB));
1715: }
1716: PetscCall(PetscLogEventEnd(PC_BDDC_LocalWork[pcbddc->current_level], pc, 0, 0, 0));
1718: /* SetUp coarse and local Neumann solvers */
1719: PetscCall(PCBDDCSetUpSolvers(pc));
1720: /* SetUp Scaling operator */
1721: if (pcbddc->use_deluxe_scaling) PetscCall(PCBDDCScalingSetUp(pc));
1723: /* mark topography as done */
1724: pcbddc->recompute_topography = PETSC_FALSE;
1726: /* wrap pcis->A_IB and pcis->A_BI if we did not change explicitly the variables on the pressures */
1727: PetscCall(PCBDDCBenignShellMat(pc, PETSC_FALSE));
1729: if (pcbddc->dbg_flag) {
1730: PetscCall(PetscViewerASCIISubtractTab(pcbddc->dbg_viewer, 2 * pcbddc->current_level));
1731: PetscCall(PetscViewerASCIIPopSynchronized(pcbddc->dbg_viewer));
1732: }
1733: PetscFunctionReturn(PETSC_SUCCESS);
1734: }
1736: // PetscClangLinter pragma disable: -fdoc-sowing-chars
1737: /*
1738: PCApply_BDDC - Applies the BDDC operator to a vector.
1740: Input Parameters:
1741: + pc - the preconditioner context
1742: - r - input vector (global)
1744: Output Parameter:
1745: . z - output vector (global)
1747: Application Interface Routine: PCApply()
1748: */
1749: static PetscErrorCode PCApply_BDDC(PC pc, Vec r, Vec z)
1750: {
1751: PC_IS *pcis = (PC_IS *)(pc->data);
1752: PC_BDDC *pcbddc = (PC_BDDC *)(pc->data);
1753: Mat lA = NULL;
1754: PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B;
1755: const PetscScalar one = 1.0;
1756: const PetscScalar m_one = -1.0;
1757: const PetscScalar zero = 0.0;
1758: /* This code is similar to that provided in nn.c for PCNN
1759: NN interface preconditioner changed to BDDC
1760: Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */
1762: PetscFunctionBegin;
1763: PetscCall(PetscCitationsRegister(citation, &cited));
1764: if (pcbddc->switch_static) PetscCall(MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat, &lA));
1766: if (pcbddc->ChangeOfBasisMatrix) {
1767: Vec swap;
1769: PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, r, pcbddc->work_change));
1770: swap = pcbddc->work_change;
1771: pcbddc->work_change = r;
1772: r = swap;
1773: /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */
1774: if (pcbddc->benign_apply_coarse_only && pcbddc->use_exact_dirichlet_trick && pcbddc->change_interior) {
1775: PetscCall(VecCopy(r, pcis->vec1_global));
1776: PetscCall(VecLockReadPush(pcis->vec1_global));
1777: }
1778: }
1779: if (pcbddc->benign_have_null) { /* get p0 from r */
1780: PetscCall(PCBDDCBenignGetOrSetP0(pc, r, PETSC_TRUE));
1781: }
1782: if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_DIRICHLET && !pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
1783: PetscCall(VecCopy(r, z));
1784: /* First Dirichlet solve */
1785: PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1786: PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1787: /*
1788: Assembling right hand side for BDDC operator
1789: - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE)
1790: - pcis->vec1_B the interface part of the global vector z
1791: */
1792: if (n_D) {
1793: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1794: PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D));
1795: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1796: PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D));
1797: PetscCall(VecScale(pcis->vec2_D, m_one));
1798: if (pcbddc->switch_static) {
1799: PetscCall(VecSet(pcis->vec1_N, 0.));
1800: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1801: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1802: if (!pcbddc->switch_static_change) {
1803: PetscCall(MatMult(lA, pcis->vec1_N, pcis->vec2_N));
1804: } else {
1805: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1806: PetscCall(MatMult(lA, pcis->vec2_N, pcis->vec1_N));
1807: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1808: }
1809: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
1810: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
1811: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1812: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1813: } else {
1814: PetscCall(MatMult(pcis->A_BI, pcis->vec2_D, pcis->vec1_B));
1815: }
1816: } else {
1817: PetscCall(VecSet(pcis->vec1_B, zero));
1818: }
1819: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
1820: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
1821: PetscCall(PCBDDCScalingRestriction(pc, z, pcis->vec1_B));
1822: } else {
1823: if (!pcbddc->benign_apply_coarse_only) PetscCall(PCBDDCScalingRestriction(pc, r, pcis->vec1_B));
1824: }
1825: if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_LUMP) {
1826: PetscCheck(pcbddc->switch_static, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "You forgot to pass -pc_bddc_switch_static");
1827: PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1828: PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1829: }
1831: /* Apply interface preconditioner
1832: input/output vecs: pcis->vec1_B and pcis->vec1_D */
1833: PetscCall(PCBDDCApplyInterfacePreconditioner(pc, PETSC_FALSE));
1835: /* Apply transpose of partition of unity operator */
1836: PetscCall(PCBDDCScalingExtension(pc, pcis->vec1_B, z));
1837: if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_LUMP) {
1838: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec1_D, z, INSERT_VALUES, SCATTER_REVERSE));
1839: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec1_D, z, INSERT_VALUES, SCATTER_REVERSE));
1840: PetscFunctionReturn(PETSC_SUCCESS);
1841: }
1842: /* Second Dirichlet solve and assembling of output */
1843: PetscCall(VecScatterBegin(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1844: PetscCall(VecScatterEnd(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1845: if (n_B) {
1846: if (pcbddc->switch_static) {
1847: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1848: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1849: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1850: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1851: if (!pcbddc->switch_static_change) {
1852: PetscCall(MatMult(lA, pcis->vec1_N, pcis->vec2_N));
1853: } else {
1854: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1855: PetscCall(MatMult(lA, pcis->vec2_N, pcis->vec1_N));
1856: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1857: }
1858: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
1859: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
1860: } else {
1861: PetscCall(MatMult(pcis->A_IB, pcis->vec1_B, pcis->vec3_D));
1862: }
1863: } else if (pcbddc->switch_static) { /* n_B is zero */
1864: if (!pcbddc->switch_static_change) {
1865: PetscCall(MatMult(lA, pcis->vec1_D, pcis->vec3_D));
1866: } else {
1867: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_D, pcis->vec1_N));
1868: PetscCall(MatMult(lA, pcis->vec1_N, pcis->vec2_N));
1869: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec2_N, pcis->vec3_D));
1870: }
1871: }
1872: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1873: PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec3_D, pcis->vec4_D));
1874: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1875: PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec4_D));
1877: if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
1878: if (pcbddc->switch_static) {
1879: PetscCall(VecAXPBYPCZ(pcis->vec2_D, m_one, one, m_one, pcis->vec4_D, pcis->vec1_D));
1880: } else {
1881: PetscCall(VecAXPBY(pcis->vec2_D, m_one, m_one, pcis->vec4_D));
1882: }
1883: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
1884: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
1885: } else {
1886: if (pcbddc->switch_static) {
1887: PetscCall(VecAXPBY(pcis->vec4_D, one, m_one, pcis->vec1_D));
1888: } else {
1889: PetscCall(VecScale(pcis->vec4_D, m_one));
1890: }
1891: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
1892: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
1893: }
1894: if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */
1895: if (pcbddc->benign_apply_coarse_only) PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
1896: PetscCall(PCBDDCBenignGetOrSetP0(pc, z, PETSC_FALSE));
1897: }
1899: if (pcbddc->ChangeOfBasisMatrix) {
1900: pcbddc->work_change = r;
1901: PetscCall(VecCopy(z, pcbddc->work_change));
1902: PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcbddc->work_change, z));
1903: }
1904: PetscFunctionReturn(PETSC_SUCCESS);
1905: }
1907: /*
1908: PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector.
1910: Input Parameters:
1911: + pc - the preconditioner context
1912: - r - input vector (global)
1914: Output Parameter:
1915: . z - output vector (global)
1917: Application Interface Routine: PCApplyTranspose()
1918: */
1919: static PetscErrorCode PCApplyTranspose_BDDC(PC pc, Vec r, Vec z)
1920: {
1921: PC_IS *pcis = (PC_IS *)(pc->data);
1922: PC_BDDC *pcbddc = (PC_BDDC *)(pc->data);
1923: Mat lA = NULL;
1924: PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B;
1925: const PetscScalar one = 1.0;
1926: const PetscScalar m_one = -1.0;
1927: const PetscScalar zero = 0.0;
1929: PetscFunctionBegin;
1930: PetscCall(PetscCitationsRegister(citation, &cited));
1931: if (pcbddc->switch_static) PetscCall(MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat, &lA));
1932: if (pcbddc->ChangeOfBasisMatrix) {
1933: Vec swap;
1935: PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, r, pcbddc->work_change));
1936: swap = pcbddc->work_change;
1937: pcbddc->work_change = r;
1938: r = swap;
1939: /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */
1940: if (pcbddc->benign_apply_coarse_only && pcbddc->exact_dirichlet_trick_app && pcbddc->change_interior) {
1941: PetscCall(VecCopy(r, pcis->vec1_global));
1942: PetscCall(VecLockReadPush(pcis->vec1_global));
1943: }
1944: }
1945: if (pcbddc->benign_have_null) { /* get p0 from r */
1946: PetscCall(PCBDDCBenignGetOrSetP0(pc, r, PETSC_TRUE));
1947: }
1948: if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
1949: PetscCall(VecCopy(r, z));
1950: /* First Dirichlet solve */
1951: PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1952: PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1953: /*
1954: Assembling right hand side for BDDC operator
1955: - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE)
1956: - pcis->vec1_B the interface part of the global vector z
1957: */
1958: if (n_D) {
1959: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1960: PetscCall(KSPSolveTranspose(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D));
1961: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1962: PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D));
1963: PetscCall(VecScale(pcis->vec2_D, m_one));
1964: if (pcbddc->switch_static) {
1965: PetscCall(VecSet(pcis->vec1_N, 0.));
1966: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1967: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1968: if (!pcbddc->switch_static_change) {
1969: PetscCall(MatMultTranspose(lA, pcis->vec1_N, pcis->vec2_N));
1970: } else {
1971: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1972: PetscCall(MatMultTranspose(lA, pcis->vec2_N, pcis->vec1_N));
1973: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1974: }
1975: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
1976: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
1977: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1978: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1979: } else {
1980: PetscCall(MatMultTranspose(pcis->A_IB, pcis->vec2_D, pcis->vec1_B));
1981: }
1982: } else {
1983: PetscCall(VecSet(pcis->vec1_B, zero));
1984: }
1985: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
1986: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
1987: PetscCall(PCBDDCScalingRestriction(pc, z, pcis->vec1_B));
1988: } else {
1989: PetscCall(PCBDDCScalingRestriction(pc, r, pcis->vec1_B));
1990: }
1992: /* Apply interface preconditioner
1993: input/output vecs: pcis->vec1_B and pcis->vec1_D */
1994: PetscCall(PCBDDCApplyInterfacePreconditioner(pc, PETSC_TRUE));
1996: /* Apply transpose of partition of unity operator */
1997: PetscCall(PCBDDCScalingExtension(pc, pcis->vec1_B, z));
1999: /* Second Dirichlet solve and assembling of output */
2000: PetscCall(VecScatterBegin(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2001: PetscCall(VecScatterEnd(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2002: if (n_B) {
2003: if (pcbddc->switch_static) {
2004: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2005: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2006: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2007: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2008: if (!pcbddc->switch_static_change) {
2009: PetscCall(MatMultTranspose(lA, pcis->vec1_N, pcis->vec2_N));
2010: } else {
2011: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2012: PetscCall(MatMultTranspose(lA, pcis->vec2_N, pcis->vec1_N));
2013: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2014: }
2015: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
2016: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
2017: } else {
2018: PetscCall(MatMultTranspose(pcis->A_BI, pcis->vec1_B, pcis->vec3_D));
2019: }
2020: } else if (pcbddc->switch_static) { /* n_B is zero */
2021: if (!pcbddc->switch_static_change) {
2022: PetscCall(MatMultTranspose(lA, pcis->vec1_D, pcis->vec3_D));
2023: } else {
2024: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_D, pcis->vec1_N));
2025: PetscCall(MatMultTranspose(lA, pcis->vec1_N, pcis->vec2_N));
2026: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec2_N, pcis->vec3_D));
2027: }
2028: }
2029: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2030: PetscCall(KSPSolveTranspose(pcbddc->ksp_D, pcis->vec3_D, pcis->vec4_D));
2031: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2032: PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec4_D));
2033: if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
2034: if (pcbddc->switch_static) {
2035: PetscCall(VecAXPBYPCZ(pcis->vec2_D, m_one, one, m_one, pcis->vec4_D, pcis->vec1_D));
2036: } else {
2037: PetscCall(VecAXPBY(pcis->vec2_D, m_one, m_one, pcis->vec4_D));
2038: }
2039: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
2040: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
2041: } else {
2042: if (pcbddc->switch_static) {
2043: PetscCall(VecAXPBY(pcis->vec4_D, one, m_one, pcis->vec1_D));
2044: } else {
2045: PetscCall(VecScale(pcis->vec4_D, m_one));
2046: }
2047: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
2048: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
2049: }
2050: if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */
2051: PetscCall(PCBDDCBenignGetOrSetP0(pc, z, PETSC_FALSE));
2052: }
2053: if (pcbddc->ChangeOfBasisMatrix) {
2054: pcbddc->work_change = r;
2055: PetscCall(VecCopy(z, pcbddc->work_change));
2056: PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcbddc->work_change, z));
2057: }
2058: PetscFunctionReturn(PETSC_SUCCESS);
2059: }
2061: static PetscErrorCode PCReset_BDDC(PC pc)
2062: {
2063: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
2064: PC_IS *pcis = (PC_IS *)pc->data;
2065: KSP kspD, kspR, kspC;
2067: PetscFunctionBegin;
2068: /* free BDDC custom data */
2069: PetscCall(PCBDDCResetCustomization(pc));
2070: /* destroy objects related to topography */
2071: PetscCall(PCBDDCResetTopography(pc));
2072: /* destroy objects for scaling operator */
2073: PetscCall(PCBDDCScalingDestroy(pc));
2074: /* free solvers stuff */
2075: PetscCall(PCBDDCResetSolvers(pc));
2076: /* free global vectors needed in presolve */
2077: PetscCall(VecDestroy(&pcbddc->temp_solution));
2078: PetscCall(VecDestroy(&pcbddc->original_rhs));
2079: /* free data created by PCIS */
2080: PetscCall(PCISReset(pc));
2082: /* restore defaults */
2083: kspD = pcbddc->ksp_D;
2084: kspR = pcbddc->ksp_R;
2085: kspC = pcbddc->coarse_ksp;
2086: PetscCall(PetscMemzero(pc->data, sizeof(*pcbddc)));
2087: pcis->n_neigh = -1;
2088: pcis->scaling_factor = 1.0;
2089: pcis->reusesubmatrices = PETSC_TRUE;
2090: pcbddc->use_local_adj = PETSC_TRUE;
2091: pcbddc->use_vertices = PETSC_TRUE;
2092: pcbddc->use_edges = PETSC_TRUE;
2093: pcbddc->symmetric_primal = PETSC_TRUE;
2094: pcbddc->vertex_size = 1;
2095: pcbddc->recompute_topography = PETSC_TRUE;
2096: pcbddc->coarse_size = -1;
2097: pcbddc->use_exact_dirichlet_trick = PETSC_TRUE;
2098: pcbddc->coarsening_ratio = 8;
2099: pcbddc->coarse_eqs_per_proc = 1;
2100: pcbddc->benign_compute_correction = PETSC_TRUE;
2101: pcbddc->nedfield = -1;
2102: pcbddc->nedglobal = PETSC_TRUE;
2103: pcbddc->graphmaxcount = PETSC_MAX_INT;
2104: pcbddc->sub_schurs_layers = -1;
2105: pcbddc->ksp_D = kspD;
2106: pcbddc->ksp_R = kspR;
2107: pcbddc->coarse_ksp = kspC;
2108: PetscFunctionReturn(PETSC_SUCCESS);
2109: }
2111: static PetscErrorCode PCDestroy_BDDC(PC pc)
2112: {
2113: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
2115: PetscFunctionBegin;
2116: PetscCall(PCReset_BDDC(pc));
2117: PetscCall(KSPDestroy(&pcbddc->ksp_D));
2118: PetscCall(KSPDestroy(&pcbddc->ksp_R));
2119: PetscCall(KSPDestroy(&pcbddc->coarse_ksp));
2120: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDiscreteGradient_C", NULL));
2121: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDivergenceMat_C", NULL));
2122: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetChangeOfBasisMat_C", NULL));
2123: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesLocalIS_C", NULL));
2124: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesIS_C", NULL));
2125: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesLocalIS_C", NULL));
2126: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesIS_C", NULL));
2127: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetCoarseningRatio_C", NULL));
2128: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevel_C", NULL));
2129: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetUseExactDirichlet_C", NULL));
2130: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevels_C", NULL));
2131: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundaries_C", NULL));
2132: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundariesLocal_C", NULL));
2133: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundaries_C", NULL));
2134: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundariesLocal_C", NULL));
2135: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundaries_C", NULL));
2136: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundariesLocal_C", NULL));
2137: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundaries_C", NULL));
2138: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundariesLocal_C", NULL));
2139: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplitting_C", NULL));
2140: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplittingLocal_C", NULL));
2141: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", NULL));
2142: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCCreateFETIDPOperators_C", NULL));
2143: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetRHS_C", NULL));
2144: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetSolution_C", NULL));
2145: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", NULL));
2146: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
2147: PetscCall(PetscFree(pc->data));
2148: PetscFunctionReturn(PETSC_SUCCESS);
2149: }
2151: static PetscErrorCode PCSetCoordinates_BDDC(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
2152: {
2153: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
2154: PCBDDCGraph mat_graph = pcbddc->mat_graph;
2156: PetscFunctionBegin;
2157: PetscCall(PetscFree(mat_graph->coords));
2158: PetscCall(PetscMalloc1(nloc * dim, &mat_graph->coords));
2159: PetscCall(PetscArraycpy(mat_graph->coords, coords, nloc * dim));
2160: mat_graph->cnloc = nloc;
2161: mat_graph->cdim = dim;
2162: mat_graph->cloc = PETSC_FALSE;
2163: /* flg setup */
2164: pcbddc->recompute_topography = PETSC_TRUE;
2165: pcbddc->corner_selected = PETSC_FALSE;
2166: PetscFunctionReturn(PETSC_SUCCESS);
2167: }
2169: static PetscErrorCode PCPreSolveChangeRHS_BDDC(PC pc, PetscBool *change)
2170: {
2171: PetscFunctionBegin;
2172: *change = PETSC_TRUE;
2173: PetscFunctionReturn(PETSC_SUCCESS);
2174: }
2176: static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
2177: {
2178: FETIDPMat_ctx mat_ctx;
2179: Vec work;
2180: PC_IS *pcis;
2181: PC_BDDC *pcbddc;
2183: PetscFunctionBegin;
2184: PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2185: pcis = (PC_IS *)mat_ctx->pc->data;
2186: pcbddc = (PC_BDDC *)mat_ctx->pc->data;
2188: PetscCall(VecSet(fetidp_flux_rhs, 0.0));
2189: /* copy rhs since we may change it during PCPreSolve_BDDC */
2190: if (!pcbddc->original_rhs) PetscCall(VecDuplicate(pcis->vec1_global, &pcbddc->original_rhs));
2191: if (mat_ctx->rhs_flip) {
2192: PetscCall(VecPointwiseMult(pcbddc->original_rhs, standard_rhs, mat_ctx->rhs_flip));
2193: } else {
2194: PetscCall(VecCopy(standard_rhs, pcbddc->original_rhs));
2195: }
2196: if (mat_ctx->g2g_p) {
2197: /* interface pressure rhs */
2198: PetscCall(VecScatterBegin(mat_ctx->g2g_p, fetidp_flux_rhs, pcbddc->original_rhs, INSERT_VALUES, SCATTER_REVERSE));
2199: PetscCall(VecScatterEnd(mat_ctx->g2g_p, fetidp_flux_rhs, pcbddc->original_rhs, INSERT_VALUES, SCATTER_REVERSE));
2200: PetscCall(VecScatterBegin(mat_ctx->g2g_p, standard_rhs, fetidp_flux_rhs, INSERT_VALUES, SCATTER_FORWARD));
2201: PetscCall(VecScatterEnd(mat_ctx->g2g_p, standard_rhs, fetidp_flux_rhs, INSERT_VALUES, SCATTER_FORWARD));
2202: if (!mat_ctx->rhs_flip) PetscCall(VecScale(fetidp_flux_rhs, -1.));
2203: }
2204: /*
2205: change of basis for physical rhs if needed
2206: It also changes the rhs in case of dirichlet boundaries
2207: */
2208: PetscCall(PCPreSolve_BDDC(mat_ctx->pc, NULL, pcbddc->original_rhs, NULL));
2209: if (pcbddc->ChangeOfBasisMatrix) {
2210: PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, pcbddc->original_rhs, pcbddc->work_change));
2211: work = pcbddc->work_change;
2212: } else {
2213: work = pcbddc->original_rhs;
2214: }
2215: /* store vectors for computation of fetidp final solution */
2216: PetscCall(VecScatterBegin(pcis->global_to_D, work, mat_ctx->temp_solution_D, INSERT_VALUES, SCATTER_FORWARD));
2217: PetscCall(VecScatterEnd(pcis->global_to_D, work, mat_ctx->temp_solution_D, INSERT_VALUES, SCATTER_FORWARD));
2218: /* scale rhs since it should be unassembled */
2219: /* TODO use counter scaling? (also below) */
2220: PetscCall(VecScatterBegin(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2221: PetscCall(VecScatterEnd(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2222: /* Apply partition of unity */
2223: PetscCall(VecPointwiseMult(mat_ctx->temp_solution_B, pcis->D, mat_ctx->temp_solution_B));
2224: /* PetscCall(PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B)); */
2225: if (!pcbddc->switch_static) {
2226: /* compute partially subassembled Schur complement right-hand side */
2227: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2228: PetscCall(KSPSolve(pcbddc->ksp_D, mat_ctx->temp_solution_D, pcis->vec1_D));
2229: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2230: /* Cannot propagate up error in KSPSolve() because there is no access to the PC */
2231: PetscCall(MatMult(pcis->A_BI, pcis->vec1_D, pcis->vec1_B));
2232: PetscCall(VecAXPY(mat_ctx->temp_solution_B, -1.0, pcis->vec1_B));
2233: PetscCall(VecSet(work, 0.0));
2234: PetscCall(VecScatterBegin(pcis->global_to_B, mat_ctx->temp_solution_B, work, ADD_VALUES, SCATTER_REVERSE));
2235: PetscCall(VecScatterEnd(pcis->global_to_B, mat_ctx->temp_solution_B, work, ADD_VALUES, SCATTER_REVERSE));
2236: /* PetscCall(PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B)); */
2237: PetscCall(VecScatterBegin(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2238: PetscCall(VecScatterEnd(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2239: PetscCall(VecPointwiseMult(mat_ctx->temp_solution_B, pcis->D, mat_ctx->temp_solution_B));
2240: }
2241: /* BDDC rhs */
2242: PetscCall(VecCopy(mat_ctx->temp_solution_B, pcis->vec1_B));
2243: if (pcbddc->switch_static) PetscCall(VecCopy(mat_ctx->temp_solution_D, pcis->vec1_D));
2244: /* apply BDDC */
2245: PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
2246: PetscCall(PCBDDCApplyInterfacePreconditioner(mat_ctx->pc, PETSC_FALSE));
2247: PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
2249: /* Application of B_delta and assembling of rhs for fetidp fluxes */
2250: PetscCall(MatMult(mat_ctx->B_delta, pcis->vec1_B, mat_ctx->lambda_local));
2251: PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, mat_ctx->lambda_local, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2252: PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, mat_ctx->lambda_local, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2253: /* Add contribution to interface pressures */
2254: if (mat_ctx->l2g_p) {
2255: PetscCall(MatMult(mat_ctx->B_BB, pcis->vec1_B, mat_ctx->vP));
2256: if (pcbddc->switch_static) PetscCall(MatMultAdd(mat_ctx->B_BI, pcis->vec1_D, mat_ctx->vP, mat_ctx->vP));
2257: PetscCall(VecScatterBegin(mat_ctx->l2g_p, mat_ctx->vP, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2258: PetscCall(VecScatterEnd(mat_ctx->l2g_p, mat_ctx->vP, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2259: }
2260: PetscFunctionReturn(PETSC_SUCCESS);
2261: }
2263: /*@
2264: PCBDDCMatFETIDPGetRHS - Compute the right-hand side for a FETI-DP linear system using the physical right-hand side
2266: Collective
2268: Input Parameters:
2269: + fetidp_mat - the FETI-DP matrix object obtained by a call to `PCBDDCCreateFETIDPOperators()`
2270: - standard_rhs - the right-hand side of the original linear system
2272: Output Parameter:
2273: . fetidp_flux_rhs - the right-hand side for the FETI-DP linear system
2275: Level: developer
2277: .seealso: `PCBDDC`, `PCBDDCCreateFETIDPOperators()`, `PCBDDCMatFETIDPGetSolution()`
2278: @*/
2279: PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
2280: {
2281: FETIDPMat_ctx mat_ctx;
2283: PetscFunctionBegin;
2287: PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2288: PetscUseMethod(mat_ctx->pc, "PCBDDCMatFETIDPGetRHS_C", (Mat, Vec, Vec), (fetidp_mat, standard_rhs, fetidp_flux_rhs));
2289: PetscFunctionReturn(PETSC_SUCCESS);
2290: }
2292: static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
2293: {
2294: FETIDPMat_ctx mat_ctx;
2295: PC_IS *pcis;
2296: PC_BDDC *pcbddc;
2297: Vec work;
2299: PetscFunctionBegin;
2300: PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2301: pcis = (PC_IS *)mat_ctx->pc->data;
2302: pcbddc = (PC_BDDC *)mat_ctx->pc->data;
2304: /* apply B_delta^T */
2305: PetscCall(VecSet(pcis->vec1_B, 0.));
2306: PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, fetidp_flux_sol, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
2307: PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, fetidp_flux_sol, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
2308: PetscCall(MatMultTranspose(mat_ctx->B_delta, mat_ctx->lambda_local, pcis->vec1_B));
2309: if (mat_ctx->l2g_p) {
2310: PetscCall(VecScatterBegin(mat_ctx->l2g_p, fetidp_flux_sol, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
2311: PetscCall(VecScatterEnd(mat_ctx->l2g_p, fetidp_flux_sol, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
2312: PetscCall(MatMultAdd(mat_ctx->Bt_BB, mat_ctx->vP, pcis->vec1_B, pcis->vec1_B));
2313: }
2315: /* compute rhs for BDDC application */
2316: PetscCall(VecAYPX(pcis->vec1_B, -1.0, mat_ctx->temp_solution_B));
2317: if (pcbddc->switch_static) {
2318: PetscCall(VecCopy(mat_ctx->temp_solution_D, pcis->vec1_D));
2319: if (mat_ctx->l2g_p) {
2320: PetscCall(VecScale(mat_ctx->vP, -1.));
2321: PetscCall(MatMultAdd(mat_ctx->Bt_BI, mat_ctx->vP, pcis->vec1_D, pcis->vec1_D));
2322: }
2323: }
2325: /* apply BDDC */
2326: PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
2327: PetscCall(PCBDDCApplyInterfacePreconditioner(mat_ctx->pc, PETSC_FALSE));
2329: /* put values into global vector */
2330: if (pcbddc->ChangeOfBasisMatrix) work = pcbddc->work_change;
2331: else work = standard_sol;
2332: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, work, INSERT_VALUES, SCATTER_REVERSE));
2333: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, work, INSERT_VALUES, SCATTER_REVERSE));
2334: if (!pcbddc->switch_static) {
2335: /* compute values into the interior if solved for the partially subassembled Schur complement */
2336: PetscCall(MatMult(pcis->A_IB, pcis->vec1_B, pcis->vec1_D));
2337: PetscCall(VecAYPX(pcis->vec1_D, -1.0, mat_ctx->temp_solution_D));
2338: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2339: PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec1_D));
2340: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2341: /* Cannot propagate up error in KSPSolve() because there is no access to the PC */
2342: }
2344: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec1_D, work, INSERT_VALUES, SCATTER_REVERSE));
2345: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec1_D, work, INSERT_VALUES, SCATTER_REVERSE));
2346: /* add p0 solution to final solution */
2347: PetscCall(PCBDDCBenignGetOrSetP0(mat_ctx->pc, work, PETSC_FALSE));
2348: if (pcbddc->ChangeOfBasisMatrix) PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, work, standard_sol));
2349: PetscCall(PCPostSolve_BDDC(mat_ctx->pc, NULL, NULL, standard_sol));
2350: if (mat_ctx->g2g_p) {
2351: PetscCall(VecScatterBegin(mat_ctx->g2g_p, fetidp_flux_sol, standard_sol, INSERT_VALUES, SCATTER_REVERSE));
2352: PetscCall(VecScatterEnd(mat_ctx->g2g_p, fetidp_flux_sol, standard_sol, INSERT_VALUES, SCATTER_REVERSE));
2353: }
2354: PetscFunctionReturn(PETSC_SUCCESS);
2355: }
2357: static PetscErrorCode PCView_BDDCIPC(PC pc, PetscViewer viewer)
2358: {
2359: BDDCIPC_ctx bddcipc_ctx;
2360: PetscBool isascii;
2362: PetscFunctionBegin;
2363: PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2364: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2365: if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "BDDC interface preconditioner\n"));
2366: PetscCall(PetscViewerASCIIPushTab(viewer));
2367: PetscCall(PCView(bddcipc_ctx->bddc, viewer));
2368: PetscCall(PetscViewerASCIIPopTab(viewer));
2369: PetscFunctionReturn(PETSC_SUCCESS);
2370: }
2372: static PetscErrorCode PCSetUp_BDDCIPC(PC pc)
2373: {
2374: BDDCIPC_ctx bddcipc_ctx;
2375: PetscBool isbddc;
2376: Vec vv;
2377: IS is;
2378: PC_IS *pcis;
2380: PetscFunctionBegin;
2381: PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2382: PetscCall(PetscObjectTypeCompare((PetscObject)bddcipc_ctx->bddc, PCBDDC, &isbddc));
2383: PetscCheck(isbddc, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Invalid type %s. Must be of type bddc", ((PetscObject)bddcipc_ctx->bddc)->type_name);
2384: PetscCall(PCSetUp(bddcipc_ctx->bddc));
2386: /* create interface scatter */
2387: pcis = (PC_IS *)(bddcipc_ctx->bddc->data);
2388: PetscCall(VecScatterDestroy(&bddcipc_ctx->g2l));
2389: PetscCall(MatCreateVecs(pc->pmat, &vv, NULL));
2390: PetscCall(ISRenumber(pcis->is_B_global, NULL, NULL, &is));
2391: PetscCall(VecScatterCreate(vv, is, pcis->vec1_B, NULL, &bddcipc_ctx->g2l));
2392: PetscCall(ISDestroy(&is));
2393: PetscCall(VecDestroy(&vv));
2394: PetscFunctionReturn(PETSC_SUCCESS);
2395: }
2397: static PetscErrorCode PCApply_BDDCIPC(PC pc, Vec r, Vec x)
2398: {
2399: BDDCIPC_ctx bddcipc_ctx;
2400: PC_IS *pcis;
2401: VecScatter tmps;
2403: PetscFunctionBegin;
2404: PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2405: pcis = (PC_IS *)(bddcipc_ctx->bddc->data);
2406: tmps = pcis->global_to_B;
2407: pcis->global_to_B = bddcipc_ctx->g2l;
2408: PetscCall(PCBDDCScalingRestriction(bddcipc_ctx->bddc, r, pcis->vec1_B));
2409: PetscCall(PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc, PETSC_FALSE));
2410: PetscCall(PCBDDCScalingExtension(bddcipc_ctx->bddc, pcis->vec1_B, x));
2411: pcis->global_to_B = tmps;
2412: PetscFunctionReturn(PETSC_SUCCESS);
2413: }
2415: static PetscErrorCode PCApplyTranspose_BDDCIPC(PC pc, Vec r, Vec x)
2416: {
2417: BDDCIPC_ctx bddcipc_ctx;
2418: PC_IS *pcis;
2419: VecScatter tmps;
2421: PetscFunctionBegin;
2422: PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2423: pcis = (PC_IS *)(bddcipc_ctx->bddc->data);
2424: tmps = pcis->global_to_B;
2425: pcis->global_to_B = bddcipc_ctx->g2l;
2426: PetscCall(PCBDDCScalingRestriction(bddcipc_ctx->bddc, r, pcis->vec1_B));
2427: PetscCall(PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc, PETSC_TRUE));
2428: PetscCall(PCBDDCScalingExtension(bddcipc_ctx->bddc, pcis->vec1_B, x));
2429: pcis->global_to_B = tmps;
2430: PetscFunctionReturn(PETSC_SUCCESS);
2431: }
2433: static PetscErrorCode PCDestroy_BDDCIPC(PC pc)
2434: {
2435: BDDCIPC_ctx bddcipc_ctx;
2437: PetscFunctionBegin;
2438: PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2439: PetscCall(PCDestroy(&bddcipc_ctx->bddc));
2440: PetscCall(VecScatterDestroy(&bddcipc_ctx->g2l));
2441: PetscCall(PetscFree(bddcipc_ctx));
2442: PetscFunctionReturn(PETSC_SUCCESS);
2443: }
2445: /*@
2446: PCBDDCMatFETIDPGetSolution - Compute the physical solution using the solution of the FETI-DP linear system
2448: Collective
2450: Input Parameters:
2451: + fetidp_mat - the FETI-DP matrix obtained by a call to `PCBDDCCreateFETIDPOperators()`
2452: - fetidp_flux_sol - the solution of the FETI-DP linear system`
2454: Output Parameter:
2455: . standard_sol - the solution defined on the physical domain
2457: Level: developer
2459: .seealso: `PCBDDC`, `PCBDDCCreateFETIDPOperators()`, `PCBDDCMatFETIDPGetRHS()`
2460: @*/
2461: PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
2462: {
2463: FETIDPMat_ctx mat_ctx;
2465: PetscFunctionBegin;
2469: PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2470: PetscUseMethod(mat_ctx->pc, "PCBDDCMatFETIDPGetSolution_C", (Mat, Vec, Vec), (fetidp_mat, fetidp_flux_sol, standard_sol));
2471: PetscFunctionReturn(PETSC_SUCCESS);
2472: }
2474: static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, PetscBool fully_redundant, const char *prefix, Mat *fetidp_mat, PC *fetidp_pc)
2475: {
2476: FETIDPMat_ctx fetidpmat_ctx;
2477: Mat newmat;
2478: FETIDPPC_ctx fetidppc_ctx;
2479: PC newpc;
2480: MPI_Comm comm;
2482: PetscFunctionBegin;
2483: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
2484: /* FETI-DP matrix */
2485: PetscCall(PCBDDCCreateFETIDPMatContext(pc, &fetidpmat_ctx));
2486: fetidpmat_ctx->fully_redundant = fully_redundant;
2487: PetscCall(PCBDDCSetupFETIDPMatContext(fetidpmat_ctx));
2488: PetscCall(MatCreateShell(comm, fetidpmat_ctx->n, fetidpmat_ctx->n, fetidpmat_ctx->N, fetidpmat_ctx->N, fetidpmat_ctx, &newmat));
2489: PetscCall(PetscObjectSetName((PetscObject)newmat, !fetidpmat_ctx->l2g_lambda_only ? "F" : "G"));
2490: PetscCall(MatShellSetOperation(newmat, MATOP_MULT, (void (*)(void))FETIDPMatMult));
2491: PetscCall(MatShellSetOperation(newmat, MATOP_MULT_TRANSPOSE, (void (*)(void))FETIDPMatMultTranspose));
2492: PetscCall(MatShellSetOperation(newmat, MATOP_DESTROY, (void (*)(void))PCBDDCDestroyFETIDPMat));
2493: /* propagate MatOptions */
2494: {
2495: PC_BDDC *pcbddc = (PC_BDDC *)fetidpmat_ctx->pc->data;
2496: PetscBool isset, issym;
2498: PetscCall(MatIsSymmetricKnown(pc->mat, &isset, &issym));
2499: if ((isset && issym) || pcbddc->symmetric_primal) PetscCall(MatSetOption(newmat, MAT_SYMMETRIC, PETSC_TRUE));
2500: }
2501: PetscCall(MatSetOptionsPrefix(newmat, prefix));
2502: PetscCall(MatAppendOptionsPrefix(newmat, "fetidp_"));
2503: PetscCall(MatSetUp(newmat));
2504: /* FETI-DP preconditioner */
2505: PetscCall(PCBDDCCreateFETIDPPCContext(pc, &fetidppc_ctx));
2506: PetscCall(PCBDDCSetupFETIDPPCContext(newmat, fetidppc_ctx));
2507: PetscCall(PCCreate(comm, &newpc));
2508: PetscCall(PCSetOperators(newpc, newmat, newmat));
2509: PetscCall(PCSetOptionsPrefix(newpc, prefix));
2510: PetscCall(PCAppendOptionsPrefix(newpc, "fetidp_"));
2511: PetscCall(PCSetErrorIfFailure(newpc, pc->erroriffailure));
2512: if (!fetidpmat_ctx->l2g_lambda_only) { /* standard FETI-DP */
2513: PetscCall(PCSetType(newpc, PCSHELL));
2514: PetscCall(PCShellSetName(newpc, "FETI-DP multipliers"));
2515: PetscCall(PCShellSetContext(newpc, fetidppc_ctx));
2516: PetscCall(PCShellSetApply(newpc, FETIDPPCApply));
2517: PetscCall(PCShellSetApplyTranspose(newpc, FETIDPPCApplyTranspose));
2518: PetscCall(PCShellSetView(newpc, FETIDPPCView));
2519: PetscCall(PCShellSetDestroy(newpc, PCBDDCDestroyFETIDPPC));
2520: } else { /* saddle-point FETI-DP */
2521: Mat M;
2522: PetscInt psize;
2523: PetscBool fake = PETSC_FALSE, isfieldsplit;
2525: PetscCall(ISViewFromOptions(fetidpmat_ctx->lagrange, NULL, "-lag_view"));
2526: PetscCall(ISViewFromOptions(fetidpmat_ctx->pressure, NULL, "-press_view"));
2527: PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_PPmat", (PetscObject *)&M));
2528: PetscCall(PCSetType(newpc, PCFIELDSPLIT));
2529: PetscCall(PCFieldSplitSetIS(newpc, "lag", fetidpmat_ctx->lagrange));
2530: PetscCall(PCFieldSplitSetIS(newpc, "p", fetidpmat_ctx->pressure));
2531: PetscCall(PCFieldSplitSetType(newpc, PC_COMPOSITE_SCHUR));
2532: PetscCall(PCFieldSplitSetSchurFactType(newpc, PC_FIELDSPLIT_SCHUR_FACT_DIAG));
2533: PetscCall(ISGetSize(fetidpmat_ctx->pressure, &psize));
2534: if (psize != M->rmap->N) {
2535: Mat M2;
2536: PetscInt lpsize;
2538: fake = PETSC_TRUE;
2539: PetscCall(ISGetLocalSize(fetidpmat_ctx->pressure, &lpsize));
2540: PetscCall(MatCreate(comm, &M2));
2541: PetscCall(MatSetType(M2, MATAIJ));
2542: PetscCall(MatSetSizes(M2, lpsize, lpsize, psize, psize));
2543: PetscCall(MatSetUp(M2));
2544: PetscCall(MatAssemblyBegin(M2, MAT_FINAL_ASSEMBLY));
2545: PetscCall(MatAssemblyEnd(M2, MAT_FINAL_ASSEMBLY));
2546: PetscCall(PCFieldSplitSetSchurPre(newpc, PC_FIELDSPLIT_SCHUR_PRE_USER, M2));
2547: PetscCall(MatDestroy(&M2));
2548: } else {
2549: PetscCall(PCFieldSplitSetSchurPre(newpc, PC_FIELDSPLIT_SCHUR_PRE_USER, M));
2550: }
2551: PetscCall(PCFieldSplitSetSchurScale(newpc, 1.0));
2553: /* we need to setfromoptions and setup here to access the blocks */
2554: PetscCall(PCSetFromOptions(newpc));
2555: PetscCall(PCSetUp(newpc));
2557: /* user may have changed the type (e.g. -fetidp_pc_type none) */
2558: PetscCall(PetscObjectTypeCompare((PetscObject)newpc, PCFIELDSPLIT, &isfieldsplit));
2559: if (isfieldsplit) {
2560: KSP *ksps;
2561: PC ppc, lagpc;
2562: PetscInt nn;
2563: PetscBool ismatis, matisok = PETSC_FALSE, check = PETSC_FALSE;
2565: /* set the solver for the (0,0) block */
2566: PetscCall(PCFieldSplitSchurGetSubKSP(newpc, &nn, &ksps));
2567: if (!nn) { /* not of type PC_COMPOSITE_SCHUR */
2568: PetscCall(PCFieldSplitGetSubKSP(newpc, &nn, &ksps));
2569: if (!fake) { /* pass pmat to the pressure solver */
2570: Mat F;
2572: PetscCall(KSPGetOperators(ksps[1], &F, NULL));
2573: PetscCall(KSPSetOperators(ksps[1], F, M));
2574: }
2575: } else {
2576: PetscBool issym, isset;
2577: Mat S;
2579: PetscCall(PCFieldSplitSchurGetS(newpc, &S));
2580: PetscCall(MatIsSymmetricKnown(newmat, &isset, &issym));
2581: if (isset) PetscCall(MatSetOption(S, MAT_SYMMETRIC, issym));
2582: }
2583: PetscCall(KSPGetPC(ksps[0], &lagpc));
2584: PetscCall(PCSetType(lagpc, PCSHELL));
2585: PetscCall(PCShellSetName(lagpc, "FETI-DP multipliers"));
2586: PetscCall(PCShellSetContext(lagpc, fetidppc_ctx));
2587: PetscCall(PCShellSetApply(lagpc, FETIDPPCApply));
2588: PetscCall(PCShellSetApplyTranspose(lagpc, FETIDPPCApplyTranspose));
2589: PetscCall(PCShellSetView(lagpc, FETIDPPCView));
2590: PetscCall(PCShellSetDestroy(lagpc, PCBDDCDestroyFETIDPPC));
2592: /* Olof's idea: interface Schur complement preconditioner for the mass matrix */
2593: PetscCall(KSPGetPC(ksps[1], &ppc));
2594: if (fake) {
2595: BDDCIPC_ctx bddcipc_ctx;
2596: PetscContainer c;
2598: matisok = PETSC_TRUE;
2600: /* create inner BDDC solver */
2601: PetscCall(PetscNew(&bddcipc_ctx));
2602: PetscCall(PCCreate(comm, &bddcipc_ctx->bddc));
2603: PetscCall(PCSetType(bddcipc_ctx->bddc, PCBDDC));
2604: PetscCall(PCSetOperators(bddcipc_ctx->bddc, M, M));
2605: PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_pCSR", (PetscObject *)&c));
2606: PetscCall(PetscObjectTypeCompare((PetscObject)M, MATIS, &ismatis));
2607: if (c && ismatis) {
2608: Mat lM;
2609: PetscInt *csr, n;
2611: PetscCall(MatISGetLocalMat(M, &lM));
2612: PetscCall(MatGetSize(lM, &n, NULL));
2613: PetscCall(PetscContainerGetPointer(c, (void **)&csr));
2614: PetscCall(PCBDDCSetLocalAdjacencyGraph(bddcipc_ctx->bddc, n, csr, csr + (n + 1), PETSC_COPY_VALUES));
2615: PetscCall(MatISRestoreLocalMat(M, &lM));
2616: }
2617: PetscCall(PCSetOptionsPrefix(bddcipc_ctx->bddc, ((PetscObject)ksps[1])->prefix));
2618: PetscCall(PCSetErrorIfFailure(bddcipc_ctx->bddc, pc->erroriffailure));
2619: PetscCall(PCSetFromOptions(bddcipc_ctx->bddc));
2621: /* wrap the interface application */
2622: PetscCall(PCSetType(ppc, PCSHELL));
2623: PetscCall(PCShellSetName(ppc, "FETI-DP pressure"));
2624: PetscCall(PCShellSetContext(ppc, bddcipc_ctx));
2625: PetscCall(PCShellSetSetUp(ppc, PCSetUp_BDDCIPC));
2626: PetscCall(PCShellSetApply(ppc, PCApply_BDDCIPC));
2627: PetscCall(PCShellSetApplyTranspose(ppc, PCApplyTranspose_BDDCIPC));
2628: PetscCall(PCShellSetView(ppc, PCView_BDDCIPC));
2629: PetscCall(PCShellSetDestroy(ppc, PCDestroy_BDDCIPC));
2630: }
2632: /* determine if we need to assemble M to construct a preconditioner */
2633: if (!matisok) {
2634: PetscCall(PetscObjectTypeCompare((PetscObject)M, MATIS, &ismatis));
2635: PetscCall(PetscObjectTypeCompareAny((PetscObject)ppc, &matisok, PCBDDC, PCJACOBI, PCNONE, PCMG, ""));
2636: if (ismatis && !matisok) PetscCall(MatConvert(M, MATAIJ, MAT_INPLACE_MATRIX, &M));
2637: }
2639: /* run the subproblems to check convergence */
2640: PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)newmat)->prefix, "-check_saddlepoint", &check, NULL));
2641: if (check) {
2642: PetscInt i;
2644: for (i = 0; i < nn; i++) {
2645: KSP kspC;
2646: PC npc;
2647: Mat F, pF;
2648: Vec x, y;
2649: PetscBool isschur, prec = PETSC_TRUE;
2651: PetscCall(KSPCreate(PetscObjectComm((PetscObject)ksps[i]), &kspC));
2652: PetscCall(KSPSetNestLevel(kspC, pc->kspnestlevel));
2653: PetscCall(KSPSetOptionsPrefix(kspC, ((PetscObject)ksps[i])->prefix));
2654: PetscCall(KSPAppendOptionsPrefix(kspC, "check_"));
2655: PetscCall(KSPGetOperators(ksps[i], &F, &pF));
2656: PetscCall(PetscObjectTypeCompare((PetscObject)F, MATSCHURCOMPLEMENT, &isschur));
2657: if (isschur) {
2658: KSP kspS, kspS2;
2659: Mat A00, pA00, A10, A01, A11;
2660: char prefix[256];
2662: PetscCall(MatSchurComplementGetKSP(F, &kspS));
2663: PetscCall(MatSchurComplementGetSubMatrices(F, &A00, &pA00, &A01, &A10, &A11));
2664: PetscCall(MatCreateSchurComplement(A00, pA00, A01, A10, A11, &F));
2665: PetscCall(MatSchurComplementGetKSP(F, &kspS2));
2666: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sschur_", ((PetscObject)kspC)->prefix));
2667: PetscCall(KSPSetOptionsPrefix(kspS2, prefix));
2668: PetscCall(KSPGetPC(kspS2, &npc));
2669: PetscCall(PCSetType(npc, PCKSP));
2670: PetscCall(PCKSPSetKSP(npc, kspS));
2671: PetscCall(KSPSetFromOptions(kspS2));
2672: PetscCall(KSPGetPC(kspS2, &npc));
2673: PetscCall(PCSetUseAmat(npc, PETSC_TRUE));
2674: } else {
2675: PetscCall(PetscObjectReference((PetscObject)F));
2676: }
2677: PetscCall(KSPSetFromOptions(kspC));
2678: PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)kspC)->prefix, "-preconditioned", &prec, NULL));
2679: if (prec) {
2680: PetscCall(KSPGetPC(ksps[i], &npc));
2681: PetscCall(KSPSetPC(kspC, npc));
2682: }
2683: PetscCall(KSPSetOperators(kspC, F, pF));
2684: PetscCall(MatCreateVecs(F, &x, &y));
2685: PetscCall(VecSetRandom(x, NULL));
2686: PetscCall(MatMult(F, x, y));
2687: PetscCall(KSPSolve(kspC, y, x));
2688: PetscCall(KSPCheckSolve(kspC, npc, x));
2689: PetscCall(KSPDestroy(&kspC));
2690: PetscCall(MatDestroy(&F));
2691: PetscCall(VecDestroy(&x));
2692: PetscCall(VecDestroy(&y));
2693: }
2694: }
2695: PetscCall(PetscFree(ksps));
2696: }
2697: }
2698: /* return pointers for objects created */
2699: *fetidp_mat = newmat;
2700: *fetidp_pc = newpc;
2701: PetscFunctionReturn(PETSC_SUCCESS);
2702: }
2704: /*@C
2705: PCBDDCCreateFETIDPOperators - Create FETI-DP operators
2707: Collective
2709: Input Parameters:
2710: + pc - the `PCBDDC` preconditioning context (setup should have been called before)
2711: . fully_redundant - true for a fully redundant set of Lagrange multipliers
2712: - prefix - optional options database prefix for the objects to be created (can be `NULL`)
2714: Output Parameters:
2715: + fetidp_mat - shell FETI-DP matrix object
2716: - fetidp_pc - shell Dirichlet preconditioner for FETI-DP matrix
2718: Level: developer
2720: Note:
2721: Currently the only operations provided for FETI-DP matrix are `MatMult()` and `MatMultTranspose()`
2723: .seealso: `PCBDDC`, `PCBDDCMatFETIDPGetRHS()`, `PCBDDCMatFETIDPGetSolution()`
2724: @*/
2725: PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, PetscBool fully_redundant, const char *prefix, Mat *fetidp_mat, PC *fetidp_pc)
2726: {
2727: PetscFunctionBegin;
2729: PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "You must call PCSetup_BDDC() first");
2730: PetscUseMethod(pc, "PCBDDCCreateFETIDPOperators_C", (PC, PetscBool, const char *, Mat *, PC *), (pc, fully_redundant, prefix, fetidp_mat, fetidp_pc));
2731: PetscFunctionReturn(PETSC_SUCCESS);
2732: }
2734: /*MC
2735: PCBDDC - Balancing Domain Decomposition by Constraints preconditioners
2737: Requires `MATIS` matrices (Pmat) with local matrices (inside the `MATIS`) of type `MATSEQAIJ`, `MATSEQBAIJ` or `MATSEQSBAIJ`
2739: It also works with unsymmetric and indefinite problems.
2741: Unlike 'conventional' interface preconditioners, `PCBDDC` iterates over all degrees of freedom, not just those on the interface. This allows the use
2742: of approximate solvers on the subdomains.
2744: Approximate local solvers are automatically adapted (see [1]) if the user has attached a nullspace object to the subdomain matrices, and informed
2745: `PCBDDC` of using approximate solvers (via the command line).
2747: Boundary nodes are split in vertices, edges and faces classes using information from the local to global mapping of dofs and the local connectivity graph of nodes.
2748: The latter can be customized by using `PCBDDCSetLocalAdjacencyGraph()`
2750: Additional information on dofs can be provided by using `PCBDDCSetDofsSplitting()`, `PCBDDCSetDirichletBoundaries()`, `PCBDDCSetNeumannBoundaries()`, and
2751: `PCBDDCSetPrimalVerticesIS()` and their local counterparts.
2753: Constraints can be customized by attaching a `MatNullSpace` object to the `MATIS` matrix via `MatSetNearNullSpace()`. Non-singular modes are retained via SVD.
2755: Change of basis is performed similarly to [2] when requested. When more than one constraint is present on a single connected component
2756: (i.e. an edge or a face), a robust method based on local QR factorizations is used.
2757: User defined change of basis can be passed to `PCBDDC` with `PCBDDCSetChangeOfBasisMat()`
2759: The PETSc implementation also supports multilevel `PCBDDC` [3]. Coarse grids are partitioned using a `MatPartitioning` object.
2761: Adaptive selection of primal constraints [4] is supported for SPD systems with high-contrast in the coefficients if MUMPS or MKL_PARDISO are present.
2762: Future versions of the code will also consider using PASTIX.
2764: An experimental interface to the FETI-DP method is available. FETI-DP operators could be created using `PCBDDCCreateFETIDPOperators()`.
2765: A stand-alone class for the FETI-DP method will be provided in the next releases.
2767: Options Database Keys:
2768: + -pc_bddc_use_vertices <true> - use or not vertices in primal space
2769: . -pc_bddc_use_edges <true> - use or not edges in primal space
2770: . -pc_bddc_use_faces <false> - use or not faces in primal space
2771: . -pc_bddc_symmetric <true> - symmetric computation of primal basis functions. Specify false for unsymmetric problems
2772: . -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only)
2773: . -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested
2774: . -pc_bddc_switch_static <false> - switches from M_2 (default) to M_3 operator (see reference article [1])
2775: . -pc_bddc_levels <0> - maximum number of levels for multilevel
2776: . -pc_bddc_coarsening_ratio <8> - number of subdomains which will be aggregated together at the coarser level (e.g. H/h ratio at the coarser level, significative only in the multilevel case)
2777: . -pc_bddc_coarse_redistribute <0> - size of a subset of processors where the coarse problem will be remapped (the value is ignored if not at the coarsest level)
2778: . -pc_bddc_use_deluxe_scaling <false> - use deluxe scaling
2779: . -pc_bddc_schur_layers <\-1> - select the economic version of deluxe scaling by specifying the number of layers (-1 corresponds to the original deluxe scaling)
2780: . -pc_bddc_adaptive_threshold <0.0> - when a value different than zero is specified, adaptive selection of constraints is performed on edges and faces (requires deluxe scaling and MUMPS or MKL_PARDISO installed)
2781: - -pc_bddc_check_level <0> - set verbosity level of debugging output
2783: Options for Dirichlet, Neumann or coarse solver can be set using the appropriate options prefix
2784: .vb
2785: -pc_bddc_dirichlet_
2786: -pc_bddc_neumann_
2787: -pc_bddc_coarse_
2788: .ve
2789: e.g. -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. `PCBDDC` uses by default `KSPPREONLY` and `PCLU`.
2791: When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified using the options prefix
2792: .vb
2793: -pc_bddc_dirichlet_lN_
2794: -pc_bddc_neumann_lN_
2795: -pc_bddc_coarse_lN_
2796: .ve
2797: Note that level number ranges from the finest (0) to the coarsest (N).
2798: In order to specify options for the `PCBDDC` operators at the coarser levels (and not for the solvers), prepend -pc_bddc_coarse_ or -pc_bddc_coarse_l
2799: to the option, e.g.
2800: .vb
2801: -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3
2802: .ve
2803: will use a threshold of 5 for constraints' selection at the first coarse level and will redistribute the coarse problem of the first coarse level on 3 processors
2805: References:
2806: + * - C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149--168, March 2007
2807: . * - A. Klawonn and O. B. Widlund. "Dual-Primal FETI Methods for Linear Elasticity", Communications on Pure and Applied Mathematics Volume 59, Issue 11, pages 1523--1572, November 2006
2808: . * - J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", Computing Volume 83, Issue 2--3, pages 55--85, November 2008
2809: - * - C. Pechstein and C. R. Dohrmann. "Modern domain decomposition methods BDDC, deluxe scaling, and an algebraic approach", Seminar talk, Linz, December 2013, http://people.ricam.oeaw.ac.at/c.pechstein/pechstein-bddc2013.pdf
2811: Level: intermediate
2813: Contributed by Stefano Zampini
2815: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `MATIS`, `PCLU`, `PCGAMG`, `PC`, `PCBDDCSetLocalAdjacencyGraph()`, `PCBDDCSetDofsSplitting()`,
2816: `PCBDDCSetDirichletBoundaries()`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCSetPrimalVerticesIS()`, `MatNullSpace`, `MatSetNearNullSpace()`,
2817: `PCBDDCSetChangeOfBasisMat()`, `PCBDDCCreateFETIDPOperators()`, `PCNN`
2818: M*/
2820: PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc)
2821: {
2822: PC_BDDC *pcbddc;
2824: PetscFunctionBegin;
2825: PetscCall(PetscNew(&pcbddc));
2826: pc->data = pcbddc;
2828: PetscCall(PCISInitialize(pc));
2830: /* create local graph structure */
2831: PetscCall(PCBDDCGraphCreate(&pcbddc->mat_graph));
2833: /* BDDC nonzero defaults */
2834: pcbddc->use_nnsp = PETSC_TRUE;
2835: pcbddc->use_local_adj = PETSC_TRUE;
2836: pcbddc->use_vertices = PETSC_TRUE;
2837: pcbddc->use_edges = PETSC_TRUE;
2838: pcbddc->symmetric_primal = PETSC_TRUE;
2839: pcbddc->vertex_size = 1;
2840: pcbddc->recompute_topography = PETSC_TRUE;
2841: pcbddc->coarse_size = -1;
2842: pcbddc->use_exact_dirichlet_trick = PETSC_TRUE;
2843: pcbddc->coarsening_ratio = 8;
2844: pcbddc->coarse_eqs_per_proc = 1;
2845: pcbddc->benign_compute_correction = PETSC_TRUE;
2846: pcbddc->nedfield = -1;
2847: pcbddc->nedglobal = PETSC_TRUE;
2848: pcbddc->graphmaxcount = PETSC_MAX_INT;
2849: pcbddc->sub_schurs_layers = -1;
2850: pcbddc->adaptive_threshold[0] = 0.0;
2851: pcbddc->adaptive_threshold[1] = 0.0;
2853: /* function pointers */
2854: pc->ops->apply = PCApply_BDDC;
2855: pc->ops->applytranspose = PCApplyTranspose_BDDC;
2856: pc->ops->setup = PCSetUp_BDDC;
2857: pc->ops->destroy = PCDestroy_BDDC;
2858: pc->ops->setfromoptions = PCSetFromOptions_BDDC;
2859: pc->ops->view = PCView_BDDC;
2860: pc->ops->applyrichardson = NULL;
2861: pc->ops->applysymmetricleft = NULL;
2862: pc->ops->applysymmetricright = NULL;
2863: pc->ops->presolve = PCPreSolve_BDDC;
2864: pc->ops->postsolve = PCPostSolve_BDDC;
2865: pc->ops->reset = PCReset_BDDC;
2867: /* composing function */
2868: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDiscreteGradient_C", PCBDDCSetDiscreteGradient_BDDC));
2869: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDivergenceMat_C", PCBDDCSetDivergenceMat_BDDC));
2870: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetChangeOfBasisMat_C", PCBDDCSetChangeOfBasisMat_BDDC));
2871: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesLocalIS_C", PCBDDCSetPrimalVerticesLocalIS_BDDC));
2872: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesIS_C", PCBDDCSetPrimalVerticesIS_BDDC));
2873: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesLocalIS_C", PCBDDCGetPrimalVerticesLocalIS_BDDC));
2874: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesIS_C", PCBDDCGetPrimalVerticesIS_BDDC));
2875: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetCoarseningRatio_C", PCBDDCSetCoarseningRatio_BDDC));
2876: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevel_C", PCBDDCSetLevel_BDDC));
2877: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetUseExactDirichlet_C", PCBDDCSetUseExactDirichlet_BDDC));
2878: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevels_C", PCBDDCSetLevels_BDDC));
2879: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundaries_C", PCBDDCSetDirichletBoundaries_BDDC));
2880: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundariesLocal_C", PCBDDCSetDirichletBoundariesLocal_BDDC));
2881: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundaries_C", PCBDDCSetNeumannBoundaries_BDDC));
2882: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundariesLocal_C", PCBDDCSetNeumannBoundariesLocal_BDDC));
2883: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundaries_C", PCBDDCGetDirichletBoundaries_BDDC));
2884: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundariesLocal_C", PCBDDCGetDirichletBoundariesLocal_BDDC));
2885: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundaries_C", PCBDDCGetNeumannBoundaries_BDDC));
2886: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundariesLocal_C", PCBDDCGetNeumannBoundariesLocal_BDDC));
2887: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplitting_C", PCBDDCSetDofsSplitting_BDDC));
2888: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplittingLocal_C", PCBDDCSetDofsSplittingLocal_BDDC));
2889: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", PCBDDCSetLocalAdjacencyGraph_BDDC));
2890: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCCreateFETIDPOperators_C", PCBDDCCreateFETIDPOperators_BDDC));
2891: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetRHS_C", PCBDDCMatFETIDPGetRHS_BDDC));
2892: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetSolution_C", PCBDDCMatFETIDPGetSolution_BDDC));
2893: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", PCPreSolveChangeRHS_BDDC));
2894: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_BDDC));
2895: PetscFunctionReturn(PETSC_SUCCESS);
2896: }
2898: /*@C
2899: PCBDDCInitializePackage - This function initializes everything in the `PCBDDC` package. It is called
2900: from `PCInitializePackage()`.
2902: Level: developer
2904: .seealso: `PetscInitialize()`, `PCBDDCFinalizePackage()`
2905: @*/
2906: PetscErrorCode PCBDDCInitializePackage(void)
2907: {
2908: int i;
2910: PetscFunctionBegin;
2911: if (PCBDDCPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
2912: PCBDDCPackageInitialized = PETSC_TRUE;
2913: PetscCall(PetscRegisterFinalize(PCBDDCFinalizePackage));
2915: /* general events */
2916: PetscCall(PetscLogEventRegister("PCBDDCTopo", PC_CLASSID, &PC_BDDC_Topology[0]));
2917: PetscCall(PetscLogEventRegister("PCBDDCLKSP", PC_CLASSID, &PC_BDDC_LocalSolvers[0]));
2918: PetscCall(PetscLogEventRegister("PCBDDCLWor", PC_CLASSID, &PC_BDDC_LocalWork[0]));
2919: PetscCall(PetscLogEventRegister("PCBDDCCorr", PC_CLASSID, &PC_BDDC_CorrectionSetUp[0]));
2920: PetscCall(PetscLogEventRegister("PCBDDCASet", PC_CLASSID, &PC_BDDC_ApproxSetUp[0]));
2921: PetscCall(PetscLogEventRegister("PCBDDCAApp", PC_CLASSID, &PC_BDDC_ApproxApply[0]));
2922: PetscCall(PetscLogEventRegister("PCBDDCCSet", PC_CLASSID, &PC_BDDC_CoarseSetUp[0]));
2923: PetscCall(PetscLogEventRegister("PCBDDCCKSP", PC_CLASSID, &PC_BDDC_CoarseSolver[0]));
2924: PetscCall(PetscLogEventRegister("PCBDDCAdap", PC_CLASSID, &PC_BDDC_AdaptiveSetUp[0]));
2925: PetscCall(PetscLogEventRegister("PCBDDCScal", PC_CLASSID, &PC_BDDC_Scaling[0]));
2926: PetscCall(PetscLogEventRegister("PCBDDCSchr", PC_CLASSID, &PC_BDDC_Schurs[0]));
2927: PetscCall(PetscLogEventRegister("PCBDDCDirS", PC_CLASSID, &PC_BDDC_Solves[0][0]));
2928: PetscCall(PetscLogEventRegister("PCBDDCNeuS", PC_CLASSID, &PC_BDDC_Solves[0][1]));
2929: PetscCall(PetscLogEventRegister("PCBDDCCoaS", PC_CLASSID, &PC_BDDC_Solves[0][2]));
2930: for (i = 1; i < PETSC_PCBDDC_MAXLEVELS; i++) {
2931: char ename[32];
2933: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCTopo l%02d", i));
2934: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Topology[i]));
2935: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCLKSP l%02d", i));
2936: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_LocalSolvers[i]));
2937: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCLWor l%02d", i));
2938: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_LocalWork[i]));
2939: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCorr l%02d", i));
2940: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_CorrectionSetUp[i]));
2941: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCASet l%02d", i));
2942: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_ApproxSetUp[i]));
2943: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCAApp l%02d", i));
2944: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_ApproxApply[i]));
2945: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCSet l%02d", i));
2946: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_CoarseSetUp[i]));
2947: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCKSP l%02d", i));
2948: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_CoarseSolver[i]));
2949: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCAdap l%02d", i));
2950: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_AdaptiveSetUp[i]));
2951: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCScal l%02d", i));
2952: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Scaling[i]));
2953: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCSchr l%02d", i));
2954: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Schurs[i]));
2955: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCDirS l%02d", i));
2956: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Solves[i][0]));
2957: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCNeuS l%02d", i));
2958: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Solves[i][1]));
2959: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCoaS l%02d", i));
2960: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Solves[i][2]));
2961: }
2962: PetscFunctionReturn(PETSC_SUCCESS);
2963: }
2965: /*@C
2966: PCBDDCFinalizePackage - This function frees everything from the `PCBDDC` package. It is
2967: called from `PetscFinalize()` automatically.
2969: Level: developer
2971: .seealso: `PetscFinalize()`, `PCBDDCInitializePackage()`
2972: @*/
2973: PetscErrorCode PCBDDCFinalizePackage(void)
2974: {
2975: PetscFunctionBegin;
2976: PCBDDCPackageInitialized = PETSC_FALSE;
2977: PetscFunctionReturn(PETSC_SUCCESS);
2978: }