Actual source code: bddc.c
1: #include <petsc/private/pcbddcimpl.h>
2: #include <petsc/private/pcbddcprivateimpl.h>
3: #include <petscblaslapack.h>
5: static PetscBool PCBDDCPackageInitialized = PETSC_FALSE;
7: static PetscBool cited = PETSC_FALSE;
8: static const char citation[] = "@article{ZampiniPCBDDC,\n"
9: "author = {Stefano Zampini},\n"
10: "title = {{PCBDDC}: A Class of Robust Dual-Primal Methods in {PETS}c},\n"
11: "journal = {SIAM Journal on Scientific Computing},\n"
12: "volume = {38},\n"
13: "number = {5},\n"
14: "pages = {S282-S306},\n"
15: "year = {2016},\n"
16: "doi = {10.1137/15M1025785},\n"
17: "URL = {http://dx.doi.org/10.1137/15M1025785},\n"
18: "eprint = {http://dx.doi.org/10.1137/15M1025785}\n"
19: "}\n";
21: PetscLogEvent PC_BDDC_Topology[PETSC_PCBDDC_MAXLEVELS];
22: PetscLogEvent PC_BDDC_LocalSolvers[PETSC_PCBDDC_MAXLEVELS];
23: PetscLogEvent PC_BDDC_LocalWork[PETSC_PCBDDC_MAXLEVELS];
24: PetscLogEvent PC_BDDC_CorrectionSetUp[PETSC_PCBDDC_MAXLEVELS];
25: PetscLogEvent PC_BDDC_ApproxSetUp[PETSC_PCBDDC_MAXLEVELS];
26: PetscLogEvent PC_BDDC_ApproxApply[PETSC_PCBDDC_MAXLEVELS];
27: PetscLogEvent PC_BDDC_CoarseSetUp[PETSC_PCBDDC_MAXLEVELS];
28: PetscLogEvent PC_BDDC_CoarseSolver[PETSC_PCBDDC_MAXLEVELS];
29: PetscLogEvent PC_BDDC_AdaptiveSetUp[PETSC_PCBDDC_MAXLEVELS];
30: PetscLogEvent PC_BDDC_Scaling[PETSC_PCBDDC_MAXLEVELS];
31: PetscLogEvent PC_BDDC_Schurs[PETSC_PCBDDC_MAXLEVELS];
32: PetscLogEvent PC_BDDC_Solves[PETSC_PCBDDC_MAXLEVELS][3];
34: const char *const PCBDDCInterfaceExtTypes[] = {"DIRICHLET", "LUMP", "PCBDDCInterfaceExtType", "PC_BDDC_INTERFACE_EXT_", NULL};
36: static PetscErrorCode PCApply_BDDC(PC, Vec, Vec);
38: static PetscErrorCode PCSetFromOptions_BDDC(PC pc, PetscOptionItems PetscOptionsObject)
39: {
40: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
41: PetscInt nt, i, load_version = PETSC_DECIDE;
42: char load[PETSC_MAX_PATH_LEN] = {'\0'};
43: PetscBool flg;
45: PetscFunctionBegin;
46: PetscOptionsHeadBegin(PetscOptionsObject, "BDDC options");
47: /* Load customization from binary file (debugging) */
48: PetscCall(PetscOptionsString("-pc_bddc_load", "Load customization from file (intended for debug)", "none", load, load, sizeof(load), &flg));
49: PetscCall(PetscOptionsInt("-pc_bddc_load_version", "Version of the customization file to load", "none", load_version, &load_version, NULL));
50: if (flg) {
51: size_t len;
53: PetscCall(PetscStrlen(load, &len));
54: PetscCall(PCBDDCLoadCustomization(pc, len ? load : NULL, load_version));
55: }
56: /* Verbose debugging */
57: PetscCall(PetscOptionsInt("-pc_bddc_check_level", "Verbose output for PCBDDC (intended for debug)", "none", pcbddc->dbg_flag, &pcbddc->dbg_flag, NULL));
58: /* Approximate solvers */
59: PetscCall(PetscOptionsEnum("-pc_bddc_interface_ext_type", "Use DIRICHLET or LUMP to extend interface corrections to interior", "PCBDDCInterfaceExtType", PCBDDCInterfaceExtTypes, (PetscEnum)pcbddc->interface_extension, (PetscEnum *)&pcbddc->interface_extension, NULL));
60: if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_DIRICHLET) {
61: 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));
62: 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));
63: } else {
64: /* This flag is needed/implied by lumping */
65: pcbddc->switch_static = PETSC_TRUE;
66: }
67: 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));
68: 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));
69: /* Primal space customization */
70: 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));
71: PetscCall(PetscOptionsInt("-pc_bddc_local_mat_graph_square", "Square adjacency graph of local mat for interface analysis", "none", pcbddc->local_adj_square, &pcbddc->local_adj_square, NULL));
72: PetscCall(PetscOptionsInt("-pc_bddc_graph_maxcount", "Maximum number of shared subdomains for a connected component", "none", pcbddc->graphmaxcount, &pcbddc->graphmaxcount, NULL));
73: PetscCall(PetscOptionsBool("-pc_bddc_corner_selection", "Activates face-based corner selection", "none", pcbddc->corner_selection, &pcbddc->corner_selection, NULL));
74: PetscCall(PetscOptionsBool("-pc_bddc_use_vertices", "Use or not corner dofs in coarse space", "none", pcbddc->use_vertices, &pcbddc->use_vertices, NULL));
75: PetscCall(PetscOptionsBool("-pc_bddc_use_edges", "Use or not edge constraints in coarse space", "none", pcbddc->use_edges, &pcbddc->use_edges, NULL));
76: PetscCall(PetscOptionsBool("-pc_bddc_use_faces", "Use or not face constraints in coarse space", "none", pcbddc->use_faces, &pcbddc->use_faces, NULL));
77: 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));
78: 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));
79: 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));
80: 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));
81: /* Change of basis */
82: 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));
83: 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));
84: if (!pcbddc->use_change_of_basis) pcbddc->use_change_on_faces = PETSC_FALSE;
85: /* Switch between M_2 (default) and M_3 preconditioners (as defined by C. Dohrmann in the ref. article) */
86: PetscCall(PetscOptionsBool("-pc_bddc_switch_static", "Switch on static condensation ops around the interface preconditioner", "none", pcbddc->switch_static, &pcbddc->switch_static, NULL));
87: 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));
88: i = pcbddc->coarsening_ratio;
89: PetscCall(PetscOptionsInt("-pc_bddc_coarsening_ratio", "Set coarsening ratio used in multilevel coarsening", "PCBDDCSetCoarseningRatio", i, &i, NULL));
90: PetscCall(PCBDDCSetCoarseningRatio(pc, i));
91: i = pcbddc->max_levels;
92: PetscCall(PetscOptionsInt("-pc_bddc_levels", "Set maximum number of levels for multilevel", "PCBDDCSetLevels", i, &i, NULL));
93: PetscCall(PCBDDCSetLevels(pc, i));
94: 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));
95: PetscCall(PetscOptionsBool("-pc_bddc_use_coarse_estimates", "Use estimated eigenvalues for coarse problem", "none", pcbddc->use_coarse_estimates, &pcbddc->use_coarse_estimates, NULL));
96: PetscCall(PetscOptionsBool("-pc_bddc_use_deluxe_scaling", "Use deluxe scaling for BDDC", "none", pcbddc->use_deluxe_scaling, &pcbddc->use_deluxe_scaling, NULL));
97: 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));
98: 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));
99: 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));
100: 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));
101: 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));
102: PetscCall(PetscOptionsBool("-pc_bddc_deluxe_singlemat", "Collapse deluxe operators", "none", pcbddc->deluxe_singlemat, &pcbddc->deluxe_singlemat, NULL));
103: 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));
104: nt = 2;
105: PetscCall(PetscOptionsRealArray("-pc_bddc_adaptive_threshold", "Thresholds to be used for adaptive selection of constraints", "none", pcbddc->adaptive_threshold, &nt, NULL));
106: if (nt == 1) pcbddc->adaptive_threshold[1] = pcbddc->adaptive_threshold[0];
107: PetscCall(PetscOptionsInt("-pc_bddc_adaptive_nmin", "Minimum number of constraints per connected components", "none", pcbddc->adaptive_nmin, &pcbddc->adaptive_nmin, NULL));
108: PetscCall(PetscOptionsInt("-pc_bddc_adaptive_nmax", "Maximum number of constraints per connected components", "none", pcbddc->adaptive_nmax, &pcbddc->adaptive_nmax, NULL));
109: PetscCall(PetscOptionsBool("-pc_bddc_symmetric", "Symmetric computation of primal basis functions", "none", pcbddc->symmetric_primal, &pcbddc->symmetric_primal, NULL));
110: 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));
111: 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));
112: PetscCall(PetscOptionsBool("-pc_bddc_benign_change", "Compute the pressure change of basis explicitly", "none", pcbddc->benign_change_explicit, &pcbddc->benign_change_explicit, NULL));
113: PetscCall(PetscOptionsBool("-pc_bddc_benign_compute_correction", "Compute the benign correction during PreSolve", "none", pcbddc->benign_compute_correction, &pcbddc->benign_compute_correction, NULL));
114: PetscCall(PetscOptionsBool("-pc_bddc_nonetflux", "Automatic computation of no-net-flux quadrature weights", "none", pcbddc->compute_nonetflux, &pcbddc->compute_nonetflux, NULL));
115: PetscCall(PetscOptionsBool("-pc_bddc_detect_disconnected", "Detects disconnected subdomains", "none", pcbddc->detect_disconnected, &pcbddc->detect_disconnected, NULL));
116: 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));
117: 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));
118: PetscOptionsHeadEnd();
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: static PetscErrorCode PCView_BDDC(PC pc, PetscViewer viewer)
123: {
124: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
125: PC_IS *pcis = (PC_IS *)pc->data;
126: PetscBool isascii;
127: PetscSubcomm subcomm;
128: PetscViewer subviewer;
130: PetscFunctionBegin;
131: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
132: /* ASCII viewer */
133: if (isascii) {
134: PetscMPIInt color, rank, size;
135: PetscInt64 loc[7], gsum[6], gmax[6], gmin[6], totbenign;
136: PetscScalar interface_size;
137: PetscReal ratio1 = 0., ratio2 = 0.;
138: Vec counter;
140: if (!pc->setupcalled) PetscCall(PetscViewerASCIIPrintf(viewer, " Partial information available: preconditioner has not been setup yet\n"));
141: PetscCall(PetscViewerASCIIPrintf(viewer, " Use verbose output: %" PetscInt_FMT "\n", pcbddc->dbg_flag));
142: PetscCall(PetscViewerASCIIPrintf(viewer, " Use user-defined CSR: %d\n", !!pcbddc->mat_graph->nvtxs_csr));
143: PetscCall(PetscViewerASCIIPrintf(viewer, " Use local mat graph: %d\n", pcbddc->use_local_adj && !pcbddc->mat_graph->nvtxs_csr));
144: if (pcbddc->mat_graph->twodim) {
145: PetscCall(PetscViewerASCIIPrintf(viewer, " Connectivity graph topological dimension: 2\n"));
146: } else {
147: PetscCall(PetscViewerASCIIPrintf(viewer, " Connectivity graph topological dimension: 3\n"));
148: }
149: if (pcbddc->graphmaxcount != PETSC_INT_MAX) PetscCall(PetscViewerASCIIPrintf(viewer, " Graph max count: %" PetscInt_FMT "\n", pcbddc->graphmaxcount));
150: PetscCall(PetscViewerASCIIPrintf(viewer, " Corner selection: %d (selected %d)\n", pcbddc->corner_selection, pcbddc->corner_selected));
151: PetscCall(PetscViewerASCIIPrintf(viewer, " Use vertices: %d (vertex size %" PetscInt_FMT ")\n", pcbddc->use_vertices, pcbddc->vertex_size));
152: PetscCall(PetscViewerASCIIPrintf(viewer, " Use edges: %d\n", pcbddc->use_edges));
153: PetscCall(PetscViewerASCIIPrintf(viewer, " Use faces: %d\n", pcbddc->use_faces));
154: PetscCall(PetscViewerASCIIPrintf(viewer, " Use true near null space: %d\n", pcbddc->use_nnsp_true));
155: PetscCall(PetscViewerASCIIPrintf(viewer, " Use QR for single constraints on cc: %d\n", pcbddc->use_qr_single));
156: PetscCall(PetscViewerASCIIPrintf(viewer, " Use change of basis on local edge nodes: %d\n", pcbddc->use_change_of_basis));
157: PetscCall(PetscViewerASCIIPrintf(viewer, " Use change of basis on local face nodes: %d\n", pcbddc->use_change_on_faces));
158: PetscCall(PetscViewerASCIIPrintf(viewer, " User defined change of basis matrix: %d\n", !!pcbddc->user_ChangeOfBasisMatrix));
159: PetscCall(PetscViewerASCIIPrintf(viewer, " Has change of basis matrix: %d\n", !!pcbddc->ChangeOfBasisMatrix));
160: PetscCall(PetscViewerASCIIPrintf(viewer, " Eliminate dirichlet boundary dofs: %d\n", pcbddc->eliminate_dirdofs));
161: PetscCall(PetscViewerASCIIPrintf(viewer, " Switch on static condensation ops around the interface preconditioner: %d\n", pcbddc->switch_static));
162: PetscCall(PetscViewerASCIIPrintf(viewer, " Use exact dirichlet trick: %d\n", pcbddc->use_exact_dirichlet_trick));
163: PetscCall(PetscViewerASCIIPrintf(viewer, " Interface extension: %s\n", PCBDDCInterfaceExtTypes[pcbddc->interface_extension]));
164: PetscCall(PetscViewerASCIIPrintf(viewer, " Multilevel max levels: %" PetscInt_FMT "\n", pcbddc->max_levels));
165: PetscCall(PetscViewerASCIIPrintf(viewer, " Multilevel coarsening ratio: %" PetscInt_FMT "\n", pcbddc->coarsening_ratio));
166: PetscCall(PetscViewerASCIIPrintf(viewer, " Use estimated eigs for coarse problem: %d\n", pcbddc->use_coarse_estimates));
167: PetscCall(PetscViewerASCIIPrintf(viewer, " Use deluxe scaling: %d\n", pcbddc->use_deluxe_scaling));
168: PetscCall(PetscViewerASCIIPrintf(viewer, " Use deluxe zerorows: %d\n", pcbddc->deluxe_zerorows));
169: PetscCall(PetscViewerASCIIPrintf(viewer, " Use deluxe singlemat: %d\n", pcbddc->deluxe_singlemat));
170: PetscCall(PetscViewerASCIIPrintf(viewer, " Rebuild interface graph for Schur principal minors: %d\n", pcbddc->sub_schurs_rebuild));
171: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of dofs' layers for the computation of principal minors: %" PetscInt_FMT "\n", pcbddc->sub_schurs_layers));
172: PetscCall(PetscViewerASCIIPrintf(viewer, " Use user CSR graph to compute successive layers: %d\n", pcbddc->sub_schurs_use_useradj));
173: if (pcbddc->adaptive_threshold[1] != pcbddc->adaptive_threshold[0]) {
174: 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]));
175: } else {
176: PetscCall(PetscViewerASCIIPrintf(viewer, " Adaptive constraint selection threshold (active %d, userdefined %d): %g\n", pcbddc->adaptive_selection, pcbddc->adaptive_userdefined, (double)pcbddc->adaptive_threshold[0]));
177: }
178: PetscCall(PetscViewerASCIIPrintf(viewer, " Min constraints / connected component: %" PetscInt_FMT "\n", pcbddc->adaptive_nmin));
179: PetscCall(PetscViewerASCIIPrintf(viewer, " Max constraints / connected component: %" PetscInt_FMT "\n", pcbddc->adaptive_nmax));
180: PetscCall(PetscViewerASCIIPrintf(viewer, " Invert exact Schur complement for adaptive selection: %d\n", pcbddc->sub_schurs_exact_schur));
181: PetscCall(PetscViewerASCIIPrintf(viewer, " Symmetric computation of primal basis functions: %d\n", pcbddc->symmetric_primal));
182: PetscCall(PetscViewerASCIIPrintf(viewer, " Num. Procs. to map coarse adjacency list: %" PetscInt_FMT "\n", pcbddc->coarse_adj_red));
183: PetscCall(PetscViewerASCIIPrintf(viewer, " Coarse eqs per proc (significant at the coarsest level): %" PetscInt_FMT "\n", pcbddc->coarse_eqs_per_proc));
184: PetscCall(PetscViewerASCIIPrintf(viewer, " Detect disconnected: %d (filter %d)\n", pcbddc->detect_disconnected, pcbddc->detect_disconnected_filter));
185: PetscCall(PetscViewerASCIIPrintf(viewer, " Benign subspace trick: %d (change explicit %d)\n", pcbddc->benign_saddle_point, pcbddc->benign_change_explicit));
186: PetscCall(PetscViewerASCIIPrintf(viewer, " Benign subspace trick is active: %d\n", pcbddc->benign_have_null));
187: PetscCall(PetscViewerASCIIPrintf(viewer, " Algebraic computation of no-net-flux: %d\n", pcbddc->compute_nonetflux));
188: if (!pc->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
190: /* compute interface size */
191: PetscCall(VecSet(pcis->vec1_B, 1.0));
192: PetscCall(MatCreateVecs(pc->pmat, &counter, NULL));
193: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, counter, INSERT_VALUES, SCATTER_REVERSE));
194: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, counter, INSERT_VALUES, SCATTER_REVERSE));
195: PetscCall(VecSum(counter, &interface_size));
196: PetscCall(VecDestroy(&counter));
198: /* compute some statistics on the domain decomposition */
199: gsum[0] = 1;
200: gsum[1] = gsum[2] = gsum[3] = gsum[4] = gsum[5] = 0;
201: loc[0] = !!pcis->n;
202: loc[1] = pcis->n - pcis->n_B;
203: loc[2] = pcis->n_B;
204: loc[3] = pcbddc->local_primal_size;
205: loc[4] = pcis->n;
206: loc[5] = pcbddc->n_local_subs > 0 ? pcbddc->n_local_subs : (pcis->n ? 1 : 0);
207: loc[6] = pcbddc->benign_n;
208: PetscCallMPI(MPI_Reduce(loc, gsum, 6, MPIU_INT64, MPI_SUM, 0, PetscObjectComm((PetscObject)pc)));
209: if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = loc[5] = -1;
210: PetscCallMPI(MPI_Reduce(loc, gmax, 6, MPIU_INT64, MPI_MAX, 0, PetscObjectComm((PetscObject)pc)));
211: if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = loc[5] = PETSC_INT_MAX;
212: PetscCallMPI(MPI_Reduce(loc, gmin, 6, MPIU_INT64, MPI_MIN, 0, PetscObjectComm((PetscObject)pc)));
213: PetscCallMPI(MPI_Reduce(&loc[6], &totbenign, 1, MPIU_INT64, MPI_SUM, 0, PetscObjectComm((PetscObject)pc)));
214: if (pcbddc->coarse_size) {
215: ratio1 = pc->pmat->rmap->N / (1. * pcbddc->coarse_size);
216: ratio2 = PetscRealPart(interface_size) / pcbddc->coarse_size;
217: }
218: PetscCall(PetscViewerASCIIPrintf(viewer, "********************************** STATISTICS AT LEVEL %" PetscInt_FMT " **********************************\n", pcbddc->current_level));
219: 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));
220: PetscCall(PetscViewerASCIIPrintf(viewer, " Coarsening ratios: all/coarse %" PetscInt_FMT " interface/coarse %" PetscInt_FMT "\n", (PetscInt)ratio1, (PetscInt)ratio2));
221: PetscCall(PetscViewerASCIIPrintf(viewer, " Active processes : %" PetscInt64_FMT "\n", gsum[0]));
222: PetscCall(PetscViewerASCIIPrintf(viewer, " Total subdomains : %" PetscInt64_FMT "\n", gsum[5]));
223: if (pcbddc->benign_have_null) PetscCall(PetscViewerASCIIPrintf(viewer, " Benign subs : %" PetscInt64_FMT "\n", totbenign));
224: PetscCall(PetscViewerASCIIPrintf(viewer, " Dofs type :\tMIN\tMAX\tMEAN\n"));
225: PetscCall(PetscViewerASCIIPrintf(viewer, " Interior dofs :\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\n", gmin[1], gmax[1], gsum[1] / gsum[0]));
226: PetscCall(PetscViewerASCIIPrintf(viewer, " Interface dofs :\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\n", gmin[2], gmax[2], gsum[2] / gsum[0]));
227: PetscCall(PetscViewerASCIIPrintf(viewer, " Primal dofs :\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\n", gmin[3], gmax[3], gsum[3] / gsum[0]));
228: PetscCall(PetscViewerASCIIPrintf(viewer, " Local dofs :\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\n", gmin[4], gmax[4], gsum[4] / gsum[0]));
229: PetscCall(PetscViewerASCIIPrintf(viewer, " Local subs :\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\n", gmin[5], gmax[5]));
230: PetscCall(PetscViewerFlush(viewer));
232: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
234: /* local solvers */
235: PetscCall(PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)pcbddc->ksp_D), &subviewer));
236: if (rank == 0) {
237: PetscCall(PetscViewerASCIIPrintf(subviewer, "--- Interior solver (rank 0)\n"));
238: PetscCall(PetscViewerASCIIPushTab(subviewer));
239: PetscCall(KSPView(pcbddc->ksp_D, subviewer));
240: PetscCall(PetscViewerASCIIPopTab(subviewer));
241: PetscCall(PetscViewerASCIIPrintf(subviewer, "--- Correction solver (rank 0)\n"));
242: PetscCall(PetscViewerASCIIPushTab(subviewer));
243: PetscCall(KSPView(pcbddc->ksp_R, subviewer));
244: PetscCall(PetscViewerASCIIPopTab(subviewer));
245: PetscCall(PetscViewerFlush(subviewer));
246: }
247: PetscCall(PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)pcbddc->ksp_D), &subviewer));
248: /* the coarse problem can be handled by a different communicator */
249: if (pcbddc->coarse_ksp) color = 1;
250: else color = 0;
251: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
252: PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)pc), &subcomm));
253: PetscCall(PetscSubcommSetNumber(subcomm, PetscMin(size, 2)));
254: PetscCall(PetscSubcommSetTypeGeneral(subcomm, color, rank));
255: PetscCall(PetscViewerGetSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer));
256: if (color == 1) {
257: PetscCall(PetscViewerASCIIPrintf(subviewer, "--- Coarse solver\n"));
258: PetscCall(PetscViewerASCIIPushTab(subviewer));
259: PetscCall(KSPView(pcbddc->coarse_ksp, subviewer));
260: PetscCall(PetscViewerASCIIPopTab(subviewer));
261: PetscCall(PetscViewerFlush(subviewer));
262: }
263: PetscCall(PetscViewerRestoreSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer));
264: PetscCall(PetscSubcommDestroy(&subcomm));
265: PetscCall(PetscViewerFlush(viewer));
266: }
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: static PetscErrorCode PCBDDCSetDiscreteGradient_BDDC(PC pc, Mat G, PetscInt order, PetscInt field, PetscBool global, PetscBool conforming)
271: {
272: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
274: PetscFunctionBegin;
275: PetscCall(PetscObjectReference((PetscObject)G));
276: PetscCall(MatDestroy(&pcbddc->discretegradient));
277: pcbddc->discretegradient = G;
278: pcbddc->nedorder = order > 0 ? order : -order;
279: pcbddc->nedfield = field;
280: pcbddc->nedglobal = global;
281: pcbddc->conforming = conforming;
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: /*@
286: PCBDDCSetDiscreteGradient - Sets the discrete gradient to be used by the `PCBDDC` preconditioner
288: Collective
290: Input Parameters:
291: + pc - the preconditioning context
292: . G - the discrete gradient matrix (in `MATAIJ` format)
293: . order - the order of the Nedelec space (1 for the lowest order, 0 for variable order)
294: . field - the field index of the Nedelec degrees of freedom, or `PETSC_DECIDE` to infer the field
295: . global - `PETSC_TRUE` if the rows of `G` use the global numbering of all degrees of freedom, `PETSC_FALSE` for the Nedelec field only
296: - conforming - `PETSC_TRUE` if the mesh is conforming
298: Level: advanced
300: Notes:
301: The discrete gradient matrix `G` is used to analyze the subdomain edges and should not contain explicitly stored zero entries.
303: If `PCBDDCSetPrimalVerticesIS()` or `PCBDDCSetPrimalVerticesLocalIS()` specifies a Nedelec degree of freedom as primal,
304: all degrees of freedom on the same mesh edge are made primal before the analysis. These degrees of freedom retain
305: their original coordinates in the generated change of basis.
307: If `global` is `PETSC_FALSE`, the numbering of the Nedelec field must preserve the relative order of its degrees of freedom
308: in the global numbering of all fields. That is, `gid[i] < gid[j]` if and only if `geid[i] < geid[j]`, where `gid` is the global
309: numbering of all degrees of freedom and `geid` is the global numbering of the Nedelec field.
311: The `field` index is not used if no field splitting has been specified.
312: If `field` is `PETSC_DECIDE`, `global` must be `PETSC_TRUE`; the Nedelec field is inferred from the rows of `G` with more than one nonzero.
314: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDofsSplitting()`, `PCBDDCSetDofsSplittingLocal()`, `MATAIJ`, `PCBDDCSetDivergenceMat()`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`
315: @*/
316: PetscErrorCode PCBDDCSetDiscreteGradient(PC pc, Mat G, PetscInt order, PetscInt field, PetscBool global, PetscBool conforming)
317: {
318: PetscFunctionBegin;
325: PetscCheckSameComm(pc, 1, G, 2);
326: PetscTryMethod(pc, "PCBDDCSetDiscreteGradient_C", (PC, Mat, PetscInt, PetscInt, PetscBool, PetscBool), (pc, G, order, field, global, conforming));
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: static PetscErrorCode PCBDDCSetDivergenceMat_BDDC(PC pc, Mat divudotp, PetscBool trans, IS vl2l)
331: {
332: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
334: PetscFunctionBegin;
335: PetscCall(PetscObjectReference((PetscObject)divudotp));
336: PetscCall(MatDestroy(&pcbddc->divudotp));
337: pcbddc->divudotp = divudotp;
338: pcbddc->divudotp_trans = trans;
339: pcbddc->compute_nonetflux = (PetscBool)(divudotp != NULL);
340: PetscCall(PetscObjectReference((PetscObject)vl2l));
341: PetscCall(ISDestroy(&pcbddc->divudotp_vl2l));
342: pcbddc->divudotp_vl2l = vl2l;
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: /*@
347: PCBDDCSetDivergenceMat - Sets the discrete divergence operator used by `PCBDDC` to compute net-flux constraints
349: Collective
351: Input Parameters:
352: + pc - the preconditioning context
353: . divudotp - the matrix (must be of type `MATIS`)
354: . trans - `PETSC_FALSE` if pressures are in the test space and velocities are in the trial space, `PETSC_TRUE` for the transpose
355: - vl2l - optional index set mapping local velocity indices in `divudotp` to local indices in the preconditioning matrix, or `NULL`
357: Level: advanced
359: Notes:
360: The operator represents $\int_\Omega (\nabla \cdot \mathbf{u}) p\,dx$ and is used to compute quadrature weights
361: representing the net flux across subdomain boundaries. See {cite}`zampinitu2017` for their use in mixed formulations of Darcy flow.
363: Local indices refer to the local matrices inside the `MATIS` objects. If `vl2l` is `NULL`, the local velocity numbering in
364: `divudotp` must match that of the preconditioning matrix.
366: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDiscreteGradient()`
367: @*/
368: PetscErrorCode PCBDDCSetDivergenceMat(PC pc, Mat divudotp, PetscBool trans, IS vl2l)
369: {
370: PetscBool ismatis;
372: PetscFunctionBegin;
375: PetscCheckSameComm(pc, 1, divudotp, 2);
378: PetscCall(PetscObjectTypeCompare((PetscObject)divudotp, MATIS, &ismatis));
379: PetscCheck(ismatis, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Divergence matrix needs to be of type MATIS");
380: PetscTryMethod(pc, "PCBDDCSetDivergenceMat_C", (PC, Mat, PetscBool, IS), (pc, divudotp, trans, vl2l));
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: static PetscErrorCode PCBDDCSetChangeOfBasisMat_BDDC(PC pc, Mat change, PetscBool interior)
385: {
386: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
388: PetscFunctionBegin;
389: PetscCall(PetscObjectReference((PetscObject)change));
390: PetscCall(MatDestroy(&pcbddc->user_ChangeOfBasisMatrix));
391: pcbddc->user_ChangeOfBasisMatrix = change;
392: pcbddc->change_interior = interior;
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: /*@
397: PCBDDCSetChangeOfBasisMat - Sets a user-defined change of basis for `PCBDDC`
399: Collective
401: Input Parameters:
402: + pc - the preconditioning context
403: . change - the change-of-basis matrix, with the same global and local sizes as the operator
404: - interior - `PETSC_TRUE` if the change of basis modifies interior degrees of freedom
406: Level: intermediate
408: .seealso: [](ch_ksp), `PCBDDC`
409: @*/
410: PetscErrorCode PCBDDCSetChangeOfBasisMat(PC pc, Mat change, PetscBool interior)
411: {
412: PetscFunctionBegin;
415: PetscCheckSameComm(pc, 1, change, 2);
416: if (pc->mat) {
417: PetscInt rows_c, cols_c, rows, cols;
418: PetscCall(MatGetSize(pc->mat, &rows, &cols));
419: PetscCall(MatGetSize(change, &rows_c, &cols_c));
420: 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);
421: 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);
422: PetscCall(MatGetLocalSize(pc->mat, &rows, &cols));
423: PetscCall(MatGetLocalSize(change, &rows_c, &cols_c));
424: 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);
425: 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);
426: }
427: PetscTryMethod(pc, "PCBDDCSetChangeOfBasisMat_C", (PC, Mat, PetscBool), (pc, change, interior));
428: PetscFunctionReturn(PETSC_SUCCESS);
429: }
431: static PetscErrorCode PCBDDCSetPrimalVerticesIS_BDDC(PC pc, IS PrimalVertices)
432: {
433: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
434: PetscBool isequal = PETSC_FALSE;
436: PetscFunctionBegin;
437: PetscCall(PetscObjectReference((PetscObject)PrimalVertices));
438: if (pcbddc->user_primal_vertices) PetscCall(ISEqual(PrimalVertices, pcbddc->user_primal_vertices, &isequal));
439: PetscCall(ISDestroy(&pcbddc->user_primal_vertices));
440: PetscCall(ISDestroy(&pcbddc->user_primal_vertices_local));
441: pcbddc->user_primal_vertices = PrimalVertices;
442: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }
446: /*@
447: PCBDDCSetPrimalVerticesIS - Sets additional user-defined primal vertices in global numbering for `PCBDDC`
449: Collective
451: Input Parameters:
452: + pc - the preconditioning context
453: - PrimalVertices - index set of primal vertices in global numbering (can be empty)
455: Level: intermediate
457: Note:
458: Any process can list any global degree of freedom.
460: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
461: @*/
462: PetscErrorCode PCBDDCSetPrimalVerticesIS(PC pc, IS PrimalVertices)
463: {
464: PetscFunctionBegin;
467: PetscCheckSameComm(pc, 1, PrimalVertices, 2);
468: PetscTryMethod(pc, "PCBDDCSetPrimalVerticesIS_C", (PC, IS), (pc, PrimalVertices));
469: PetscFunctionReturn(PETSC_SUCCESS);
470: }
472: static PetscErrorCode PCBDDCGetPrimalVerticesIS_BDDC(PC pc, IS *is)
473: {
474: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
476: PetscFunctionBegin;
477: *is = pcbddc->user_primal_vertices;
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: /*@
482: PCBDDCGetPrimalVerticesIS - Gets the user-defined primal vertices in global numbering
484: Not Collective
486: Input Parameter:
487: . pc - the preconditioning context
489: Output Parameter:
490: . is - index set of primal vertices in global numbering (`NULL` if not set)
492: Level: intermediate
494: Note:
495: The returned `IS` is owned by `pc`; the caller must not destroy it.
497: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
498: @*/
499: PetscErrorCode PCBDDCGetPrimalVerticesIS(PC pc, IS *is)
500: {
501: PetscFunctionBegin;
503: PetscAssertPointer(is, 2);
504: PetscUseMethod(pc, "PCBDDCGetPrimalVerticesIS_C", (PC, IS *), (pc, is));
505: PetscFunctionReturn(PETSC_SUCCESS);
506: }
508: static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices)
509: {
510: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
511: PetscBool isequal = PETSC_FALSE;
513: PetscFunctionBegin;
514: PetscCall(PetscObjectReference((PetscObject)PrimalVertices));
515: if (pcbddc->user_primal_vertices_local) PetscCall(ISEqual(PrimalVertices, pcbddc->user_primal_vertices_local, &isequal));
516: PetscCall(ISDestroy(&pcbddc->user_primal_vertices));
517: PetscCall(ISDestroy(&pcbddc->user_primal_vertices_local));
518: pcbddc->user_primal_vertices_local = PrimalVertices;
519: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: /*@
524: PCBDDCSetPrimalVerticesLocalIS - Sets additional user-defined primal vertices in local numbering for `PCBDDC`
526: Collective
528: Input Parameters:
529: + pc - the preconditioning context
530: - PrimalVertices - index set of primal vertices in the numbering of the local `MATIS` matrix (can be empty)
532: Level: intermediate
534: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
535: @*/
536: PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices)
537: {
538: PetscFunctionBegin;
541: PetscCheckSameComm(pc, 1, PrimalVertices, 2);
542: PetscTryMethod(pc, "PCBDDCSetPrimalVerticesLocalIS_C", (PC, IS), (pc, PrimalVertices));
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: static PetscErrorCode PCBDDCGetPrimalVerticesLocalIS_BDDC(PC pc, IS *is)
547: {
548: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
550: PetscFunctionBegin;
551: *is = pcbddc->user_primal_vertices_local;
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
555: /*@
556: PCBDDCGetPrimalVerticesLocalIS - Gets the user-defined primal vertices in local numbering
558: Not Collective
560: Input Parameter:
561: . pc - the preconditioning context
563: Output Parameter:
564: . is - index set of primal vertices in the numbering of the local `MATIS` matrix, or `NULL` if unavailable
566: Level: intermediate
568: Notes:
569: The index set is supplied by `PCBDDCSetPrimalVerticesLocalIS()` or obtained from `PCBDDCSetPrimalVerticesIS()` during `PCSetUp()`.
571: The returned `IS` is owned by `pc`; the caller must not destroy it.
573: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`
574: @*/
575: PetscErrorCode PCBDDCGetPrimalVerticesLocalIS(PC pc, IS *is)
576: {
577: PetscFunctionBegin;
579: PetscAssertPointer(is, 2);
580: PetscUseMethod(pc, "PCBDDCGetPrimalVerticesLocalIS_C", (PC, IS *), (pc, is));
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
584: static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc, PetscInt k)
585: {
586: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
588: PetscFunctionBegin;
589: pcbddc->coarsening_ratio = k;
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: /*@
594: PCBDDCSetCoarseningRatio - Sets the coarsening ratio used by multilevel `PCBDDC`
596: Logically Collective
598: Input Parameters:
599: + pc - the preconditioning context
600: - k - target number of process subdomains or local elements per aggregate
602: Options Database Key:
603: . -pc_bddc_coarsening_ratio k - set the coarsening ratio used in multilevel coarsening
605: Level: intermediate
607: Note:
608: Approximately `k` subdomains at the finer level are aggregated into a single subdomain at the coarser level.
609: When a local `MATIS` matrix stores multiple elements, `k` is the target number of local elements per aggregate.
611: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetLevels()`
612: @*/
613: PetscErrorCode PCBDDCSetCoarseningRatio(PC pc, PetscInt k)
614: {
615: PetscFunctionBegin;
618: PetscTryMethod(pc, "PCBDDCSetCoarseningRatio_C", (PC, PetscInt), (pc, k));
619: PetscFunctionReturn(PETSC_SUCCESS);
620: }
622: /* The following functions (PCBDDCSetUseExactDirichlet PCBDDCSetLevel) are not public */
623: static PetscErrorCode PCBDDCSetUseExactDirichlet_BDDC(PC pc, PetscBool flg)
624: {
625: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
627: PetscFunctionBegin;
628: pcbddc->use_exact_dirichlet_trick = flg;
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc, PetscBool flg)
633: {
634: PetscFunctionBegin;
637: PetscTryMethod(pc, "PCBDDCSetUseExactDirichlet_C", (PC, PetscBool), (pc, flg));
638: PetscFunctionReturn(PETSC_SUCCESS);
639: }
641: static PetscErrorCode PCBDDCSetLevel_BDDC(PC pc, PetscInt level)
642: {
643: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
645: PetscFunctionBegin;
646: pcbddc->current_level = level;
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: PetscErrorCode PCBDDCSetLevel(PC pc, PetscInt level)
651: {
652: PetscFunctionBegin;
655: PetscTryMethod(pc, "PCBDDCSetLevel_C", (PC, PetscInt), (pc, level));
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: static PetscErrorCode PCBDDCSetLevels_BDDC(PC pc, PetscInt levels)
660: {
661: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
663: PetscFunctionBegin;
664: PetscCheck(levels < PETSC_PCBDDC_MAXLEVELS, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Maximum number of additional levels for BDDC is %d", PETSC_PCBDDC_MAXLEVELS - 1);
665: pcbddc->max_levels = levels;
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: /*@
670: PCBDDCSetLevels - Sets the maximum number of additional levels allowed for multilevel `PCBDDC`
672: Logically Collective
674: Input Parameters:
675: + pc - the preconditioning context
676: - levels - the maximum number of additional levels
678: Options Database Key:
679: . -pc_bddc_levels levels - set the maximum number of additional levels for multilevel BDDC
681: Level: intermediate
683: Note:
684: The default value is 0, which gives the classical two-level BDDC algorithm.
686: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetCoarseningRatio()`
687: @*/
688: PetscErrorCode PCBDDCSetLevels(PC pc, PetscInt levels)
689: {
690: PetscFunctionBegin;
693: PetscTryMethod(pc, "PCBDDCSetLevels_C", (PC, PetscInt), (pc, levels));
694: PetscFunctionReturn(PETSC_SUCCESS);
695: }
697: #define PCBDDC_CUSTOMIZATION_VERSION_LEGACY 0
698: #define PCBDDC_CUSTOMIZATION_VERSION_LATEST 1
699: #define PCBDDC_CUSTOMIZATION_HEADER_SIZE_LEGACY 11
700: #define PCBDDC_CUSTOMIZATION_HEADER_SIZE 32
702: static PetscErrorCode PCBDDCLoadOrSaveCustomization_Private(PC pc, PetscBool load, const char *outfile, PetscInt version)
703: {
704: PetscInt header_storage[PCBDDC_CUSTOMIZATION_HEADER_SIZE] = {0};
705: PetscInt *header;
706: PetscInt nheader;
707: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
708: PetscViewer viewer;
709: MPI_Comm comm = PetscObjectComm((PetscObject)pc);
711: PetscFunctionBegin;
712: if (!load && version == PETSC_DECIDE) version = PCBDDC_CUSTOMIZATION_VERSION_LATEST;
713: if (version == PCBDDC_CUSTOMIZATION_VERSION_LEGACY) {
714: header = header_storage;
715: nheader = PCBDDC_CUSTOMIZATION_HEADER_SIZE_LEGACY;
716: } else {
717: header = header_storage + 1;
718: nheader = PCBDDC_CUSTOMIZATION_HEADER_SIZE;
719: }
720: PetscCall(PetscViewerBinaryOpen(comm, outfile ? outfile : "bddc_dump.dat", load ? FILE_MODE_READ : FILE_MODE_WRITE, &viewer));
721: if (load) {
722: IS is;
723: Mat A;
725: PetscCall(PetscViewerBinaryRead(viewer, header_storage, nheader, NULL, PETSC_INT));
726: if (version == PETSC_DECIDE) version = header_storage[0];
727: PetscCheck(header[0] == 0 || header[0] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
728: PetscCheck(header[1] == 0 || header[1] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
729: PetscCheck(header[2] >= 0, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
730: PetscCheck(header[3] == 0 || header[3] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
731: PetscCheck(header[4] == 0 || header[4] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
732: PetscCheck(header[5] >= 0, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
733: PetscCheck(header[7] == 0 || header[7] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
734: PetscCheck(header[8] == 0 || header[8] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
735: PetscCheck(header[9] == 0 || header[9] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
736: PetscCheck(header[10] == 0 || header[10] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
737: if (version >= 1) {
738: PetscCheck(header[11] == 0 || header[11] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
739: PetscCheck(header[12] == 0 || header[12] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
740: }
741: if (header[0]) {
742: PetscCall(ISCreate(comm, &is));
743: PetscCall(ISLoad(is, viewer));
744: PetscCall(PCBDDCSetDirichletBoundaries(pc, is));
745: PetscCall(ISDestroy(&is));
746: }
747: if (header[1]) {
748: PetscCall(ISCreate(comm, &is));
749: PetscCall(ISLoad(is, viewer));
750: PetscCall(PCBDDCSetNeumannBoundaries(pc, is));
751: PetscCall(ISDestroy(&is));
752: }
753: if (header[2]) {
754: IS *isarray;
756: PetscCall(PetscMalloc1(header[2], &isarray));
757: for (PetscInt i = 0; i < header[2]; i++) {
758: PetscCall(ISCreate(comm, &isarray[i]));
759: PetscCall(ISLoad(isarray[i], viewer));
760: }
761: PetscCall(PCBDDCSetDofsSplitting(pc, header[2], isarray));
762: for (PetscInt i = 0; i < header[2]; i++) PetscCall(ISDestroy(&isarray[i]));
763: PetscCall(PetscFree(isarray));
764: }
765: if (header[3]) {
766: PetscCall(ISCreate(comm, &is));
767: PetscCall(ISLoad(is, viewer));
768: PetscCall(PCBDDCSetPrimalVerticesIS(pc, is));
769: PetscCall(ISDestroy(&is));
770: }
771: if (header[4]) {
772: PetscCall(MatCreate(comm, &A));
773: PetscCall(MatSetType(A, MATAIJ));
774: PetscCall(MatLoad(A, viewer));
775: PetscCall(PCBDDCSetDiscreteGradient(pc, A, header[5], header[6], (PetscBool)header[7], (PetscBool)header[8]));
776: PetscCall(MatDestroy(&A));
777: }
778: if (header[9]) {
779: PetscCall(MatCreate(comm, &A));
780: PetscCall(MatSetType(A, MATIS));
781: PetscCall(MatLoad(A, viewer));
782: PetscCall(PCBDDCSetDivergenceMat(pc, A, (PetscBool)header[10], NULL));
783: PetscCall(MatDestroy(&A));
784: }
785: if (header[11]) {
786: PetscCheck(pcbddc->discretegradient, comm, PETSC_ERR_ARG_CORRUPT, "Missing discrete gradient");
787: PetscCall(ISCreate(comm, &is));
788: PetscCall(ISLoad(is, viewer));
789: PetscCall(PetscObjectCompose((PetscObject)pcbddc->discretegradient, "_elements_corners", (PetscObject)is));
790: PetscCall(ISDestroy(&is));
791: }
792: if (header[12]) {
793: MatNullSpace nsp;
795: PetscCheck(pcbddc->discretegradient, comm, PETSC_ERR_ARG_CORRUPT, "Missing discrete gradient");
796: PetscCall(MatNullSpaceLoad(viewer, &nsp));
797: PetscCall(MatSetNullSpace(pcbddc->discretegradient, nsp));
798: PetscCall(MatNullSpaceDestroy(&nsp));
799: }
800: if (header[13]) {
801: PetscReal *coords;
802: PetscInt cdim, nl;
803: PetscCount nc;
805: PetscCheck(pc->pmat, comm, PETSC_ERR_ORDER, "Need to set the matrix first with PCSetOperators()");
806: PetscCall(PetscLayoutGetLocalSize(pc->pmat->rmap, &nl));
807: cdim = header[13];
808: nc = (PetscCount)cdim * nl;
810: PetscCall(PetscMalloc1(nc, &coords));
811: PetscCall(PetscViewerBinaryReadAll(viewer, coords, nc, PETSC_DECIDE, PETSC_DECIDE, PETSC_REAL));
812: PetscCall(PCSetCoordinates(pc, cdim, nl, coords));
813: PetscCall(PetscFree(coords));
814: }
815: } else {
816: if (version != PCBDDC_CUSTOMIZATION_VERSION_LEGACY) header_storage[0] = version;
817: header[0] = (PetscInt)!!pcbddc->DirichletBoundariesLocal;
818: header[1] = (PetscInt)!!pcbddc->NeumannBoundariesLocal;
819: header[2] = pcbddc->n_ISForDofsLocal;
820: header[3] = (PetscInt)!!pcbddc->user_primal_vertices_local;
821: header[4] = (PetscInt)!!pcbddc->discretegradient;
822: header[5] = pcbddc->nedorder;
823: header[6] = pcbddc->nedfield;
824: header[7] = (PetscInt)pcbddc->nedglobal;
825: header[8] = (PetscInt)pcbddc->conforming;
826: header[9] = (PetscInt)!!pcbddc->divudotp;
827: header[10] = (PetscInt)pcbddc->divudotp_trans;
828: if (header[4]) header[3] = 0;
830: if (version >= 1) {
831: if (pcbddc->discretegradient) {
832: IS is;
833: MatNullSpace nsp;
835: PetscCall(PetscObjectQuery((PetscObject)pcbddc->discretegradient, "_elements_corners", (PetscObject *)&is));
836: header[11] = (PetscBool)!!is;
837: PetscCall(MatGetNullSpace(pcbddc->discretegradient, &nsp));
838: header[12] = (PetscBool)!!nsp;
839: }
840: header[13] = pcbddc->mat_graph->cdim;
841: }
843: PetscCall(PetscViewerBinaryWrite(viewer, header_storage, nheader, PETSC_INT));
844: if (header[0]) PetscCall(PCBDDCViewGlobalIS(pc, pcbddc->DirichletBoundariesLocal, viewer));
845: if (header[1]) PetscCall(PCBDDCViewGlobalIS(pc, pcbddc->NeumannBoundariesLocal, viewer));
846: for (PetscInt i = 0; i < header[2]; i++) PetscCall(PCBDDCViewGlobalIS(pc, pcbddc->ISForDofsLocal[i], viewer));
847: if (header[3]) PetscCall(PCBDDCViewGlobalIS(pc, pcbddc->user_primal_vertices_local, viewer));
848: if (header[4]) PetscCall(MatView(pcbddc->discretegradient, viewer));
849: if (header[9]) PetscCall(MatView(pcbddc->divudotp, viewer));
850: if (header[11]) {
851: IS is;
853: PetscCall(PetscObjectQuery((PetscObject)pcbddc->discretegradient, "_elements_corners", (PetscObject *)&is));
854: PetscCall(ISView(is, viewer));
855: }
856: if (header[12]) {
857: MatNullSpace nsp;
859: PetscCall(MatGetNullSpace(pcbddc->discretegradient, &nsp));
860: PetscCall(MatNullSpaceView(nsp, viewer));
861: }
862: if (header[13]) {
863: PetscReal *coords = pcbddc->mat_graph->coords;
864: PetscCount nc = (PetscCount)pcbddc->mat_graph->cdim * pc->pmat->rmap->n;
866: if (pcbddc->mat_graph->cloc) {
867: Mat_IS *matis = (Mat_IS *)pc->pmat->data;
868: PetscMPIInt cdimi;
869: MPI_Datatype dimrealtype;
871: PetscCall(PetscMalloc1(nc, &coords));
872: PetscCall(PetscMPIIntCast(pcbddc->mat_graph->cdim, &cdimi));
873: PetscCallMPI(MPI_Type_contiguous(cdimi, MPIU_REAL, &dimrealtype));
874: PetscCallMPI(MPI_Type_commit(&dimrealtype));
875: PetscCall(PetscSFReduceBegin(matis->sf, dimrealtype, pcbddc->mat_graph->coords, coords, MPI_REPLACE));
876: PetscCall(PetscSFReduceEnd(matis->sf, dimrealtype, pcbddc->mat_graph->coords, coords, MPI_REPLACE));
877: PetscCallMPI(MPI_Type_free(&dimrealtype));
878: }
879: PetscCall(PetscViewerBinaryWriteAll(viewer, coords, nc, PETSC_DECIDE, PETSC_DECIDE, PETSC_REAL));
880: if (pcbddc->mat_graph->cloc) PetscCall(PetscFree(coords));
881: }
882: }
883: PetscCall(PetscViewerDestroy(&viewer));
884: PetscFunctionReturn(PETSC_SUCCESS);
885: }
887: static PetscErrorCode PCBDDCLoadCustomization_BDDC(PC pc, const char filename[], PetscInt version)
888: {
889: PetscFunctionBegin;
890: PetscCall(PCBDDCLoadOrSaveCustomization_Private(pc, PETSC_TRUE, filename, version));
891: PetscFunctionReturn(PETSC_SUCCESS);
892: }
894: /*@
895: PCBDDCLoadCustomization - Loads user-defined customization data for `PCBDDC` from a binary file
897: Collective
899: Input Parameters:
900: + pc - the preconditioning context
901: . filename - path to the binary file, or `NULL` for `bddc_dump.dat`
902: - version - file format version, or `PETSC_DECIDE` to detect the version from the file
904: Level: advanced
906: Note:
907: This routine is normally called before `PCSetUp()`.
909: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSaveCustomization()`, `PCSetUp()`
910: @*/
911: PetscErrorCode PCBDDCLoadCustomization(PC pc, const char filename[], PetscInt version)
912: {
913: PetscFunctionBegin;
915: if (filename) PetscAssertPointer(filename, 2);
917: PetscTryMethod(pc, "PCBDDCLoadCustomization_C", (PC, const char[], PetscInt), (pc, filename, version));
918: PetscFunctionReturn(PETSC_SUCCESS);
919: }
921: static PetscErrorCode PCBDDCSaveCustomization_BDDC(PC pc, const char filename[], PetscInt version)
922: {
923: PetscFunctionBegin;
924: PetscCall(PCBDDCLoadOrSaveCustomization_Private(pc, PETSC_FALSE, filename, version));
925: PetscFunctionReturn(PETSC_SUCCESS);
926: }
928: /*@
929: PCBDDCSaveCustomization - Saves user-defined customization data for `PCBDDC` to a binary file
931: Collective
933: Input Parameters:
934: + pc - the preconditioning context
935: . filename - path to the binary file, or `NULL` for `bddc_dump.dat`
936: - version - file format version, or `PETSC_DECIDE` to use the latest version
938: Level: advanced
940: Note:
941: Call `PCSetUp()` before this routine so that global customization data has been converted to the local representation stored in the file.
943: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCLoadCustomization()`, `PCSetUp()`
944: @*/
945: PetscErrorCode PCBDDCSaveCustomization(PC pc, const char filename[], PetscInt version)
946: {
947: PetscFunctionBegin;
949: if (filename) PetscAssertPointer(filename, 2);
951: PetscTryMethod(pc, "PCBDDCSaveCustomization_C", (PC, const char[], PetscInt), (pc, filename, version));
952: PetscFunctionReturn(PETSC_SUCCESS);
953: }
955: static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc, IS DirichletBoundaries)
956: {
957: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
958: PetscBool isequal = PETSC_FALSE;
960: PetscFunctionBegin;
961: PetscCall(PetscObjectReference((PetscObject)DirichletBoundaries));
962: if (pcbddc->DirichletBoundaries) PetscCall(ISEqual(DirichletBoundaries, pcbddc->DirichletBoundaries, &isequal));
963: /* last user setting takes precedence -> destroy any other customization */
964: PetscCall(ISDestroy(&pcbddc->DirichletBoundariesLocal));
965: PetscCall(ISDestroy(&pcbddc->DirichletBoundaries));
966: pcbddc->DirichletBoundaries = DirichletBoundaries;
967: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
968: PetscFunctionReturn(PETSC_SUCCESS);
969: }
971: /*@
972: PCBDDCSetDirichletBoundaries - Sets the Dirichlet boundary degrees of freedom in global numbering
974: Collective
976: Input Parameters:
977: + pc - the preconditioning context
978: - DirichletBoundaries - index set of Dirichlet boundary degrees of freedom in global numbering
980: Level: intermediate
982: Note:
983: Provide this information when Dirichlet conditions have been imposed with `MatZeroRows()` or `MatZeroRowsColumns()`.
984: Any process can list any global degree of freedom.
986: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCGetDirichletBoundaries()`, `PCBDDCSetDirichletBoundariesLocal()`, `MatZeroRows()`, `MatZeroRowsColumns()`
987: @*/
988: PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc, IS DirichletBoundaries)
989: {
990: PetscFunctionBegin;
993: PetscCheckSameComm(pc, 1, DirichletBoundaries, 2);
994: PetscTryMethod(pc, "PCBDDCSetDirichletBoundaries_C", (PC, IS), (pc, DirichletBoundaries));
995: PetscFunctionReturn(PETSC_SUCCESS);
996: }
998: static PetscErrorCode PCBDDCSetDirichletBoundariesLocal_BDDC(PC pc, IS DirichletBoundaries)
999: {
1000: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1001: PetscBool isequal = PETSC_FALSE;
1003: PetscFunctionBegin;
1004: PetscCall(PetscObjectReference((PetscObject)DirichletBoundaries));
1005: if (pcbddc->DirichletBoundariesLocal) PetscCall(ISEqual(DirichletBoundaries, pcbddc->DirichletBoundariesLocal, &isequal));
1006: /* last user setting takes precedence -> destroy any other customization */
1007: PetscCall(ISDestroy(&pcbddc->DirichletBoundariesLocal));
1008: PetscCall(ISDestroy(&pcbddc->DirichletBoundaries));
1009: pcbddc->DirichletBoundariesLocal = DirichletBoundaries;
1010: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1011: PetscFunctionReturn(PETSC_SUCCESS);
1012: }
1014: /*@
1015: PCBDDCSetDirichletBoundariesLocal - Sets the Dirichlet boundary degrees of freedom in local numbering
1017: Collective
1019: Input Parameters:
1020: + pc - the preconditioning context
1021: - DirichletBoundaries - index set of Dirichlet boundary degrees of freedom in the numbering of the local `MATIS` matrix
1023: Level: intermediate
1025: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCGetDirichletBoundariesLocal()`, `PCBDDCSetDirichletBoundaries()`, `MatZeroRows()`, `MatZeroRowsColumns()`
1026: @*/
1027: PetscErrorCode PCBDDCSetDirichletBoundariesLocal(PC pc, IS DirichletBoundaries)
1028: {
1029: PetscFunctionBegin;
1032: PetscCheckSameComm(pc, 1, DirichletBoundaries, 2);
1033: PetscTryMethod(pc, "PCBDDCSetDirichletBoundariesLocal_C", (PC, IS), (pc, DirichletBoundaries));
1034: PetscFunctionReturn(PETSC_SUCCESS);
1035: }
1037: static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc, IS NeumannBoundaries)
1038: {
1039: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1040: PetscBool isequal = PETSC_FALSE;
1042: PetscFunctionBegin;
1043: PetscCall(PetscObjectReference((PetscObject)NeumannBoundaries));
1044: if (pcbddc->NeumannBoundaries) PetscCall(ISEqual(NeumannBoundaries, pcbddc->NeumannBoundaries, &isequal));
1045: /* last user setting takes precedence -> destroy any other customization */
1046: PetscCall(ISDestroy(&pcbddc->NeumannBoundariesLocal));
1047: PetscCall(ISDestroy(&pcbddc->NeumannBoundaries));
1048: pcbddc->NeumannBoundaries = NeumannBoundaries;
1049: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1050: PetscFunctionReturn(PETSC_SUCCESS);
1051: }
1053: /*@
1054: PCBDDCSetNeumannBoundaries - Sets the Neumann boundary degrees of freedom in global numbering
1056: Collective
1058: Input Parameters:
1059: + pc - the preconditioning context
1060: - NeumannBoundaries - index set of Neumann boundary degrees of freedom in global numbering
1062: Level: intermediate
1064: Note:
1065: Any process can list any global degree of freedom.
1067: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCGetNeumannBoundaries()`, `PCBDDCSetNeumannBoundariesLocal()`
1068: @*/
1069: PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc, IS NeumannBoundaries)
1070: {
1071: PetscFunctionBegin;
1074: PetscCheckSameComm(pc, 1, NeumannBoundaries, 2);
1075: PetscTryMethod(pc, "PCBDDCSetNeumannBoundaries_C", (PC, IS), (pc, NeumannBoundaries));
1076: PetscFunctionReturn(PETSC_SUCCESS);
1077: }
1079: static PetscErrorCode PCBDDCSetNeumannBoundariesLocal_BDDC(PC pc, IS NeumannBoundaries)
1080: {
1081: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1082: PetscBool isequal = PETSC_FALSE;
1084: PetscFunctionBegin;
1085: PetscCall(PetscObjectReference((PetscObject)NeumannBoundaries));
1086: if (pcbddc->NeumannBoundariesLocal) PetscCall(ISEqual(NeumannBoundaries, pcbddc->NeumannBoundariesLocal, &isequal));
1087: /* last user setting takes precedence -> destroy any other customization */
1088: PetscCall(ISDestroy(&pcbddc->NeumannBoundariesLocal));
1089: PetscCall(ISDestroy(&pcbddc->NeumannBoundaries));
1090: pcbddc->NeumannBoundariesLocal = NeumannBoundaries;
1091: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1092: PetscFunctionReturn(PETSC_SUCCESS);
1093: }
1095: /*@
1096: PCBDDCSetNeumannBoundariesLocal - Sets the Neumann boundary degrees of freedom in local numbering
1098: Collective
1100: Input Parameters:
1101: + pc - the preconditioning context
1102: - NeumannBoundaries - index set of Neumann boundary degrees of freedom in the numbering of the local `MATIS` matrix
1104: Level: intermediate
1106: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCGetNeumannBoundariesLocal()`, `PCBDDCSetNeumannBoundaries()`
1107: @*/
1108: PetscErrorCode PCBDDCSetNeumannBoundariesLocal(PC pc, IS NeumannBoundaries)
1109: {
1110: PetscFunctionBegin;
1113: PetscCheckSameComm(pc, 1, NeumannBoundaries, 2);
1114: PetscTryMethod(pc, "PCBDDCSetNeumannBoundariesLocal_C", (PC, IS), (pc, NeumannBoundaries));
1115: PetscFunctionReturn(PETSC_SUCCESS);
1116: }
1118: static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc, IS *DirichletBoundaries)
1119: {
1120: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1122: PetscFunctionBegin;
1123: *DirichletBoundaries = pcbddc->DirichletBoundaries;
1124: PetscFunctionReturn(PETSC_SUCCESS);
1125: }
1127: /*@
1128: PCBDDCGetDirichletBoundaries - Gets the Dirichlet boundary degrees of freedom in global numbering
1130: Not Collective
1132: Input Parameter:
1133: . pc - the preconditioning context
1135: Output Parameter:
1136: . DirichletBoundaries - index set of Dirichlet boundary degrees of freedom in global numbering, or `NULL` if not set
1138: Level: intermediate
1140: Notes:
1141: The returned `IS`, if any, is the one supplied to `PCBDDCSetDirichletBoundaries()`.
1143: The returned `IS` is owned by `pc`; the caller must not destroy it.
1145: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDirichletBoundaries()`, `PCBDDCGetDirichletBoundariesLocal()`
1146: @*/
1147: PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc, IS *DirichletBoundaries)
1148: {
1149: PetscFunctionBegin;
1151: PetscUseMethod(pc, "PCBDDCGetDirichletBoundaries_C", (PC, IS *), (pc, DirichletBoundaries));
1152: PetscFunctionReturn(PETSC_SUCCESS);
1153: }
1155: static PetscErrorCode PCBDDCGetDirichletBoundariesLocal_BDDC(PC pc, IS *DirichletBoundaries)
1156: {
1157: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1159: PetscFunctionBegin;
1160: *DirichletBoundaries = pcbddc->DirichletBoundariesLocal;
1161: PetscFunctionReturn(PETSC_SUCCESS);
1162: }
1164: /*@
1165: PCBDDCGetDirichletBoundariesLocal - Gets the Dirichlet boundary degrees of freedom in local numbering
1167: Not Collective
1169: Input Parameter:
1170: . pc - the preconditioning context
1172: Output Parameter:
1173: . DirichletBoundaries - index set of Dirichlet boundary degrees of freedom in the numbering of the local `MATIS` matrix, or `NULL` if unavailable
1175: Level: intermediate
1177: Notes:
1178: The index set is supplied by `PCBDDCSetDirichletBoundariesLocal()` or obtained from `PCBDDCSetDirichletBoundaries()` during `PCSetUp()`.
1179: Setup may update the local index set to make the boundary information consistent across shared degrees of freedom.
1181: The returned `IS` is owned by `pc`; the caller must not destroy it.
1183: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDirichletBoundariesLocal()`, `PCBDDCGetDirichletBoundaries()`, `PCBDDCSetDirichletBoundaries()`
1184: @*/
1185: PetscErrorCode PCBDDCGetDirichletBoundariesLocal(PC pc, IS *DirichletBoundaries)
1186: {
1187: PetscFunctionBegin;
1189: PetscUseMethod(pc, "PCBDDCGetDirichletBoundariesLocal_C", (PC, IS *), (pc, DirichletBoundaries));
1190: PetscFunctionReturn(PETSC_SUCCESS);
1191: }
1193: static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc, IS *NeumannBoundaries)
1194: {
1195: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1197: PetscFunctionBegin;
1198: *NeumannBoundaries = pcbddc->NeumannBoundaries;
1199: PetscFunctionReturn(PETSC_SUCCESS);
1200: }
1202: /*@
1203: PCBDDCGetNeumannBoundaries - Gets the Neumann boundary degrees of freedom in global numbering
1205: Not Collective
1207: Input Parameter:
1208: . pc - the preconditioning context
1210: Output Parameter:
1211: . NeumannBoundaries - index set of Neumann boundary degrees of freedom in global numbering, or `NULL` if not set
1213: Level: intermediate
1215: Notes:
1216: The returned `IS`, if any, is the one supplied to `PCBDDCSetNeumannBoundaries()`.
1218: The returned `IS` is owned by `pc`; the caller must not destroy it.
1220: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCGetNeumannBoundariesLocal()`
1221: @*/
1222: PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc, IS *NeumannBoundaries)
1223: {
1224: PetscFunctionBegin;
1226: PetscUseMethod(pc, "PCBDDCGetNeumannBoundaries_C", (PC, IS *), (pc, NeumannBoundaries));
1227: PetscFunctionReturn(PETSC_SUCCESS);
1228: }
1230: static PetscErrorCode PCBDDCGetNeumannBoundariesLocal_BDDC(PC pc, IS *NeumannBoundaries)
1231: {
1232: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1234: PetscFunctionBegin;
1235: *NeumannBoundaries = pcbddc->NeumannBoundariesLocal;
1236: PetscFunctionReturn(PETSC_SUCCESS);
1237: }
1239: /*@
1240: PCBDDCGetNeumannBoundariesLocal - Gets the Neumann boundary degrees of freedom in local numbering
1242: Not Collective
1244: Input Parameter:
1245: . pc - the preconditioning context
1247: Output Parameter:
1248: . NeumannBoundaries - index set of Neumann boundary degrees of freedom in the numbering of the local `MATIS` matrix, or `NULL` if unavailable
1250: Level: intermediate
1252: Notes:
1253: The index set is supplied by `PCBDDCSetNeumannBoundariesLocal()` or obtained from `PCBDDCSetNeumannBoundaries()` during `PCSetUp()`.
1254: Setup may update the local index set to make the boundary information consistent across shared degrees of freedom.
1256: The returned `IS` is owned by `pc`; the caller must not destroy it.
1258: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCSetNeumannBoundariesLocal()`, `PCBDDCGetNeumannBoundaries()`
1259: @*/
1260: PetscErrorCode PCBDDCGetNeumannBoundariesLocal(PC pc, IS *NeumannBoundaries)
1261: {
1262: PetscFunctionBegin;
1264: PetscUseMethod(pc, "PCBDDCGetNeumannBoundariesLocal_C", (PC, IS *), (pc, NeumannBoundaries));
1265: PetscFunctionReturn(PETSC_SUCCESS);
1266: }
1268: static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs, const PetscInt xadj[], const PetscInt adjncy[], PetscCopyMode copymode)
1269: {
1270: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1271: PCBDDCGraph mat_graph = pcbddc->mat_graph;
1272: PetscBool same_data = PETSC_FALSE;
1274: PetscFunctionBegin;
1275: if (!nvtxs) {
1276: if (copymode == PETSC_OWN_POINTER) {
1277: PetscCall(PetscFree(xadj));
1278: PetscCall(PetscFree(adjncy));
1279: }
1280: PetscCall(PCBDDCGraphResetCSR(mat_graph));
1281: PetscFunctionReturn(PETSC_SUCCESS);
1282: }
1283: if (mat_graph->nvtxs == nvtxs && mat_graph->freecsr) { /* we own the data */
1284: if (mat_graph->xadj == xadj && mat_graph->adjncy == adjncy) same_data = PETSC_TRUE;
1285: if (!same_data && mat_graph->xadj[nvtxs] == xadj[nvtxs]) {
1286: PetscCall(PetscArraycmp(xadj, mat_graph->xadj, nvtxs + 1, &same_data));
1287: if (same_data) PetscCall(PetscArraycmp(adjncy, mat_graph->adjncy, xadj[nvtxs], &same_data));
1288: }
1289: }
1290: if (!same_data) {
1291: /* free old CSR */
1292: PetscCall(PCBDDCGraphResetCSR(mat_graph));
1293: /* get CSR into graph structure */
1294: if (copymode == PETSC_COPY_VALUES) {
1295: PetscCall(PetscMalloc1(nvtxs + 1, &mat_graph->xadj));
1296: PetscCall(PetscMalloc1(xadj[nvtxs], &mat_graph->adjncy));
1297: PetscCall(PetscArraycpy(mat_graph->xadj, xadj, nvtxs + 1));
1298: PetscCall(PetscArraycpy(mat_graph->adjncy, adjncy, xadj[nvtxs]));
1299: mat_graph->freecsr = PETSC_TRUE;
1300: } else if (copymode == PETSC_OWN_POINTER) {
1301: mat_graph->xadj = (PetscInt *)xadj;
1302: mat_graph->adjncy = (PetscInt *)adjncy;
1303: mat_graph->freecsr = PETSC_TRUE;
1304: } else if (copymode == PETSC_USE_POINTER) {
1305: mat_graph->xadj = (PetscInt *)xadj;
1306: mat_graph->adjncy = (PetscInt *)adjncy;
1307: mat_graph->freecsr = PETSC_FALSE;
1308: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported copy mode %d", copymode);
1309: mat_graph->nvtxs_csr = nvtxs;
1310: pcbddc->recompute_topography = PETSC_TRUE;
1311: }
1312: PetscFunctionReturn(PETSC_SUCCESS);
1313: }
1315: /*@
1316: PCBDDCSetLocalAdjacencyGraph - Sets the adjacency graph of the local degrees of freedom in CSR format
1318: Not Collective
1320: Input Parameters:
1321: + pc - the preconditioning context
1322: . nvtxs - number of local graph vertices, equal to the number of local degrees of freedom
1323: . xadj - CSR row offsets, of length `nvtxs` + 1
1324: . adjncy - CSR column indices, of length `xadj[nvtxs]`
1325: - copymode - `PETSC_COPY_VALUES`, `PETSC_USE_POINTER`, or `PETSC_OWN_POINTER`
1327: Level: intermediate
1329: Notes:
1330: Graph vertices use the numbering of the local `MATIS` matrix. The CSR arrays use zero-based indexing.
1331: A degree of freedom `i` is considered connected to all local degrees of freedom if `xadj[i+1] - xadj[i] == 1`
1332: and `adjncy[xadj[i]]` is negative.
1334: With `PETSC_COPY_VALUES`, the arrays are copied. With `PETSC_OWN_POINTER`, ownership of the arrays is transferred to `pc`;
1335: the arrays must have been allocated with `PetscMalloc()`. With `PETSC_USE_POINTER`, the caller must keep the arrays valid
1336: and unchanged until the graph is replaced or the preconditioner is reset or destroyed.
1338: Passing `nvtxs` equal to 0 clears the user-defined graph.
1340: .seealso: [](ch_ksp), `PCBDDC`, `PetscCopyMode`
1341: @*/
1342: PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc, PetscInt nvtxs, const PetscInt xadj[], const PetscInt adjncy[], PetscCopyMode copymode)
1343: {
1344: PetscBool f = PETSC_FALSE;
1346: PetscFunctionBegin;
1348: if (nvtxs) {
1349: PetscAssertPointer(xadj, 3);
1350: if (xadj[nvtxs]) PetscAssertPointer(adjncy, 4);
1351: }
1352: PetscTryMethod(pc, "PCBDDCSetLocalAdjacencyGraph_C", (PC, PetscInt, const PetscInt[], const PetscInt[], PetscCopyMode), (pc, nvtxs, xadj, adjncy, copymode));
1353: /* free arrays if PCBDDC is not the PC type */
1354: PetscCall(PetscObjectHasFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", &f));
1355: if (!f && copymode == PETSC_OWN_POINTER) {
1356: PetscCall(PetscFree(xadj));
1357: PetscCall(PetscFree(adjncy));
1358: }
1359: PetscFunctionReturn(PETSC_SUCCESS);
1360: }
1362: static PetscErrorCode PCBDDCSetDofsSplittingLocal_BDDC(PC pc, PetscInt n_is, IS ISForDofs[])
1363: {
1364: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1365: PetscInt i;
1366: PetscBool isequal = PETSC_FALSE;
1368: PetscFunctionBegin;
1369: if (pcbddc->n_ISForDofsLocal == n_is) {
1370: for (i = 0; i < n_is; i++) {
1371: PetscBool isequalt;
1372: PetscCall(ISEqual(ISForDofs[i], pcbddc->ISForDofsLocal[i], &isequalt));
1373: if (!isequalt) break;
1374: }
1375: if (i == n_is) isequal = PETSC_TRUE;
1376: }
1377: for (i = 0; i < n_is; i++) PetscCall(PetscObjectReference((PetscObject)ISForDofs[i]));
1378: /* Destroy ISes if they were already set */
1379: for (i = 0; i < pcbddc->n_ISForDofsLocal; i++) PetscCall(ISDestroy(&pcbddc->ISForDofsLocal[i]));
1380: PetscCall(PetscFree(pcbddc->ISForDofsLocal));
1381: /* last user setting takes precedence -> destroy any other customization */
1382: for (i = 0; i < pcbddc->n_ISForDofs; i++) PetscCall(ISDestroy(&pcbddc->ISForDofs[i]));
1383: PetscCall(PetscFree(pcbddc->ISForDofs));
1384: pcbddc->n_ISForDofs = 0;
1385: /* allocate space then set */
1386: if (n_is) PetscCall(PetscMalloc1(n_is, &pcbddc->ISForDofsLocal));
1387: for (i = 0; i < n_is; i++) pcbddc->ISForDofsLocal[i] = ISForDofs[i];
1388: pcbddc->n_ISForDofsLocal = n_is;
1389: if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
1390: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1391: PetscFunctionReturn(PETSC_SUCCESS);
1392: }
1394: /*@
1395: PCBDDCSetDofsSplittingLocal - Sets the fields of the local subdomain matrix
1397: Collective
1399: Input Parameters:
1400: + pc - the preconditioning context
1401: . n_is - number of index sets defining the fields, must be the same on all MPI processes
1402: - ISForDofs - array of `n_is` index sets describing the fields in the numbering of the local `MATIS` matrix
1404: Level: intermediate
1406: Note:
1407: Degrees of freedom not listed in any index set belong to the complement field.
1409: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDofsSplitting()`
1410: @*/
1411: PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc, PetscInt n_is, IS ISForDofs[])
1412: {
1413: PetscInt i;
1415: PetscFunctionBegin;
1418: for (i = 0; i < n_is; i++) {
1419: PetscCheckSameComm(pc, 1, ISForDofs[i], 3);
1421: }
1422: PetscTryMethod(pc, "PCBDDCSetDofsSplittingLocal_C", (PC, PetscInt, IS[]), (pc, n_is, ISForDofs));
1423: PetscFunctionReturn(PETSC_SUCCESS);
1424: }
1426: static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc, PetscInt n_is, IS ISForDofs[])
1427: {
1428: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1429: PetscInt i;
1430: PetscBool isequal = PETSC_FALSE;
1432: PetscFunctionBegin;
1433: if (pcbddc->n_ISForDofs == n_is) {
1434: for (i = 0; i < n_is; i++) {
1435: PetscBool isequalt;
1436: PetscCall(ISEqual(ISForDofs[i], pcbddc->ISForDofs[i], &isequalt));
1437: if (!isequalt) break;
1438: }
1439: if (i == n_is) isequal = PETSC_TRUE;
1440: }
1441: for (i = 0; i < n_is; i++) PetscCall(PetscObjectReference((PetscObject)ISForDofs[i]));
1442: /* Destroy ISes if they were already set */
1443: for (i = 0; i < pcbddc->n_ISForDofs; i++) PetscCall(ISDestroy(&pcbddc->ISForDofs[i]));
1444: PetscCall(PetscFree(pcbddc->ISForDofs));
1445: /* last user setting takes precedence -> destroy any other customization */
1446: for (i = 0; i < pcbddc->n_ISForDofsLocal; i++) PetscCall(ISDestroy(&pcbddc->ISForDofsLocal[i]));
1447: PetscCall(PetscFree(pcbddc->ISForDofsLocal));
1448: pcbddc->n_ISForDofsLocal = 0;
1449: /* allocate space then set */
1450: if (n_is) PetscCall(PetscMalloc1(n_is, &pcbddc->ISForDofs));
1451: for (i = 0; i < n_is; i++) pcbddc->ISForDofs[i] = ISForDofs[i];
1452: pcbddc->n_ISForDofs = n_is;
1453: if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
1454: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1455: PetscFunctionReturn(PETSC_SUCCESS);
1456: }
1458: /*@
1459: PCBDDCSetDofsSplitting - Sets the fields of the global matrix
1461: Collective
1463: Input Parameters:
1464: + pc - the preconditioning context
1465: . n_is - number of index sets defining the fields, must be the same on all MPI processes
1466: - ISForDofs - array of `n_is` index sets describing the fields in global numbering
1468: Level: intermediate
1470: Note:
1471: Any process can list any global degree of freedom. Degrees of freedom not listed in any index set belong to the complement field.
1473: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDofsSplittingLocal()`
1474: @*/
1475: PetscErrorCode PCBDDCSetDofsSplitting(PC pc, PetscInt n_is, IS ISForDofs[])
1476: {
1477: PetscInt i;
1479: PetscFunctionBegin;
1482: for (i = 0; i < n_is; i++) {
1484: PetscCheckSameComm(pc, 1, ISForDofs[i], 3);
1485: }
1486: PetscTryMethod(pc, "PCBDDCSetDofsSplitting_C", (PC, PetscInt, IS[]), (pc, n_is, ISForDofs));
1487: PetscFunctionReturn(PETSC_SUCCESS);
1488: }
1490: static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
1491: {
1492: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1493: PC_IS *pcis = (PC_IS *)pc->data;
1494: Vec used_vec;
1495: PetscBool iscg, save_rhs = PETSC_TRUE, benign_correction_computed;
1497: PetscFunctionBegin;
1498: /* if we are working with CG, one dirichlet solve can be avoided during Krylov iterations */
1499: if (ksp) {
1500: PetscCall(PetscObjectTypeCompareAny((PetscObject)ksp, &iscg, KSPCG, KSPGROPPCG, KSPPIPECG, KSPPIPELCG, KSPPIPECGRR, ""));
1501: if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || !iscg || pc->mat != pc->pmat) PetscCall(PCBDDCSetUseExactDirichlet(pc, PETSC_FALSE));
1502: }
1503: if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || pc->mat != pc->pmat) PetscCall(PCBDDCSetUseExactDirichlet(pc, PETSC_FALSE));
1505: /* Creates parallel work vectors used in presolve */
1506: if (!pcbddc->original_rhs) PetscCall(VecDuplicate(pcis->vec1_global, &pcbddc->original_rhs));
1507: if (!pcbddc->temp_solution) PetscCall(VecDuplicate(pcis->vec1_global, &pcbddc->temp_solution));
1509: pcbddc->temp_solution_used = PETSC_FALSE;
1510: if (x) {
1511: PetscCall(PetscObjectReference((PetscObject)x));
1512: used_vec = x;
1513: } else { /* it can only happen when calling PCBDDCMatFETIDPGetRHS */
1514: PetscCall(PetscObjectReference((PetscObject)pcbddc->temp_solution));
1515: used_vec = pcbddc->temp_solution;
1516: PetscCall(VecSet(used_vec, 0.0));
1517: pcbddc->temp_solution_used = PETSC_TRUE;
1518: PetscCall(VecCopy(rhs, pcbddc->original_rhs));
1519: save_rhs = PETSC_FALSE;
1520: pcbddc->eliminate_dirdofs = PETSC_TRUE;
1521: }
1523: /* hack into ksp data structure since PCPreSolve comes earlier than setting to zero the guess in src/ksp/ksp/interface/itfunc.c */
1524: if (ksp) {
1525: /* store the flag for the initial guess since it will be restored back during PCPostSolve_BDDC */
1526: PetscCall(KSPGetInitialGuessNonzero(ksp, &pcbddc->ksp_guess_nonzero));
1527: if (!pcbddc->ksp_guess_nonzero) PetscCall(VecSet(used_vec, 0.0));
1528: }
1530: pcbddc->rhs_change = PETSC_FALSE;
1531: /* Take into account zeroed rows -> change rhs and store solution removed */
1532: if (rhs && pcbddc->eliminate_dirdofs) {
1533: IS dirIS = NULL;
1535: /* DirichletBoundariesLocal may not be consistent among neighbours; gets a dirichlet dofs IS from graph (may be cached) */
1536: PetscCall(PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph, &dirIS));
1537: if (dirIS) {
1538: Mat_IS *matis = (Mat_IS *)pc->pmat->data;
1539: PetscInt dirsize, i, *is_indices;
1540: PetscScalar *array_x;
1541: const PetscScalar *array_diagonal;
1543: PetscCall(MatGetDiagonal(pc->pmat, pcis->vec1_global));
1544: PetscCall(VecPointwiseDivide(pcis->vec1_global, rhs, pcis->vec1_global));
1545: PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_global, pcis->vec2_N, INSERT_VALUES, SCATTER_FORWARD));
1546: PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_global, pcis->vec2_N, INSERT_VALUES, SCATTER_FORWARD));
1547: PetscCall(VecScatterBegin(matis->rctx, used_vec, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD));
1548: PetscCall(VecScatterEnd(matis->rctx, used_vec, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD));
1549: PetscCall(ISGetLocalSize(dirIS, &dirsize));
1550: PetscCall(VecGetArray(pcis->vec1_N, &array_x));
1551: PetscCall(VecGetArrayRead(pcis->vec2_N, &array_diagonal));
1552: PetscCall(ISGetIndices(dirIS, (const PetscInt **)&is_indices));
1553: for (i = 0; i < dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]];
1554: PetscCall(ISRestoreIndices(dirIS, (const PetscInt **)&is_indices));
1555: PetscCall(VecRestoreArrayRead(pcis->vec2_N, &array_diagonal));
1556: PetscCall(VecRestoreArray(pcis->vec1_N, &array_x));
1557: PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_N, used_vec, INSERT_VALUES, SCATTER_REVERSE));
1558: PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_N, used_vec, INSERT_VALUES, SCATTER_REVERSE));
1559: pcbddc->rhs_change = PETSC_TRUE;
1560: PetscCall(ISDestroy(&dirIS));
1561: }
1562: }
1564: /* remove the computed solution or the initial guess from the rhs */
1565: if (pcbddc->rhs_change || (ksp && pcbddc->ksp_guess_nonzero)) {
1566: /* save the original rhs */
1567: if (save_rhs) {
1568: PetscCall(VecSwap(rhs, pcbddc->original_rhs));
1569: save_rhs = PETSC_FALSE;
1570: }
1571: pcbddc->rhs_change = PETSC_TRUE;
1572: PetscCall(VecScale(used_vec, -1.0));
1573: PetscCall(MatMultAdd(pc->mat, used_vec, pcbddc->original_rhs, rhs));
1574: PetscCall(VecScale(used_vec, -1.0));
1575: PetscCall(VecCopy(used_vec, pcbddc->temp_solution));
1576: pcbddc->temp_solution_used = PETSC_TRUE;
1577: if (ksp) PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_FALSE));
1578: }
1579: PetscCall(VecDestroy(&used_vec));
1581: /* compute initial vector in benign space if needed
1582: and remove non-benign solution from the rhs */
1583: benign_correction_computed = PETSC_FALSE;
1584: if (rhs && pcbddc->benign_compute_correction && (pcbddc->benign_have_null || pcbddc->benign_apply_coarse_only)) {
1585: /* compute u^*_h using ideas similar to those in Xuemin Tu's PhD thesis (see Section 4.8.1)
1586: Recursively apply BDDC in the multilevel case */
1587: if (!pcbddc->benign_vec) PetscCall(VecDuplicate(rhs, &pcbddc->benign_vec));
1588: /* keep applying coarse solver unless we no longer have benign subdomains */
1589: pcbddc->benign_apply_coarse_only = pcbddc->benign_have_null ? PETSC_TRUE : PETSC_FALSE;
1590: if (!pcbddc->benign_skip_correction) {
1591: PetscCall(PCApply_BDDC(pc, rhs, pcbddc->benign_vec));
1592: benign_correction_computed = PETSC_TRUE;
1593: if (pcbddc->temp_solution_used) PetscCall(VecAXPY(pcbddc->temp_solution, 1.0, pcbddc->benign_vec));
1594: PetscCall(VecScale(pcbddc->benign_vec, -1.0));
1595: /* store the original rhs if not done earlier */
1596: if (save_rhs) PetscCall(VecSwap(rhs, pcbddc->original_rhs));
1597: if (pcbddc->rhs_change) {
1598: PetscCall(MatMultAdd(pc->mat, pcbddc->benign_vec, rhs, rhs));
1599: } else {
1600: PetscCall(MatMultAdd(pc->mat, pcbddc->benign_vec, pcbddc->original_rhs, rhs));
1601: }
1602: pcbddc->rhs_change = PETSC_TRUE;
1603: }
1604: pcbddc->benign_apply_coarse_only = PETSC_FALSE;
1605: } else {
1606: PetscCall(VecDestroy(&pcbddc->benign_vec));
1607: }
1609: /* dbg output */
1610: if (pcbddc->dbg_flag && benign_correction_computed) {
1611: Vec v;
1613: PetscCall(VecDuplicate(pcis->vec1_global, &v));
1614: if (pcbddc->ChangeOfBasisMatrix) {
1615: PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, rhs, v));
1616: } else {
1617: PetscCall(VecCopy(rhs, v));
1618: }
1619: PetscCall(PCBDDCBenignGetOrSetP0(pc, v, PETSC_TRUE));
1620: PetscCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "LEVEL %" PetscInt_FMT ": is the correction benign?\n", pcbddc->current_level));
1621: PetscCall(PetscScalarView(pcbddc->benign_n, pcbddc->benign_p0, pcbddc->dbg_viewer));
1622: PetscCall(PetscViewerFlush(pcbddc->dbg_viewer));
1623: PetscCall(VecDestroy(&v));
1624: }
1626: /* set initial guess if using PCG */
1627: pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
1628: if (x && pcbddc->use_exact_dirichlet_trick) {
1629: PetscCall(VecSet(x, 0.0));
1630: if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) {
1631: if (benign_correction_computed) { /* we have already saved the changed rhs */
1632: PetscCall(VecLockReadPop(pcis->vec1_global));
1633: } else {
1634: PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, rhs, pcis->vec1_global));
1635: }
1636: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec1_global, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1637: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec1_global, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1638: } else {
1639: PetscCall(VecScatterBegin(pcis->global_to_D, rhs, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1640: PetscCall(VecScatterEnd(pcis->global_to_D, rhs, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1641: }
1642: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1643: PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D));
1644: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1645: PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D));
1646: if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) {
1647: PetscCall(VecSet(pcis->vec1_global, 0.));
1648: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, pcis->vec1_global, INSERT_VALUES, SCATTER_REVERSE));
1649: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, pcis->vec1_global, INSERT_VALUES, SCATTER_REVERSE));
1650: PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcis->vec1_global, x));
1651: } else {
1652: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, x, INSERT_VALUES, SCATTER_REVERSE));
1653: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, x, INSERT_VALUES, SCATTER_REVERSE));
1654: }
1655: if (ksp) PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
1656: pcbddc->exact_dirichlet_trick_app = PETSC_TRUE;
1657: } else if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior && benign_correction_computed && pcbddc->use_exact_dirichlet_trick) {
1658: PetscCall(VecLockReadPop(pcis->vec1_global));
1659: }
1660: PetscFunctionReturn(PETSC_SUCCESS);
1661: }
1663: static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
1664: {
1665: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1667: PetscFunctionBegin;
1668: /* add solution removed in presolve */
1669: if (x && pcbddc->rhs_change) {
1670: if (pcbddc->temp_solution_used) PetscCall(VecAXPY(x, 1.0, pcbddc->temp_solution));
1671: else if (pcbddc->benign_compute_correction && pcbddc->benign_vec) PetscCall(VecAXPY(x, -1.0, pcbddc->benign_vec));
1672: /* restore to original state (not for FETI-DP) */
1673: if (ksp) pcbddc->temp_solution_used = PETSC_FALSE;
1674: }
1676: /* restore rhs to its original state (not needed for FETI-DP) */
1677: if (rhs && pcbddc->rhs_change) {
1678: PetscCall(VecSwap(rhs, pcbddc->original_rhs));
1679: pcbddc->rhs_change = PETSC_FALSE;
1680: }
1681: /* restore ksp guess state */
1682: if (ksp) {
1683: PetscCall(KSPSetInitialGuessNonzero(ksp, pcbddc->ksp_guess_nonzero));
1684: /* reset flag for exact dirichlet trick */
1685: pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
1686: }
1687: PetscFunctionReturn(PETSC_SUCCESS);
1688: }
1690: static PetscErrorCode PCSetUp_BDDC(PC pc)
1691: {
1692: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1693: PCBDDCSubSchurs sub_schurs;
1694: Mat_IS *matis;
1695: MatNullSpace nearnullspace;
1696: Mat lA;
1697: IS lP, zerodiag = NULL;
1698: PetscInt nrows, ncols;
1699: PetscMPIInt size;
1700: PetscBool computesubschurs;
1701: PetscBool computeconstraintsmatrix;
1702: PetscBool new_nearnullspace_provided, ismatis;
1703: PetscBool isset, issym, isspd;
1705: PetscFunctionBegin;
1706: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATIS, &ismatis));
1707: PetscCheck(ismatis, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PCBDDC preconditioner requires matrix of type MATIS");
1708: PetscCall(MatGetSize(pc->pmat, &nrows, &ncols));
1709: PetscCheck(nrows == ncols, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCBDDC preconditioner requires a square matrix for constructing the preconditioner");
1710: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
1712: matis = (Mat_IS *)pc->pmat->data;
1713: /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */
1714: /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetUp
1715: Also, BDDC builds its own KSP for the Dirichlet problem */
1716: if (!pc->setupcalled || pc->flag == DIFFERENT_NONZERO_PATTERN) pcbddc->recompute_topography = PETSC_TRUE;
1717: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &pcbddc->recompute_topography, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)pc)));
1718: if (pcbddc->recompute_topography) {
1719: pcbddc->graphanalyzed = PETSC_FALSE;
1720: computeconstraintsmatrix = PETSC_TRUE;
1721: } else {
1722: computeconstraintsmatrix = PETSC_FALSE;
1723: }
1725: /* check parameters' compatibility */
1726: if (!pcbddc->use_deluxe_scaling) pcbddc->deluxe_zerorows = PETSC_FALSE;
1727: pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_threshold[0] != 0.0 || pcbddc->adaptive_threshold[1] != 0.0);
1728: pcbddc->use_deluxe_scaling = (PetscBool)(pcbddc->use_deluxe_scaling && (size > 1 || matis->allow_repeated));
1729: pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_selection && (size > 1 || matis->allow_repeated));
1730: pcbddc->adaptive_userdefined = (PetscBool)(pcbddc->adaptive_selection && pcbddc->adaptive_userdefined);
1731: if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE;
1733: computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling);
1735: /* activate all connected components if the netflux has been requested */
1736: if (pcbddc->compute_nonetflux) {
1737: pcbddc->use_vertices = PETSC_TRUE;
1738: pcbddc->use_edges = PETSC_TRUE;
1739: pcbddc->use_faces = PETSC_TRUE;
1740: }
1742: /* Get stdout for dbg */
1743: if (pcbddc->dbg_flag) {
1744: if (!pcbddc->dbg_viewer) pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc));
1745: PetscCall(PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer));
1746: PetscCall(PetscViewerASCIIAddTab(pcbddc->dbg_viewer, 2 * pcbddc->current_level));
1747: }
1749: /* process topology information */
1750: PetscCall(PetscLogEventBegin(PC_BDDC_Topology[pcbddc->current_level], pc, 0, 0, 0));
1751: if (pcbddc->recompute_topography) {
1752: PetscCall(PCBDDCComputeLocalTopologyInfo(pc));
1753: if (pcbddc->discretegradient) PetscCall(PCBDDCNedelecSupport(pc));
1754: }
1755: if (pcbddc->corner_selected) pcbddc->use_vertices = PETSC_TRUE;
1757: /* change basis if requested by the user */
1758: if (pcbddc->user_ChangeOfBasisMatrix) {
1759: /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */
1760: pcbddc->use_change_of_basis = PETSC_FALSE;
1761: PetscCall(PCBDDCComputeLocalMatrix(pc, pcbddc->user_ChangeOfBasisMatrix));
1762: } else {
1763: PetscCall(MatDestroy(&pcbddc->local_mat));
1764: PetscCall(PetscObjectReference((PetscObject)matis->A));
1765: pcbddc->local_mat = matis->A;
1766: }
1768: /*
1769: Compute change of basis on local pressures (aka zerodiag dofs) with the benign trick
1770: This should come earlier than PCISSetUp for extracting the correct subdomain matrices
1771: */
1772: PetscCall(PCBDDCBenignShellMat(pc, PETSC_TRUE));
1773: if (pcbddc->benign_saddle_point) {
1774: PC_IS *pcis = (PC_IS *)pc->data;
1776: if (pcbddc->user_ChangeOfBasisMatrix || pcbddc->use_change_of_basis || !computesubschurs) pcbddc->benign_change_explicit = PETSC_TRUE;
1777: /* detect local saddle point and change the basis in pcbddc->local_mat */
1778: PetscCall(PCBDDCBenignDetectSaddlePoint(pc, (PetscBool)(!pcbddc->recompute_topography), &zerodiag));
1779: /* pop B0 mat from local mat */
1780: PetscCall(PCBDDCBenignPopOrPushB0(pc, PETSC_TRUE));
1781: /* give pcis a hint to not reuse submatrices during PCISCreate */
1782: if (pc->flag == SAME_NONZERO_PATTERN && pcis->reusesubmatrices == PETSC_TRUE) {
1783: if (pcbddc->benign_n && (pcbddc->benign_change_explicit || pcbddc->dbg_flag)) {
1784: pcis->reusesubmatrices = PETSC_FALSE;
1785: } else {
1786: pcis->reusesubmatrices = PETSC_TRUE;
1787: }
1788: } else {
1789: pcis->reusesubmatrices = PETSC_FALSE;
1790: }
1791: }
1793: /* propagate relevant information */
1794: PetscCall(MatIsSymmetricKnown(matis->A, &isset, &issym));
1795: if (isset) PetscCall(MatSetOption(pcbddc->local_mat, MAT_SYMMETRIC, issym));
1796: PetscCall(MatIsSPDKnown(matis->A, &isset, &isspd));
1797: if (isset) PetscCall(MatSetOption(pcbddc->local_mat, MAT_SPD, isspd));
1799: /* Set up all the "iterative substructuring" common block without computing solvers */
1800: {
1801: Mat temp_mat;
1803: temp_mat = matis->A;
1804: matis->A = pcbddc->local_mat;
1805: PetscCall(PCISSetUp(pc, PETSC_TRUE, PETSC_FALSE));
1806: pcbddc->local_mat = matis->A;
1807: matis->A = temp_mat;
1808: }
1810: /* Analyze interface */
1811: if (!pcbddc->graphanalyzed) {
1812: PetscCall(PCBDDCAnalyzeInterface(pc));
1813: computeconstraintsmatrix = PETSC_TRUE;
1814: 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");
1815: if (pcbddc->compute_nonetflux) {
1816: MatNullSpace nnfnnsp;
1818: PetscCheck(pcbddc->divudotp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Missing divudotp operator");
1819: PetscCall(PCBDDCComputeNoNetFlux(pc->pmat, pcbddc->divudotp, pcbddc->divudotp_trans, pcbddc->divudotp_vl2l, pcbddc->mat_graph, &nnfnnsp));
1820: /* TODO what if a nearnullspace is already attached? */
1821: if (nnfnnsp) {
1822: PetscCall(MatSetNearNullSpace(pc->pmat, nnfnnsp));
1823: PetscCall(MatNullSpaceDestroy(&nnfnnsp));
1824: }
1825: }
1826: }
1827: PetscCall(PetscLogEventEnd(PC_BDDC_Topology[pcbddc->current_level], pc, 0, 0, 0));
1829: /* check existence of a divergence free extension, i.e.
1830: b(v_I,p_0) = 0 for all v_I (raise error if not).
1831: Also, check that PCBDDCBenignGetOrSetP0 works */
1832: if (pcbddc->benign_saddle_point && pcbddc->dbg_flag > 1) PetscCall(PCBDDCBenignCheck(pc, zerodiag));
1833: PetscCall(ISDestroy(&zerodiag));
1835: /* Setup local dirichlet solver ksp_D and sub_schurs solvers */
1836: if (computesubschurs && pcbddc->recompute_topography) PetscCall(PCBDDCInitSubSchurs(pc));
1837: /* SetUp Scaling operator (scaling matrices could be needed in SubSchursSetUp)*/
1838: if (!pcbddc->use_deluxe_scaling) PetscCall(PCBDDCScalingSetUp(pc));
1840: /* finish setup solvers and do adaptive selection of constraints */
1841: sub_schurs = pcbddc->sub_schurs;
1842: if (sub_schurs && sub_schurs->schur_explicit) {
1843: if (computesubschurs) PetscCall(PCBDDCSetUpSubSchurs(pc));
1844: PetscCall(PCBDDCSetUpLocalSolvers(pc, PETSC_TRUE, PETSC_FALSE));
1845: } else {
1846: PetscCall(PCBDDCSetUpLocalSolvers(pc, PETSC_TRUE, PETSC_FALSE));
1847: if (computesubschurs) PetscCall(PCBDDCSetUpSubSchurs(pc));
1848: }
1849: if (pcbddc->adaptive_selection) {
1850: PetscCall(PCBDDCAdaptiveSelection(pc));
1851: computeconstraintsmatrix = PETSC_TRUE;
1852: }
1854: /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */
1855: new_nearnullspace_provided = PETSC_FALSE;
1856: PetscCall(MatGetNearNullSpace(pc->pmat, &nearnullspace));
1857: if (pcbddc->onearnullspace) { /* already used nearnullspace */
1858: if (!nearnullspace) { /* near null space attached to mat has been destroyed */
1859: new_nearnullspace_provided = PETSC_TRUE;
1860: } else {
1861: /* determine if the two nullspaces are different (should be lightweight) */
1862: if (nearnullspace != pcbddc->onearnullspace) {
1863: new_nearnullspace_provided = PETSC_TRUE;
1864: } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */
1865: const Vec *nearnullvecs;
1866: PetscObjectState state;
1867: PetscInt nnsp_size;
1868: PetscCall(MatNullSpaceGetVecs(nearnullspace, NULL, &nnsp_size, &nearnullvecs));
1869: for (PetscInt i = 0; i < nnsp_size; i++) {
1870: PetscCall(PetscObjectStateGet((PetscObject)nearnullvecs[i], &state));
1871: if (pcbddc->onearnullvecs_state[i] != state) {
1872: new_nearnullspace_provided = PETSC_TRUE;
1873: break;
1874: }
1875: }
1876: }
1877: }
1878: } else {
1879: if (!nearnullspace) { /* both nearnullspaces are null */
1880: new_nearnullspace_provided = PETSC_FALSE;
1881: } else { /* nearnullspace attached later */
1882: new_nearnullspace_provided = PETSC_TRUE;
1883: }
1884: }
1886: /* Setup constraints and related work vectors */
1887: /* reset primal space flags */
1888: PetscCall(PetscLogEventBegin(PC_BDDC_LocalWork[pcbddc->current_level], pc, 0, 0, 0));
1889: pcbddc->new_primal_space = PETSC_FALSE;
1890: pcbddc->new_primal_space_local = PETSC_FALSE;
1891: if (computeconstraintsmatrix || new_nearnullspace_provided) {
1892: /* It also sets the primal space flags */
1893: PetscCall(PCBDDCConstraintsSetUp(pc));
1894: }
1895: /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */
1896: PetscCall(PCBDDCSetUpLocalWorkVectors(pc));
1898: if (pcbddc->use_change_of_basis) {
1899: PC_IS *pcis = (PC_IS *)pc->data;
1901: PetscCall(PCBDDCComputeLocalMatrix(pc, pcbddc->ChangeOfBasisMatrix));
1902: if (pcbddc->benign_change) {
1903: PetscCall(MatDestroy(&pcbddc->benign_B0));
1904: /* pop B0 from pcbddc->local_mat */
1905: PetscCall(PCBDDCBenignPopOrPushB0(pc, PETSC_TRUE));
1906: }
1907: /* get submatrices */
1908: PetscCall(MatDestroy(&pcis->A_IB));
1909: PetscCall(MatDestroy(&pcis->A_BI));
1910: PetscCall(MatDestroy(&pcis->A_BB));
1911: PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcis->is_B_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &pcis->A_BB));
1912: PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcis->is_I_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &pcis->A_IB));
1913: PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcis->is_B_local, pcis->is_I_local, MAT_INITIAL_MATRIX, &pcis->A_BI));
1914: /* set flag in pcis to not reuse submatrices during PCISCreate */
1915: pcis->reusesubmatrices = PETSC_FALSE;
1916: } else if (!pcbddc->user_ChangeOfBasisMatrix && !pcbddc->benign_change) {
1917: PetscCall(MatDestroy(&pcbddc->local_mat));
1918: PetscCall(PetscObjectReference((PetscObject)matis->A));
1919: pcbddc->local_mat = matis->A;
1920: }
1922: /* interface pressure block row for B_C */
1923: PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_lP", (PetscObject *)&lP));
1924: PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_lA", (PetscObject *)&lA));
1925: if (lA && lP) {
1926: PC_IS *pcis = (PC_IS *)pc->data;
1927: Mat B_BI, B_BB, Bt_BI, Bt_BB;
1928: PetscBool issym;
1930: PetscCall(MatIsSymmetric(lA, PETSC_SMALL, &issym));
1931: if (issym) {
1932: PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_I_local, MAT_INITIAL_MATRIX, &B_BI));
1933: PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_B_local, MAT_INITIAL_MATRIX, &B_BB));
1934: PetscCall(MatCreateTranspose(B_BI, &Bt_BI));
1935: PetscCall(MatCreateTranspose(B_BB, &Bt_BB));
1936: } else {
1937: PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_I_local, MAT_INITIAL_MATRIX, &B_BI));
1938: PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_B_local, MAT_INITIAL_MATRIX, &B_BB));
1939: PetscCall(MatCreateSubMatrix(lA, pcis->is_I_local, lP, MAT_INITIAL_MATRIX, &Bt_BI));
1940: PetscCall(MatCreateSubMatrix(lA, pcis->is_B_local, lP, MAT_INITIAL_MATRIX, &Bt_BB));
1941: }
1942: PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_B_BI", (PetscObject)B_BI));
1943: PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_B_BB", (PetscObject)B_BB));
1944: PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_Bt_BI", (PetscObject)Bt_BI));
1945: PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_Bt_BB", (PetscObject)Bt_BB));
1946: PetscCall(MatDestroy(&B_BI));
1947: PetscCall(MatDestroy(&B_BB));
1948: PetscCall(MatDestroy(&Bt_BI));
1949: PetscCall(MatDestroy(&Bt_BB));
1950: }
1951: PetscCall(PetscLogEventEnd(PC_BDDC_LocalWork[pcbddc->current_level], pc, 0, 0, 0));
1953: /* SetUp coarse and local Neumann solvers */
1954: PetscCall(PCBDDCSetUpSolvers(pc));
1955: /* SetUp Scaling operator */
1956: if (pcbddc->use_deluxe_scaling) PetscCall(PCBDDCScalingSetUp(pc));
1958: /* mark topography as done */
1959: pcbddc->recompute_topography = PETSC_FALSE;
1961: /* wrap pcis->A_IB and pcis->A_BI if we did not change explicitly the variables on the pressures */
1962: PetscCall(PCBDDCBenignShellMat(pc, PETSC_FALSE));
1964: if (pcbddc->dbg_flag) {
1965: PetscCall(PetscViewerASCIISubtractTab(pcbddc->dbg_viewer, 2 * pcbddc->current_level));
1966: PetscCall(PetscViewerASCIIPopSynchronized(pcbddc->dbg_viewer));
1967: }
1969: { /* Dump customization */
1970: PetscInt save_version = PETSC_DECIDE;
1971: PetscBool flg;
1972: char save[PETSC_MAX_PATH_LEN] = {'\0'};
1974: PetscCall(PetscOptionsGetString(NULL, ((PetscObject)pc)->prefix, "-pc_bddc_save", save, sizeof(save), &flg));
1975: PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)pc)->prefix, "-pc_bddc_save_version", &save_version, NULL));
1976: if (flg) {
1977: size_t len;
1979: PetscCall(PetscStrlen(save, &len));
1980: PetscCall(PCBDDCSaveCustomization(pc, len ? save : NULL, save_version));
1981: }
1982: }
1983: PetscFunctionReturn(PETSC_SUCCESS);
1984: }
1986: static PetscErrorCode PCApply_BDDC(PC pc, Vec r, Vec z)
1987: {
1988: PC_IS *pcis = (PC_IS *)pc->data;
1989: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1990: Mat lA = NULL;
1991: PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B;
1992: const PetscScalar one = 1.0;
1993: const PetscScalar m_one = -1.0;
1994: const PetscScalar zero = 0.0;
1995: /* This code is similar to that provided in nn.c for PCNN
1996: NN interface preconditioner changed to BDDC
1997: Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */
1999: PetscFunctionBegin;
2000: PetscCall(PetscCitationsRegister(citation, &cited));
2001: if (pcbddc->switch_static) PetscCall(MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat, &lA));
2003: if (pcbddc->ChangeOfBasisMatrix) {
2004: Vec swap;
2006: PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, r, pcbddc->work_change));
2007: swap = pcbddc->work_change;
2008: pcbddc->work_change = r;
2009: r = swap;
2010: /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */
2011: if (pcbddc->benign_apply_coarse_only && pcbddc->use_exact_dirichlet_trick && pcbddc->change_interior) {
2012: PetscCall(VecCopy(r, pcis->vec1_global));
2013: PetscCall(VecLockReadPush(pcis->vec1_global));
2014: }
2015: }
2016: if (pcbddc->benign_have_null) { /* get p0 from r */
2017: PetscCall(PCBDDCBenignGetOrSetP0(pc, r, PETSC_TRUE));
2018: }
2019: if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_DIRICHLET && !pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
2020: PetscCall(VecCopy(r, z));
2021: /* First Dirichlet solve */
2022: PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
2023: PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
2024: /*
2025: Assembling right-hand side for BDDC operator
2026: - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE)
2027: - pcis->vec1_B the interface part of the global vector z
2028: */
2029: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2030: if (n_D) {
2031: PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D));
2032: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2033: PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D));
2034: PetscCall(VecScale(pcis->vec2_D, m_one));
2035: if (pcbddc->switch_static) {
2036: PetscCall(VecSet(pcis->vec1_N, 0.));
2037: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2038: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2039: if (!pcbddc->switch_static_change) PetscCall(MatMult(lA, pcis->vec1_N, pcis->vec2_N));
2040: else {
2041: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2042: PetscCall(MatMult(lA, pcis->vec2_N, pcis->vec1_N));
2043: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2044: }
2045: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
2046: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
2047: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2048: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2049: } else {
2050: PetscCall(MatMult(pcis->A_BI, pcis->vec2_D, pcis->vec1_B));
2051: }
2052: } else {
2053: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2054: PetscCall(VecSet(pcis->vec1_B, zero));
2055: }
2056: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
2057: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
2058: PetscCall(PCBDDCScalingRestriction(pc, z, pcis->vec1_B));
2059: } else {
2060: if (!pcbddc->benign_apply_coarse_only) PetscCall(PCBDDCScalingRestriction(pc, r, pcis->vec1_B));
2061: }
2062: if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_LUMP) {
2063: PetscCheck(pcbddc->switch_static, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "You forgot to pass -pc_bddc_switch_static");
2064: PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
2065: PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
2066: }
2068: /* Apply interface preconditioner
2069: input/output vecs: pcis->vec1_B and pcis->vec1_D */
2070: PetscCall(PCBDDCApplyInterfacePreconditioner(pc, PETSC_FALSE));
2072: /* Apply transpose of partition of unity operator */
2073: PetscCall(PCBDDCScalingExtension(pc, pcis->vec1_B, z));
2074: if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_LUMP) {
2075: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec1_D, z, INSERT_VALUES, SCATTER_REVERSE));
2076: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec1_D, z, INSERT_VALUES, SCATTER_REVERSE));
2077: PetscFunctionReturn(PETSC_SUCCESS);
2078: }
2079: /* Second Dirichlet solve and assembling of output */
2080: PetscCall(VecScatterBegin(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2081: PetscCall(VecScatterEnd(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2082: if (n_B) {
2083: if (pcbddc->switch_static) {
2084: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2085: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2086: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2087: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2088: if (!pcbddc->switch_static_change) PetscCall(MatMult(lA, pcis->vec1_N, pcis->vec2_N));
2089: else {
2090: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2091: PetscCall(MatMult(lA, pcis->vec2_N, pcis->vec1_N));
2092: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2093: }
2094: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
2095: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
2096: } else {
2097: PetscCall(MatMult(pcis->A_IB, pcis->vec1_B, pcis->vec3_D));
2098: }
2099: } else if (pcbddc->switch_static) { /* n_B is zero */
2100: if (!pcbddc->switch_static_change) PetscCall(MatMult(lA, pcis->vec1_D, pcis->vec3_D));
2101: else {
2102: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_D, pcis->vec1_N));
2103: PetscCall(MatMult(lA, pcis->vec1_N, pcis->vec2_N));
2104: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec2_N, pcis->vec3_D));
2105: }
2106: }
2107: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2108: PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec3_D, pcis->vec4_D));
2109: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2110: PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec4_D));
2112: if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
2113: if (pcbddc->switch_static) PetscCall(VecAXPBYPCZ(pcis->vec2_D, m_one, one, m_one, pcis->vec4_D, pcis->vec1_D));
2114: else PetscCall(VecAXPBY(pcis->vec2_D, m_one, m_one, pcis->vec4_D));
2115: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
2116: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
2117: } else {
2118: if (pcbddc->switch_static) PetscCall(VecAXPBY(pcis->vec4_D, one, m_one, pcis->vec1_D));
2119: else PetscCall(VecScale(pcis->vec4_D, m_one));
2120: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
2121: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
2122: }
2123: if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */
2124: if (pcbddc->benign_apply_coarse_only) PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
2125: PetscCall(PCBDDCBenignGetOrSetP0(pc, z, PETSC_FALSE));
2126: }
2128: if (pcbddc->ChangeOfBasisMatrix) {
2129: pcbddc->work_change = r;
2130: PetscCall(VecCopy(z, pcbddc->work_change));
2131: PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcbddc->work_change, z));
2132: }
2133: PetscFunctionReturn(PETSC_SUCCESS);
2134: }
2136: static PetscErrorCode PCApplyTranspose_BDDC(PC pc, Vec r, Vec z)
2137: {
2138: PC_IS *pcis = (PC_IS *)pc->data;
2139: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
2140: Mat lA = NULL;
2141: PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B;
2142: const PetscScalar one = 1.0;
2143: const PetscScalar m_one = -1.0;
2144: const PetscScalar zero = 0.0;
2146: PetscFunctionBegin;
2147: PetscCall(PetscCitationsRegister(citation, &cited));
2148: if (pcbddc->switch_static) PetscCall(MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat, &lA));
2149: if (pcbddc->ChangeOfBasisMatrix) {
2150: Vec swap;
2152: PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, r, pcbddc->work_change));
2153: swap = pcbddc->work_change;
2154: pcbddc->work_change = r;
2155: r = swap;
2156: /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */
2157: if (pcbddc->benign_apply_coarse_only && pcbddc->exact_dirichlet_trick_app && pcbddc->change_interior) {
2158: PetscCall(VecCopy(r, pcis->vec1_global));
2159: PetscCall(VecLockReadPush(pcis->vec1_global));
2160: }
2161: }
2162: if (pcbddc->benign_have_null) { /* get p0 from r */
2163: PetscCall(PCBDDCBenignGetOrSetP0(pc, r, PETSC_TRUE));
2164: }
2165: if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
2166: PetscCall(VecCopy(r, z));
2167: /* First Dirichlet solve */
2168: PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
2169: PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
2170: /*
2171: Assembling right-hand side for BDDC operator
2172: - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE)
2173: - pcis->vec1_B the interface part of the global vector z
2174: */
2175: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2176: if (n_D) {
2177: PetscCall(KSPSolveTranspose(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D));
2178: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2179: PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D));
2180: PetscCall(VecScale(pcis->vec2_D, m_one));
2181: if (pcbddc->switch_static) {
2182: PetscCall(VecSet(pcis->vec1_N, 0.));
2183: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2184: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2185: if (!pcbddc->switch_static_change) {
2186: PetscCall(MatMultTranspose(lA, pcis->vec1_N, pcis->vec2_N));
2187: } else {
2188: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2189: PetscCall(MatMultTranspose(lA, pcis->vec2_N, pcis->vec1_N));
2190: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2191: }
2192: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
2193: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
2194: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2195: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2196: } else {
2197: PetscCall(MatMultTranspose(pcis->A_IB, pcis->vec2_D, pcis->vec1_B));
2198: }
2199: } else {
2200: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2201: PetscCall(VecSet(pcis->vec1_B, zero));
2202: }
2203: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
2204: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
2205: PetscCall(PCBDDCScalingRestriction(pc, z, pcis->vec1_B));
2206: } else {
2207: PetscCall(PCBDDCScalingRestriction(pc, r, pcis->vec1_B));
2208: }
2210: /* Apply interface preconditioner
2211: input/output vecs: pcis->vec1_B and pcis->vec1_D */
2212: PetscCall(PCBDDCApplyInterfacePreconditioner(pc, PETSC_TRUE));
2214: /* Apply transpose of partition of unity operator */
2215: PetscCall(PCBDDCScalingExtension(pc, pcis->vec1_B, z));
2217: /* Second Dirichlet solve and assembling of output */
2218: PetscCall(VecScatterBegin(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2219: PetscCall(VecScatterEnd(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2220: if (n_B) {
2221: if (pcbddc->switch_static) {
2222: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2223: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2224: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2225: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2226: if (!pcbddc->switch_static_change) {
2227: PetscCall(MatMultTranspose(lA, pcis->vec1_N, pcis->vec2_N));
2228: } else {
2229: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2230: PetscCall(MatMultTranspose(lA, pcis->vec2_N, pcis->vec1_N));
2231: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2232: }
2233: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
2234: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
2235: } else {
2236: PetscCall(MatMultTranspose(pcis->A_BI, pcis->vec1_B, pcis->vec3_D));
2237: }
2238: } else if (pcbddc->switch_static) { /* n_B is zero */
2239: if (!pcbddc->switch_static_change) {
2240: PetscCall(MatMultTranspose(lA, pcis->vec1_D, pcis->vec3_D));
2241: } else {
2242: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_D, pcis->vec1_N));
2243: PetscCall(MatMultTranspose(lA, pcis->vec1_N, pcis->vec2_N));
2244: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec2_N, pcis->vec3_D));
2245: }
2246: }
2247: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2248: PetscCall(KSPSolveTranspose(pcbddc->ksp_D, pcis->vec3_D, pcis->vec4_D));
2249: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2250: PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec4_D));
2251: if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
2252: if (pcbddc->switch_static) PetscCall(VecAXPBYPCZ(pcis->vec2_D, m_one, one, m_one, pcis->vec4_D, pcis->vec1_D));
2253: else PetscCall(VecAXPBY(pcis->vec2_D, m_one, m_one, pcis->vec4_D));
2254: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
2255: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
2256: } else {
2257: if (pcbddc->switch_static) PetscCall(VecAXPBY(pcis->vec4_D, one, m_one, pcis->vec1_D));
2258: else PetscCall(VecScale(pcis->vec4_D, m_one));
2259: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
2260: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
2261: }
2262: if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */
2263: PetscCall(PCBDDCBenignGetOrSetP0(pc, z, PETSC_FALSE));
2264: }
2265: if (pcbddc->ChangeOfBasisMatrix) {
2266: pcbddc->work_change = r;
2267: PetscCall(VecCopy(z, pcbddc->work_change));
2268: PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcbddc->work_change, z));
2269: }
2270: PetscFunctionReturn(PETSC_SUCCESS);
2271: }
2273: static PetscErrorCode PCReset_BDDC(PC pc)
2274: {
2275: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
2276: PC_IS *pcis = (PC_IS *)pc->data;
2277: KSP kspD, kspR, kspC;
2279: PetscFunctionBegin;
2280: /* free BDDC custom data */
2281: PetscCall(PCBDDCResetCustomization(pc));
2282: /* destroy objects related to topography */
2283: PetscCall(PCBDDCResetTopography(pc));
2284: /* destroy objects for scaling operator */
2285: PetscCall(PCBDDCScalingDestroy(pc));
2286: /* free solvers stuff */
2287: PetscCall(PCBDDCResetSolvers(pc));
2288: /* free global vectors needed in presolve */
2289: PetscCall(VecDestroy(&pcbddc->temp_solution));
2290: PetscCall(VecDestroy(&pcbddc->original_rhs));
2291: /* free data created by PCIS */
2292: PetscCall(PCISReset(pc));
2294: /* restore defaults */
2295: kspD = pcbddc->ksp_D;
2296: kspR = pcbddc->ksp_R;
2297: kspC = pcbddc->coarse_ksp;
2298: PetscCall(PetscMemzero(pc->data, sizeof(*pcbddc)));
2299: pcis->n_neigh = -1;
2300: pcis->scaling_factor = 1.0;
2301: pcis->reusesubmatrices = PETSC_TRUE;
2302: pcbddc->use_local_adj = PETSC_TRUE;
2303: pcbddc->use_vertices = PETSC_TRUE;
2304: pcbddc->use_edges = PETSC_TRUE;
2305: pcbddc->symmetric_primal = PETSC_TRUE;
2306: pcbddc->vertex_size = 1;
2307: pcbddc->recompute_topography = PETSC_TRUE;
2308: pcbddc->coarse_size = -1;
2309: pcbddc->use_exact_dirichlet_trick = PETSC_TRUE;
2310: pcbddc->coarsening_ratio = 8;
2311: pcbddc->coarse_eqs_per_proc = 1;
2312: pcbddc->benign_compute_correction = PETSC_TRUE;
2313: pcbddc->nedfield = -1;
2314: pcbddc->nedglobal = PETSC_TRUE;
2315: pcbddc->graphmaxcount = PETSC_INT_MAX;
2316: pcbddc->sub_schurs_layers = -1;
2317: pcbddc->ksp_D = kspD;
2318: pcbddc->ksp_R = kspR;
2319: pcbddc->coarse_ksp = kspC;
2320: PetscFunctionReturn(PETSC_SUCCESS);
2321: }
2323: static PetscErrorCode PCDestroy_BDDC(PC pc)
2324: {
2325: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
2327: PetscFunctionBegin;
2328: PetscCall(PCReset_BDDC(pc));
2329: PetscCall(KSPDestroy(&pcbddc->ksp_D));
2330: PetscCall(KSPDestroy(&pcbddc->ksp_R));
2331: PetscCall(KSPDestroy(&pcbddc->coarse_ksp));
2332: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDiscreteGradient_C", NULL));
2333: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDivergenceMat_C", NULL));
2334: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetChangeOfBasisMat_C", NULL));
2335: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesLocalIS_C", NULL));
2336: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesIS_C", NULL));
2337: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesLocalIS_C", NULL));
2338: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesIS_C", NULL));
2339: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetCoarseningRatio_C", NULL));
2340: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevel_C", NULL));
2341: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetUseExactDirichlet_C", NULL));
2342: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevels_C", NULL));
2343: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCLoadCustomization_C", NULL));
2344: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSaveCustomization_C", NULL));
2345: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundaries_C", NULL));
2346: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundariesLocal_C", NULL));
2347: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundaries_C", NULL));
2348: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundariesLocal_C", NULL));
2349: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundaries_C", NULL));
2350: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundariesLocal_C", NULL));
2351: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundaries_C", NULL));
2352: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundariesLocal_C", NULL));
2353: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplitting_C", NULL));
2354: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplittingLocal_C", NULL));
2355: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", NULL));
2356: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCCreateFETIDPOperators_C", NULL));
2357: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetRHS_C", NULL));
2358: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetSolution_C", NULL));
2359: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", NULL));
2360: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
2361: PetscCall(PetscFree(pc->data));
2362: PetscFunctionReturn(PETSC_SUCCESS);
2363: }
2365: static PetscErrorCode PCSetCoordinates_BDDC(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
2366: {
2367: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
2368: PCBDDCGraph mat_graph = pcbddc->mat_graph;
2370: PetscFunctionBegin;
2371: PetscCall(PetscFree(mat_graph->coords));
2372: PetscCall(PetscMalloc1(nloc * dim, &mat_graph->coords));
2373: PetscCall(PetscArraycpy(mat_graph->coords, coords, nloc * dim));
2374: mat_graph->cnloc = nloc;
2375: mat_graph->cdim = dim;
2376: mat_graph->cloc = PETSC_FALSE;
2377: /* flg setup */
2378: pcbddc->recompute_topography = PETSC_TRUE;
2379: pcbddc->corner_selected = PETSC_FALSE;
2380: PetscFunctionReturn(PETSC_SUCCESS);
2381: }
2383: static PetscErrorCode PCPreSolveChangeRHS_BDDC(PC pc, PetscBool *change)
2384: {
2385: PetscFunctionBegin;
2386: *change = PETSC_TRUE;
2387: PetscFunctionReturn(PETSC_SUCCESS);
2388: }
2390: static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
2391: {
2392: FETIDPMat_ctx mat_ctx;
2393: Vec work;
2394: PC_IS *pcis;
2395: PC_BDDC *pcbddc;
2397: PetscFunctionBegin;
2398: PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2399: pcis = (PC_IS *)mat_ctx->pc->data;
2400: pcbddc = (PC_BDDC *)mat_ctx->pc->data;
2402: PetscCall(VecSet(fetidp_flux_rhs, 0.0));
2403: /* copy rhs since we may change it during PCPreSolve_BDDC */
2404: if (!pcbddc->original_rhs) PetscCall(VecDuplicate(pcis->vec1_global, &pcbddc->original_rhs));
2405: if (mat_ctx->rhs_flip) {
2406: PetscCall(VecPointwiseMult(pcbddc->original_rhs, standard_rhs, mat_ctx->rhs_flip));
2407: } else {
2408: PetscCall(VecCopy(standard_rhs, pcbddc->original_rhs));
2409: }
2410: if (mat_ctx->g2g_p) {
2411: /* interface pressure rhs */
2412: PetscCall(VecScatterBegin(mat_ctx->g2g_p, fetidp_flux_rhs, pcbddc->original_rhs, INSERT_VALUES, SCATTER_REVERSE));
2413: PetscCall(VecScatterEnd(mat_ctx->g2g_p, fetidp_flux_rhs, pcbddc->original_rhs, INSERT_VALUES, SCATTER_REVERSE));
2414: PetscCall(VecScatterBegin(mat_ctx->g2g_p, standard_rhs, fetidp_flux_rhs, INSERT_VALUES, SCATTER_FORWARD));
2415: PetscCall(VecScatterEnd(mat_ctx->g2g_p, standard_rhs, fetidp_flux_rhs, INSERT_VALUES, SCATTER_FORWARD));
2416: if (!mat_ctx->rhs_flip) PetscCall(VecScale(fetidp_flux_rhs, -1.));
2417: }
2418: /*
2419: change of basis for physical rhs if needed
2420: It also changes the rhs in case of dirichlet boundaries
2421: */
2422: PetscCall(PCPreSolve_BDDC(mat_ctx->pc, NULL, pcbddc->original_rhs, NULL));
2423: if (pcbddc->ChangeOfBasisMatrix) {
2424: PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, pcbddc->original_rhs, pcbddc->work_change));
2425: work = pcbddc->work_change;
2426: } else {
2427: work = pcbddc->original_rhs;
2428: }
2429: /* store vectors for computation of fetidp final solution */
2430: PetscCall(VecScatterBegin(pcis->global_to_D, work, mat_ctx->temp_solution_D, INSERT_VALUES, SCATTER_FORWARD));
2431: PetscCall(VecScatterEnd(pcis->global_to_D, work, mat_ctx->temp_solution_D, INSERT_VALUES, SCATTER_FORWARD));
2432: /* scale rhs since it should be unassembled */
2433: /* TODO use counter scaling? (also below) */
2434: PetscCall(VecScatterBegin(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2435: PetscCall(VecScatterEnd(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2436: /* Apply partition of unity */
2437: PetscCall(VecPointwiseMult(mat_ctx->temp_solution_B, pcis->D, mat_ctx->temp_solution_B));
2438: /* PetscCall(PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B)); */
2439: if (!pcbddc->switch_static) {
2440: /* compute partially subassembled Schur complement right-hand side */
2441: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2442: PetscCall(KSPSolve(pcbddc->ksp_D, mat_ctx->temp_solution_D, pcis->vec1_D));
2443: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2444: /* Cannot propagate up error in KSPSolve() because there is no access to the PC */
2445: PetscCall(MatMult(pcis->A_BI, pcis->vec1_D, pcis->vec1_B));
2446: PetscCall(VecAXPY(mat_ctx->temp_solution_B, -1.0, pcis->vec1_B));
2447: PetscCall(VecSet(work, 0.0));
2448: PetscCall(VecScatterBegin(pcis->global_to_B, mat_ctx->temp_solution_B, work, ADD_VALUES, SCATTER_REVERSE));
2449: PetscCall(VecScatterEnd(pcis->global_to_B, mat_ctx->temp_solution_B, work, ADD_VALUES, SCATTER_REVERSE));
2450: /* PetscCall(PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B)); */
2451: PetscCall(VecScatterBegin(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2452: PetscCall(VecScatterEnd(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2453: PetscCall(VecPointwiseMult(mat_ctx->temp_solution_B, pcis->D, mat_ctx->temp_solution_B));
2454: }
2455: /* BDDC rhs */
2456: PetscCall(VecCopy(mat_ctx->temp_solution_B, pcis->vec1_B));
2457: if (pcbddc->switch_static) PetscCall(VecCopy(mat_ctx->temp_solution_D, pcis->vec1_D));
2458: /* apply BDDC */
2459: PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
2460: PetscCall(PCBDDCApplyInterfacePreconditioner(mat_ctx->pc, PETSC_FALSE));
2461: PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
2463: /* Application of B_delta and assembling of rhs for fetidp fluxes */
2464: PetscCall(MatMult(mat_ctx->B_delta, pcis->vec1_B, mat_ctx->lambda_local));
2465: PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, mat_ctx->lambda_local, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2466: PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, mat_ctx->lambda_local, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2467: /* Add contribution to interface pressures */
2468: if (mat_ctx->l2g_p) {
2469: PetscCall(VecISSet(pcis->vec1_B, mat_ctx->lP_B, 0));
2470: PetscCall(MatMult(mat_ctx->B_BB, pcis->vec1_B, mat_ctx->vP));
2471: if (pcbddc->switch_static) {
2472: PetscCall(VecISSet(pcis->vec1_D, mat_ctx->lP_I, 0));
2473: PetscCall(MatMultAdd(mat_ctx->B_BI, pcis->vec1_D, mat_ctx->vP, mat_ctx->vP));
2474: }
2475: PetscCall(VecScatterBegin(mat_ctx->l2g_p, mat_ctx->vP, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2476: PetscCall(VecScatterEnd(mat_ctx->l2g_p, mat_ctx->vP, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2477: }
2478: PetscFunctionReturn(PETSC_SUCCESS);
2479: }
2481: /*@
2482: PCBDDCMatFETIDPGetRHS - Computes the right-hand side of a FETI-DP system from the physical right-hand side
2484: Collective
2486: Input Parameters:
2487: + fetidp_mat - the FETI-DP matrix obtained from `PCBDDCCreateFETIDPOperators()`
2488: - standard_rhs - the right-hand side of the original linear system
2490: Output Parameter:
2491: . fetidp_flux_rhs - vector in which to store the right-hand side of the FETI-DP system
2493: Level: developer
2495: Note:
2496: Most users should employ the `KSP` interface for linear solvers and create a solver of type `KSPFETIDP`.
2498: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCCreateFETIDPOperators()`, `PCBDDCMatFETIDPGetSolution()`
2499: @*/
2500: PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
2501: {
2502: FETIDPMat_ctx mat_ctx;
2504: PetscFunctionBegin;
2508: PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2509: PetscUseMethod(mat_ctx->pc, "PCBDDCMatFETIDPGetRHS_C", (Mat, Vec, Vec), (fetidp_mat, standard_rhs, fetidp_flux_rhs));
2510: PetscFunctionReturn(PETSC_SUCCESS);
2511: }
2513: static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
2514: {
2515: FETIDPMat_ctx mat_ctx;
2516: PC_IS *pcis;
2517: PC_BDDC *pcbddc;
2518: Vec work;
2520: PetscFunctionBegin;
2521: PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2522: pcis = (PC_IS *)mat_ctx->pc->data;
2523: pcbddc = (PC_BDDC *)mat_ctx->pc->data;
2525: /* apply B_delta^T */
2526: PetscCall(VecSet(pcis->vec1_B, 0.));
2527: PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, fetidp_flux_sol, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
2528: PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, fetidp_flux_sol, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
2529: PetscCall(MatMultTranspose(mat_ctx->B_delta, mat_ctx->lambda_local, pcis->vec1_B));
2530: if (mat_ctx->l2g_p) {
2531: PetscCall(VecScatterBegin(mat_ctx->l2g_p, fetidp_flux_sol, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
2532: PetscCall(VecScatterEnd(mat_ctx->l2g_p, fetidp_flux_sol, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
2533: PetscCall(MatMultAdd(mat_ctx->Bt_BB, mat_ctx->vP, pcis->vec1_B, pcis->vec1_B));
2534: }
2536: /* compute rhs for BDDC application */
2537: PetscCall(VecAYPX(pcis->vec1_B, -1.0, mat_ctx->temp_solution_B));
2538: if (pcbddc->switch_static) {
2539: PetscCall(VecCopy(mat_ctx->temp_solution_D, pcis->vec1_D));
2540: if (mat_ctx->l2g_p) {
2541: PetscCall(VecScale(mat_ctx->vP, -1.));
2542: PetscCall(MatMultAdd(mat_ctx->Bt_BI, mat_ctx->vP, pcis->vec1_D, pcis->vec1_D));
2543: }
2544: }
2546: /* apply BDDC */
2547: PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
2548: PetscCall(PCBDDCApplyInterfacePreconditioner(mat_ctx->pc, PETSC_FALSE));
2550: /* put values into global vector */
2551: if (pcbddc->ChangeOfBasisMatrix) work = pcbddc->work_change;
2552: else work = standard_sol;
2553: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, work, INSERT_VALUES, SCATTER_REVERSE));
2554: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, work, INSERT_VALUES, SCATTER_REVERSE));
2555: if (!pcbddc->switch_static) {
2556: /* compute values into the interior if solved for the partially subassembled Schur complement */
2557: PetscCall(MatMult(pcis->A_IB, pcis->vec1_B, pcis->vec1_D));
2558: PetscCall(VecAYPX(pcis->vec1_D, -1.0, mat_ctx->temp_solution_D));
2559: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2560: PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec1_D));
2561: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2562: /* Cannot propagate up error in KSPSolve() because there is no access to the PC */
2563: }
2565: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec1_D, work, INSERT_VALUES, SCATTER_REVERSE));
2566: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec1_D, work, INSERT_VALUES, SCATTER_REVERSE));
2567: /* add p0 solution to final solution */
2568: PetscCall(PCBDDCBenignGetOrSetP0(mat_ctx->pc, work, PETSC_FALSE));
2569: if (pcbddc->ChangeOfBasisMatrix) PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, work, standard_sol));
2570: PetscCall(PCPostSolve_BDDC(mat_ctx->pc, NULL, NULL, standard_sol));
2571: if (mat_ctx->g2g_p) {
2572: PetscCall(VecScatterBegin(mat_ctx->g2g_p, fetidp_flux_sol, standard_sol, INSERT_VALUES, SCATTER_REVERSE));
2573: PetscCall(VecScatterEnd(mat_ctx->g2g_p, fetidp_flux_sol, standard_sol, INSERT_VALUES, SCATTER_REVERSE));
2574: }
2575: PetscFunctionReturn(PETSC_SUCCESS);
2576: }
2578: static PetscErrorCode PCView_BDDCIPC(PC pc, PetscViewer viewer)
2579: {
2580: BDDCIPC_ctx bddcipc_ctx;
2581: PetscBool isascii;
2583: PetscFunctionBegin;
2584: PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2585: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2586: if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "BDDC interface preconditioner\n"));
2587: PetscCall(PetscViewerASCIIPushTab(viewer));
2588: PetscCall(PCView(bddcipc_ctx->bddc, viewer));
2589: PetscCall(PetscViewerASCIIPopTab(viewer));
2590: PetscFunctionReturn(PETSC_SUCCESS);
2591: }
2593: static PetscErrorCode PCSetUp_BDDCIPC(PC pc)
2594: {
2595: BDDCIPC_ctx bddcipc_ctx;
2596: PetscBool isbddc;
2597: Vec vv;
2598: IS is;
2599: PC_IS *pcis;
2601: PetscFunctionBegin;
2602: PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2603: PetscCall(PetscObjectTypeCompare((PetscObject)bddcipc_ctx->bddc, PCBDDC, &isbddc));
2604: PetscCheck(isbddc, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Invalid type %s. Must be of type bddc", ((PetscObject)bddcipc_ctx->bddc)->type_name);
2605: PetscCall(PCSetUp(bddcipc_ctx->bddc));
2607: /* create interface scatter */
2608: pcis = (PC_IS *)bddcipc_ctx->bddc->data;
2609: PetscCall(VecScatterDestroy(&bddcipc_ctx->g2l));
2610: PetscCall(MatCreateVecs(pc->pmat, &vv, NULL));
2611: PetscCall(ISRenumber(pcis->is_B_global, NULL, NULL, &is));
2612: PetscCall(VecScatterCreate(vv, is, pcis->vec1_B, NULL, &bddcipc_ctx->g2l));
2613: PetscCall(ISDestroy(&is));
2614: PetscCall(VecDestroy(&vv));
2615: PetscFunctionReturn(PETSC_SUCCESS);
2616: }
2618: static PetscErrorCode PCApply_BDDCIPC(PC pc, Vec r, Vec x)
2619: {
2620: BDDCIPC_ctx bddcipc_ctx;
2621: PC_IS *pcis;
2622: VecScatter tmps;
2624: PetscFunctionBegin;
2625: PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2626: pcis = (PC_IS *)bddcipc_ctx->bddc->data;
2627: tmps = pcis->global_to_B;
2628: pcis->global_to_B = bddcipc_ctx->g2l;
2629: PetscCall(PCBDDCScalingRestriction(bddcipc_ctx->bddc, r, pcis->vec1_B));
2630: PetscCall(PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc, PETSC_FALSE));
2631: PetscCall(PCBDDCScalingExtension(bddcipc_ctx->bddc, pcis->vec1_B, x));
2632: pcis->global_to_B = tmps;
2633: PetscFunctionReturn(PETSC_SUCCESS);
2634: }
2636: static PetscErrorCode PCApplyTranspose_BDDCIPC(PC pc, Vec r, Vec x)
2637: {
2638: BDDCIPC_ctx bddcipc_ctx;
2639: PC_IS *pcis;
2640: VecScatter tmps;
2642: PetscFunctionBegin;
2643: PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2644: pcis = (PC_IS *)bddcipc_ctx->bddc->data;
2645: tmps = pcis->global_to_B;
2646: pcis->global_to_B = bddcipc_ctx->g2l;
2647: PetscCall(PCBDDCScalingRestriction(bddcipc_ctx->bddc, r, pcis->vec1_B));
2648: PetscCall(PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc, PETSC_TRUE));
2649: PetscCall(PCBDDCScalingExtension(bddcipc_ctx->bddc, pcis->vec1_B, x));
2650: pcis->global_to_B = tmps;
2651: PetscFunctionReturn(PETSC_SUCCESS);
2652: }
2654: static PetscErrorCode PCDestroy_BDDCIPC(PC pc)
2655: {
2656: BDDCIPC_ctx bddcipc_ctx;
2658: PetscFunctionBegin;
2659: PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2660: PetscCall(PCDestroy(&bddcipc_ctx->bddc));
2661: PetscCall(VecScatterDestroy(&bddcipc_ctx->g2l));
2662: PetscCall(PetscFree(bddcipc_ctx));
2663: PetscFunctionReturn(PETSC_SUCCESS);
2664: }
2666: /*@
2667: PCBDDCMatFETIDPGetSolution - Computes the physical solution from the solution of a FETI-DP system
2669: Collective
2671: Input Parameters:
2672: + fetidp_mat - the FETI-DP matrix obtained from `PCBDDCCreateFETIDPOperators()`
2673: - fetidp_flux_sol - the solution of the FETI-DP linear system
2675: Output Parameter:
2676: . standard_sol - vector in which to store the solution on the physical domain
2678: Level: developer
2680: Note:
2681: Most users should employ the `KSP` interface for linear solvers and create a solver of type `KSPFETIDP`.
2683: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCCreateFETIDPOperators()`, `PCBDDCMatFETIDPGetRHS()`
2684: @*/
2685: PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
2686: {
2687: FETIDPMat_ctx mat_ctx;
2689: PetscFunctionBegin;
2693: PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2694: PetscUseMethod(mat_ctx->pc, "PCBDDCMatFETIDPGetSolution_C", (Mat, Vec, Vec), (fetidp_mat, fetidp_flux_sol, standard_sol));
2695: PetscFunctionReturn(PETSC_SUCCESS);
2696: }
2698: static PetscErrorCode MatISSubMatrixEmbedLocalIS(Mat A, IS oldis, IS *newis)
2699: {
2700: Mat_IS *matis = (Mat_IS *)A->data;
2701: ISLocalToGlobalMapping ltog;
2702: IS is;
2704: PetscFunctionBegin;
2705: PetscCheck(matis->getsub_ris, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing getsub IS");
2706: PetscCall(ISLocalToGlobalMappingCreateIS(matis->getsub_ris, <og));
2707: PetscCall(ISGlobalToLocalMappingApplyIS(ltog, IS_GTOLM_DROP, oldis, &is));
2708: PetscCall(ISOnComm(is, PetscObjectComm((PetscObject)A), PETSC_COPY_VALUES, newis));
2709: PetscCall(ISLocalToGlobalMappingDestroy(<og));
2710: PetscCall(ISDestroy(&is));
2711: PetscFunctionReturn(PETSC_SUCCESS);
2712: }
2714: static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, PetscBool fully_redundant, const char *prefix, Mat *fetidp_mat, PC *fetidp_pc)
2715: {
2716: FETIDPMat_ctx fetidpmat_ctx;
2717: Mat newmat;
2718: FETIDPPC_ctx fetidppc_ctx;
2719: PC newpc;
2720: MPI_Comm comm;
2722: PetscFunctionBegin;
2723: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
2724: /* FETI-DP matrix */
2725: PetscCall(PCBDDCCreateFETIDPMatContext(pc, &fetidpmat_ctx));
2726: fetidpmat_ctx->fully_redundant = fully_redundant;
2727: PetscCall(PCBDDCSetupFETIDPMatContext(fetidpmat_ctx));
2728: PetscCall(MatCreateShell(comm, fetidpmat_ctx->n, fetidpmat_ctx->n, fetidpmat_ctx->N, fetidpmat_ctx->N, fetidpmat_ctx, &newmat));
2729: PetscCall(PetscObjectSetName((PetscObject)newmat, !fetidpmat_ctx->l2g_lambda_only ? "F" : "G"));
2730: PetscCall(MatShellSetOperation(newmat, MATOP_MULT, (PetscErrorCodeFn *)FETIDPMatMult));
2731: PetscCall(MatShellSetOperation(newmat, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)FETIDPMatMultTranspose));
2732: PetscCall(MatShellSetOperation(newmat, MATOP_DESTROY, (PetscErrorCodeFn *)PCBDDCDestroyFETIDPMat));
2733: /* propagate MatOptions */
2734: {
2735: PC_BDDC *pcbddc = (PC_BDDC *)fetidpmat_ctx->pc->data;
2736: PetscBool isset, issym;
2738: PetscCall(MatIsSymmetricKnown(pc->mat, &isset, &issym));
2739: if ((isset && issym) || pcbddc->symmetric_primal) PetscCall(MatSetOption(newmat, MAT_SYMMETRIC, PETSC_TRUE));
2740: }
2741: PetscCall(MatSetOptionsPrefix(newmat, prefix));
2742: PetscCall(MatAppendOptionsPrefix(newmat, "fetidp_"));
2743: PetscCall(MatSetUp(newmat));
2744: /* FETI-DP preconditioner */
2745: PetscCall(PCBDDCCreateFETIDPPCContext(pc, &fetidppc_ctx));
2746: PetscCall(PCBDDCSetupFETIDPPCContext(newmat, fetidppc_ctx));
2747: PetscCall(PCCreate(comm, &newpc));
2748: PetscCall(PCSetOperators(newpc, newmat, newmat));
2749: PetscCall(PCSetOptionsPrefix(newpc, prefix));
2750: PetscCall(PCAppendOptionsPrefix(newpc, "fetidp_"));
2751: PetscCall(PCSetErrorIfFailure(newpc, pc->erroriffailure));
2752: if (!fetidpmat_ctx->l2g_lambda_only) { /* standard FETI-DP */
2753: PetscCall(PCSetType(newpc, PCSHELL));
2754: PetscCall(PCShellSetName(newpc, "FETI-DP multipliers"));
2755: PetscCall(PCShellSetContext(newpc, fetidppc_ctx));
2756: PetscCall(PCShellSetApply(newpc, FETIDPPCApply));
2757: PetscCall(PCShellSetApplyTranspose(newpc, FETIDPPCApplyTranspose));
2758: PetscCall(PCShellSetView(newpc, FETIDPPCView));
2759: PetscCall(PCShellSetDestroy(newpc, PCBDDCDestroyFETIDPPC));
2760: } else { /* saddle-point FETI-DP */
2761: Mat M;
2762: PetscInt psize;
2763: PetscBool fake = PETSC_FALSE, isfieldsplit;
2765: PetscCall(ISViewFromOptions(fetidpmat_ctx->lagrange, NULL, "-lag_view"));
2766: PetscCall(ISViewFromOptions(fetidpmat_ctx->pressure, NULL, "-press_view"));
2767: PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_PPmat", (PetscObject *)&M));
2768: PetscCall(PCSetType(newpc, PCFIELDSPLIT));
2769: PetscCall(PCFieldSplitSetIS(newpc, "lag", fetidpmat_ctx->lagrange));
2770: PetscCall(PCFieldSplitSetIS(newpc, "p", fetidpmat_ctx->pressure));
2771: PetscCall(PCFieldSplitSetType(newpc, PC_COMPOSITE_SCHUR));
2772: PetscCall(PCFieldSplitSetSchurFactType(newpc, PC_FIELDSPLIT_SCHUR_FACT_DIAG));
2773: PetscCall(ISGetSize(fetidpmat_ctx->pressure, &psize));
2774: if (psize != M->rmap->N) {
2775: Mat M2;
2776: PetscInt lpsize;
2778: fake = PETSC_TRUE;
2779: PetscCall(ISGetLocalSize(fetidpmat_ctx->pressure, &lpsize));
2780: PetscCall(MatCreate(comm, &M2));
2781: PetscCall(MatSetType(M2, MATAIJ));
2782: PetscCall(MatSetSizes(M2, lpsize, lpsize, psize, psize));
2783: PetscCall(MatSetUp(M2));
2784: PetscCall(MatAssemblyBegin(M2, MAT_FINAL_ASSEMBLY));
2785: PetscCall(MatAssemblyEnd(M2, MAT_FINAL_ASSEMBLY));
2786: PetscCall(PCFieldSplitSetSchurPre(newpc, PC_FIELDSPLIT_SCHUR_PRE_USER, M2));
2787: PetscCall(MatDestroy(&M2));
2788: } else {
2789: PetscCall(PCFieldSplitSetSchurPre(newpc, PC_FIELDSPLIT_SCHUR_PRE_USER, M));
2790: }
2791: PetscCall(PCFieldSplitSetSchurScale(newpc, 1.0));
2793: /* we need to setfromoptions and setup here to access the blocks */
2794: PetscCall(PCSetFromOptions(newpc));
2795: PetscCall(PCSetUp(newpc));
2797: /* user may have changed the type (e.g. -fetidp_pc_type none) */
2798: PetscCall(PetscObjectTypeCompare((PetscObject)newpc, PCFIELDSPLIT, &isfieldsplit));
2799: if (isfieldsplit) {
2800: KSP *ksps;
2801: PC ppc, lagpc;
2802: PetscInt nn;
2803: PetscBool ismatis, matisok = PETSC_FALSE, check = PETSC_FALSE;
2805: /* set the solver for the (0,0) block */
2806: PetscCall(PCFieldSplitSchurGetSubKSP(newpc, &nn, &ksps));
2807: if (!nn) { /* not of type PC_COMPOSITE_SCHUR */
2808: PetscCall(PCFieldSplitGetSubKSP(newpc, &nn, &ksps));
2809: if (!fake) { /* pass pmat to the pressure solver */
2810: Mat F;
2812: PetscCall(KSPGetOperators(ksps[1], &F, NULL));
2813: PetscCall(KSPSetOperators(ksps[1], F, M));
2814: }
2815: } else {
2816: PetscBool issym, isset;
2817: Mat S;
2819: PetscCall(PCFieldSplitSchurGetS(newpc, &S));
2820: PetscCall(MatIsSymmetricKnown(newmat, &isset, &issym));
2821: if (isset) PetscCall(MatSetOption(S, MAT_SYMMETRIC, issym));
2822: }
2823: PetscCall(KSPGetPC(ksps[0], &lagpc));
2824: PetscCall(PCSetType(lagpc, PCSHELL));
2825: PetscCall(PCShellSetName(lagpc, "FETI-DP multipliers"));
2826: PetscCall(PCShellSetContext(lagpc, fetidppc_ctx));
2827: PetscCall(PCShellSetApply(lagpc, FETIDPPCApply));
2828: PetscCall(PCShellSetApplyTranspose(lagpc, FETIDPPCApplyTranspose));
2829: PetscCall(PCShellSetView(lagpc, FETIDPPCView));
2830: PetscCall(PCShellSetDestroy(lagpc, PCBDDCDestroyFETIDPPC));
2832: /* Olof's idea: interface Schur complement preconditioner for the mass matrix */
2833: PetscCall(KSPGetPC(ksps[1], &ppc));
2834: if (fake) {
2835: PC_BDDC *pcbddc = (PC_BDDC *)fetidpmat_ctx->pc->data;
2836: BDDCIPC_ctx bddcipc_ctx;
2837: PetscContainer c;
2839: matisok = PETSC_TRUE;
2841: /* create inner BDDC solver */
2842: PetscCall(PetscNew(&bddcipc_ctx));
2843: PetscCall(PCCreate(comm, &bddcipc_ctx->bddc));
2844: PetscCall(PCSetType(bddcipc_ctx->bddc, PCBDDC));
2845: PetscCall(PCSetOperators(bddcipc_ctx->bddc, M, M));
2846: PetscCall(PetscObjectTypeCompare((PetscObject)M, MATIS, &ismatis));
2847: PetscCheck(ismatis, comm, PETSC_ERR_PLIB, "Matrix type %s not of type MATIS", ((PetscObject)M)->type_name);
2848: /* the inner bddc for FETI-DP is already setup, we have local info available */
2849: if (pcbddc->user_primal_vertices_local || pcbddc->n_ISForDofsLocal > 2) {
2850: if (pcbddc->user_primal_vertices_local) {
2851: IS primals;
2853: PetscCall(MatISSubMatrixEmbedLocalIS(M, pcbddc->user_primal_vertices_local, &primals));
2854: PetscCall(PCBDDCSetPrimalVerticesLocalIS(bddcipc_ctx->bddc, primals));
2855: PetscCall(ISDestroy(&primals));
2856: }
2857: if (pcbddc->n_ISForDofsLocal > 2) { /* no need to propagate info if nfields < 3 */
2858: IS *split;
2859: PetscInt i, nf;
2861: PetscCall(PetscCalloc1(pcbddc->n_ISForDofsLocal, &split));
2862: for (i = 0, nf = 0; i < pcbddc->n_ISForDofsLocal; i++) {
2863: PetscInt ns;
2865: PetscCall(MatISSubMatrixEmbedLocalIS(M, pcbddc->ISForDofsLocal[i], &split[nf]));
2866: PetscCall(ISGetSize(split[nf], &ns));
2867: if (!ns) PetscCall(ISDestroy(&split[nf]));
2868: else nf++;
2869: }
2870: PetscCall(PCBDDCSetDofsSplittingLocal(bddcipc_ctx->bddc, nf, split));
2871: for (i = 0; i < nf; i++) PetscCall(ISDestroy(&split[i]));
2872: PetscCall(PetscFree(split));
2873: }
2874: }
2875: PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_pCSR", (PetscObject *)&c));
2876: PetscCall(PetscObjectTypeCompare((PetscObject)M, MATIS, &ismatis));
2877: if (c && ismatis) {
2878: Mat lM;
2879: PetscInt *csr, n;
2881: PetscCall(MatISGetLocalMat(M, &lM));
2882: PetscCall(MatGetSize(lM, &n, NULL));
2883: PetscCall(PetscContainerGetPointer(c, &csr));
2884: PetscCall(PCBDDCSetLocalAdjacencyGraph(bddcipc_ctx->bddc, n, csr, csr + (n + 1), PETSC_COPY_VALUES));
2885: PetscCall(MatISRestoreLocalMat(M, &lM));
2886: }
2887: PetscCall(PCSetOptionsPrefix(bddcipc_ctx->bddc, ((PetscObject)ksps[1])->prefix));
2888: PetscCall(PCSetErrorIfFailure(bddcipc_ctx->bddc, pc->erroriffailure));
2889: PetscCall(PCSetFromOptions(bddcipc_ctx->bddc));
2891: /* wrap the interface application */
2892: PetscCall(PCSetType(ppc, PCSHELL));
2893: PetscCall(PCShellSetName(ppc, "FETI-DP pressure"));
2894: PetscCall(PCShellSetContext(ppc, bddcipc_ctx));
2895: PetscCall(PCShellSetSetUp(ppc, PCSetUp_BDDCIPC));
2896: PetscCall(PCShellSetApply(ppc, PCApply_BDDCIPC));
2897: PetscCall(PCShellSetApplyTranspose(ppc, PCApplyTranspose_BDDCIPC));
2898: PetscCall(PCShellSetView(ppc, PCView_BDDCIPC));
2899: PetscCall(PCShellSetDestroy(ppc, PCDestroy_BDDCIPC));
2900: }
2902: /* determine if we need to assemble M to construct a preconditioner */
2903: if (!matisok) {
2904: PetscCall(PetscObjectTypeCompare((PetscObject)M, MATIS, &ismatis));
2905: PetscCall(PetscObjectTypeCompareAny((PetscObject)ppc, &matisok, PCBDDC, PCJACOBI, PCNONE, PCMG, ""));
2906: if (ismatis && !matisok) PetscCall(MatConvert(M, MATAIJ, MAT_INPLACE_MATRIX, &M));
2907: }
2909: /* run the subproblems to check convergence */
2910: PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)newmat)->prefix, "-check_saddlepoint", &check, NULL));
2911: if (check) {
2912: for (PetscInt i = 0; i < nn; i++) {
2913: KSP kspC;
2914: PC npc;
2915: Mat F, pF;
2916: Vec x, y;
2917: PetscBool isschur, prec = PETSC_TRUE;
2919: PetscCall(KSPCreate(PetscObjectComm((PetscObject)ksps[i]), &kspC));
2920: PetscCall(KSPSetNestLevel(kspC, pc->kspnestlevel));
2921: PetscCall(KSPSetOptionsPrefix(kspC, ((PetscObject)ksps[i])->prefix));
2922: PetscCall(KSPAppendOptionsPrefix(kspC, "check_"));
2923: PetscCall(KSPGetOperators(ksps[i], &F, &pF));
2924: PetscCall(PetscObjectTypeCompare((PetscObject)F, MATSCHURCOMPLEMENT, &isschur));
2925: if (isschur) {
2926: KSP kspS, kspS2;
2927: Mat A00, pA00, A10, A01, A11;
2928: char prefix[256];
2930: PetscCall(MatSchurComplementGetKSP(F, &kspS));
2931: PetscCall(MatSchurComplementGetSubMatrices(F, &A00, &pA00, &A01, &A10, &A11));
2932: PetscCall(MatCreateSchurComplement(A00, pA00, A01, A10, A11, &F));
2933: PetscCall(MatSchurComplementGetKSP(F, &kspS2));
2934: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sschur_", ((PetscObject)kspC)->prefix));
2935: PetscCall(KSPSetOptionsPrefix(kspS2, prefix));
2936: PetscCall(KSPGetPC(kspS2, &npc));
2937: PetscCall(PCSetType(npc, PCKSP));
2938: PetscCall(PCKSPSetKSP(npc, kspS));
2939: PetscCall(KSPSetFromOptions(kspS2));
2940: PetscCall(KSPGetPC(kspS2, &npc));
2941: PetscCall(PCSetUseAmat(npc, PETSC_TRUE));
2942: } else {
2943: PetscCall(PetscObjectReference((PetscObject)F));
2944: }
2945: PetscCall(KSPSetFromOptions(kspC));
2946: PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)kspC)->prefix, "-preconditioned", &prec, NULL));
2947: if (prec) {
2948: PetscCall(KSPGetPC(ksps[i], &npc));
2949: PetscCall(KSPSetPC(kspC, npc));
2950: }
2951: PetscCall(KSPSetOperators(kspC, F, pF));
2952: PetscCall(MatCreateVecs(F, &x, &y));
2953: PetscCall(VecSetRandom(x, NULL));
2954: PetscCall(MatMult(F, x, y));
2955: PetscCall(KSPSolve(kspC, y, x));
2956: PetscCall(KSPCheckSolve(kspC, npc, x));
2957: PetscCall(KSPDestroy(&kspC));
2958: PetscCall(MatDestroy(&F));
2959: PetscCall(VecDestroy(&x));
2960: PetscCall(VecDestroy(&y));
2961: }
2962: }
2963: PetscCall(PetscFree(ksps));
2964: }
2965: }
2966: /* return pointers for objects created */
2967: *fetidp_mat = newmat;
2968: *fetidp_pc = newpc;
2969: PetscFunctionReturn(PETSC_SUCCESS);
2970: }
2972: /*@
2973: PCBDDCCreateFETIDPOperators - Creates the FETI-DP matrix and its Dirichlet preconditioner
2975: Collective
2977: Input Parameters:
2978: + pc - the `PCBDDC` preconditioning context, after `PCSetUp()` has been called
2979: . fully_redundant - `PETSC_TRUE` for a fully redundant set of Lagrange multipliers
2980: - prefix - options database prefix for the objects to be created, or `NULL`
2982: Output Parameters:
2983: + fetidp_mat - the FETI-DP shell matrix
2984: - fetidp_pc - the shell Dirichlet preconditioner for the FETI-DP matrix
2986: Level: developer
2988: Notes:
2989: Most users should employ the `KSP` interface for linear solvers and create a solver of type `KSPFETIDP`.
2990: The FETI-DP matrix supports `MatMult()` and `MatMultTranspose()`.
2992: The caller must destroy the returned objects with `MatDestroy()` and `PCDestroy()`.
2994: .seealso: [](ch_ksp), `KSPFETIDP`, `PCBDDC`, `PCBDDCMatFETIDPGetRHS()`, `PCBDDCMatFETIDPGetSolution()`
2995: @*/
2996: PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, PetscBool fully_redundant, const char *prefix, Mat *fetidp_mat, PC *fetidp_pc)
2997: {
2998: PetscFunctionBegin;
3000: PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "You must call PCSetup_BDDC() first");
3001: PetscUseMethod(pc, "PCBDDCCreateFETIDPOperators_C", (PC, PetscBool, const char *, Mat *, PC *), (pc, fully_redundant, prefix, fetidp_mat, fetidp_pc));
3002: PetscFunctionReturn(PETSC_SUCCESS);
3003: }
3005: /*MC
3006: PCBDDC - Balancing Domain Decomposition by Constraints preconditioner
3008: Options Database Keys:
3009: + -pc_bddc_use_vertices (true|false) - include vertices in the primal space
3010: . -pc_bddc_use_edges (true|false) - include edge constraints in the primal space
3011: . -pc_bddc_use_faces (true|false) - include face constraints in the primal space
3012: . -pc_bddc_vertex_size size - classify connected components of at most this size as primal vertices
3013: . -pc_bddc_corner_selection (true|false) - select corners using subdomain faces and coordinates
3014: . -pc_bddc_use_local_mat_graph (true|false) - use the local matrix adjacency graph for interface analysis
3015: . -pc_bddc_local_mat_graph_square count - number of times to square the local matrix graph before interface analysis
3016: . -pc_bddc_graph_maxcount count - classify components shared by more than this many neighboring subdomains as primal vertices
3017: . -pc_bddc_detect_disconnected (true|false) - detect disconnected local subdomains
3018: . -pc_bddc_detect_disconnected_filter (true|false) - filter small local matrix entries when detecting disconnected subdomains
3019: . -pc_bddc_monolithic (true|false) - discard information about splitting degrees of freedom into fields
3020: . -pc_bddc_use_nnsp (true|false) - use the matrix near nullspace to construct constraints
3021: . -pc_bddc_use_nnsp_true (true|false) - use the supplied near-nullspace vectors directly, without orthonormalization
3022: . -pc_bddc_constraint_near_null_space_tol tol - discard near-nullspace vectors whose restriction to a connected component has norm at most tol
3023: . -pc_bddc_constraint_singular_tol tol - relative tolerance for retaining independent constraint modes
3024: . -pc_bddc_symmetric (true|false) - compute primal basis functions assuming symmetry; use false for nonsymmetric problems
3025: . -pc_bddc_use_change_of_basis (true|false) - construct a change of basis on edges
3026: . -pc_bddc_use_change_on_faces (true|false) - construct the requested change of basis on faces
3027: . -pc_bddc_interface_ext_type (dirichlet|lump) - select how interface corrections are extended to subdomain interiors
3028: . -pc_bddc_dirichlet_approximate (true|false) - enable nullspace corrections for approximate Dirichlet solvers
3029: . -pc_bddc_dirichlet_approximate_scale (true|false) - scale the approximate Dirichlet solver when applying nullspace corrections
3030: . -pc_bddc_neumann_approximate (true|false) - enable nullspace corrections for approximate Neumann solvers
3031: . -pc_bddc_neumann_approximate_scale (true|false) - scale the approximate Neumann solver when applying nullspace corrections
3032: . -pc_bddc_switch_static (true|false) - switch from the default $M_2$ operator to $M_3$ in {cite}`dohrmann2007approximate`
3033: . -pc_bddc_levels levels - maximum number of additional levels (default 0)
3034: . -pc_bddc_coarsening_ratio ratio - target number of process subdomains or local elements per aggregate (default 8)
3035: . -pc_bddc_coarse_eqs_per_proc neq - target number of equations per process at the coarsest level; a negative value uses one process
3036: . -pc_bddc_coarse_eqs_limit neq - stop adding coarse levels when the coarse problem has at most this many equations
3037: . -pc_bddc_coarse_adj nprocs - number of processes used to hold the coarse adjacency graph for partitioning
3038: . -pc_bddc_use_coarse_estimates (true|false) - use estimated eigenvalues to configure an iterative coarse solver
3039: . -pc_bddc_use_deluxe_scaling (true|false) - use deluxe scaling
3040: . -pc_bddc_deluxe_zerorows (true|false) - zero rows and columns of deluxe operators associated with primal degrees of freedom
3041: . -pc_bddc_deluxe_singlemat (true|false) - combine the deluxe scaling operations into one matrix per interface component
3042: . -pc_bddc_schur_rebuild (true|false) - rebuild the interface graph without adjacency information for computing Schur complement principal minors
3043: . -pc_bddc_schur_layers layers - number of layers used for economic deluxe scaling; -1 uses all degrees of freedom
3044: . -pc_bddc_schur_use_useradj (true|false) - use the graph supplied with `PCBDDCSetLocalAdjacencyGraph()` to select Schur complement layers
3045: . -pc_bddc_schur_exact (true|false) - use the full Schur complement, including components of size one, for adaptive constraint selection
3046: . -pc_bddc_adaptive_threshold thresholds - one or two comma-separated thresholds for adaptive constraint selection; one value sets both thresholds
3047: . -pc_bddc_adaptive_nmin count - minimum number of adaptive constraints per connected component
3048: . -pc_bddc_adaptive_nmax count - maximum number of adaptive constraints per connected component
3049: . -pc_bddc_adaptive_userdefined (true|false) - retain constraints from `MatSetNearNullSpace()` in addition to adaptive constraints
3050: . -pc_bddc_benign_trick (true|false) - use the benign subspace approach for saddle-point problems with discontinuous spaces
3051: . -pc_bddc_nonetflux (true|false) - compute quadrature weights for no-net-flux constraints automatically
3052: . -pc_bddc_nedelec_field_primal (true|false) - make Nedelec degrees of freedom shared by more than two subdomains primal
3053: . -pc_bddc_nedelec_order order - override the Nedelec order for testing; 0 selects variable order
3054: . -pc_bddc_nedelec_print (true|false) - print Nedelec setup diagnostics
3055: . -pc_bddc_load filename - load BDDC customization from a binary file for debugging
3056: . -pc_bddc_load_version version - version of the customization file to load
3057: . -pc_bddc_save filename - save BDDC customization to a binary file after setup for debugging
3058: . -pc_bddc_save_version version - version of the customization file to write
3059: . -pc_bddc_check_level level - verbosity level of debugging output
3061: Level: intermediate
3063: Notes:
3064: `PCBDDC` requires `MATIS` matrices and supports nonsymmetric and indefinite problems.
3065: The implementation and its customization are described in {cite}`zampini2016pcbddc`; see also [](sec_bddc).
3066: See {cite}`dohrmann2007approximate`, {cite}`klawonn2006dual`, and {cite}`mandel2008multispace` for the underlying methods.
3068: `PCBDDC` acts on all degrees of freedom, including subdomain interiors. This allows approximate subdomain solvers.
3069: Approximate local solvers are automatically adapted as described in {cite}`dohrmann2007approximate` when a nullspace object
3070: is attached to the subdomain matrices and approximate solvers are selected through the options database.
3072: Interface nodes are classified as vertices, edges, or faces using the local-to-global mapping of degrees of freedom
3073: and the local connectivity graph. The graph can be customized with `PCBDDCSetLocalAdjacencyGraph()`.
3074: Additional information about degrees of freedom can be supplied with `PCBDDCSetDofsSplitting()`, `PCBDDCSetDirichletBoundaries()`,
3075: `PCBDDCSetNeumannBoundaries()`, `PCBDDCSetPrimalVerticesIS()`, and their local counterparts.
3077: Support for $H(\mathrm{div})$ and $H(\mathrm{curl})$ problems is provided through `PCBDDCSetDivergenceMat()` and `PCBDDCSetDiscreteGradient()`.
3079: Constraints can be customized by attaching a `MatNullSpace` object to the `MATIS` matrix with `MatSetNearNullSpace()`.
3080: Linearly independent modes are retained using a singular value decomposition.
3082: When requested, a change of basis is performed as in {cite}`klawonn2006dual`. Local QR factorizations are used when more than
3083: one constraint is present on a connected component, such as an edge or a face. A user-defined change of basis can be supplied
3084: with `PCBDDCSetChangeOfBasisMat()`.
3086: Multilevel `PCBDDC` is supported as described in {cite}`mandel2008multispace`. Process subdomains are partitioned using a `MatPartitioning` object.
3087: When a local `MATIS` matrix stores multiple elements, their local coarse contributions are first aggregated using a `PetscPartitioner` object.
3088: In this case, the coarsening ratio is the target number of local elements per aggregate.
3090: Adaptive selection of primal constraints is supported for symmetric positive definite systems with high contrast in the coefficients
3091: when MUMPS or MKL_PARDISO is available. See {cite}`ohwidlundzampinidohrmann2017` for adaptive deluxe methods for Raviart-Thomas fields.
3092: The benign subspace approach for saddle-point problems with discontinuous spaces is described in {cite}`zampinitu2017`.
3094: Options for the Dirichlet, Neumann, coarse solver, and aggregation objects use the following prefixes, preceded by any user prefix:
3095: .vb
3096: -pc_bddc_dirichlet_
3097: -pc_bddc_neumann_
3098: -pc_bddc_coarse_
3099: -pc_bddc_aggregator_n_
3100: .ve
3101: For example, `-pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg` selects an approximate Dirichlet solver.
3102: By default, local solvers use `KSPPREONLY` with a direct factorization.
3103: At level `n`, the aggregator prefix is `pc_bddc_aggregator_n_`. Numeric-prefix fallback allows
3104: `-pc_bddc_aggregator_mat_partitioning_type type` and `-pc_bddc_aggregator_petscpartitioner_type type` to configure all levels.
3106: For BDDC level `N` > 0, the solver prefixes are:
3107: .vb
3108: -pc_bddc_dirichlet_lN_
3109: -pc_bddc_neumann_lN_
3110: -pc_bddc_coarse_lN_
3111: .ve
3112: Level 0 is the finest level. A coarse-level `PCBDDC` inherits the corresponding coarse-solver prefix. For example,
3113: .vb
3114: -pc_bddc_coarse_pc_bddc_adaptive_threshold 5
3115: .ve
3116: sets the adaptive constraint threshold to 5 for the first coarse-level `PCBDDC`.
3118: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `MATIS`, `KSPFETIDP`, `PCLU`, `PCGAMG`, `PCBDDCSetLocalAdjacencyGraph()`, `PCBDDCSetDofsSplitting()`,
3119: `PCBDDCSetDirichletBoundaries()`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCSetPrimalVerticesIS()`, `MatNullSpace`, `MatSetNearNullSpace()`,
3120: `PCBDDCSetChangeOfBasisMat()`, `PCBDDCSetDivergenceMat()`, `PCBDDCSetDiscreteGradient()`
3121: M*/
3123: PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc)
3124: {
3125: PC_BDDC *pcbddc;
3127: PetscFunctionBegin;
3128: PetscCall(PetscNew(&pcbddc));
3129: pc->data = pcbddc;
3131: PetscCall(PCISInitialize(pc));
3133: /* create local graph structure */
3134: PetscCall(PCBDDCGraphCreate(&pcbddc->mat_graph));
3136: /* BDDC nonzero defaults */
3137: pcbddc->use_nnsp = PETSC_TRUE;
3138: pcbddc->use_local_adj = PETSC_TRUE;
3139: pcbddc->use_vertices = PETSC_TRUE;
3140: pcbddc->use_edges = PETSC_TRUE;
3141: pcbddc->symmetric_primal = PETSC_TRUE;
3142: pcbddc->vertex_size = 1;
3143: pcbddc->recompute_topography = PETSC_TRUE;
3144: pcbddc->coarse_size = -1;
3145: pcbddc->use_exact_dirichlet_trick = PETSC_TRUE;
3146: pcbddc->coarsening_ratio = 8;
3147: pcbddc->coarse_eqs_per_proc = 1;
3148: pcbddc->benign_compute_correction = PETSC_TRUE;
3149: pcbddc->nedfield = -1;
3150: pcbddc->nedglobal = PETSC_TRUE;
3151: pcbddc->graphmaxcount = PETSC_INT_MAX;
3152: pcbddc->sub_schurs_layers = -1;
3153: pcbddc->adaptive_threshold[0] = 0.0;
3154: pcbddc->adaptive_threshold[1] = 0.0;
3156: /* function pointers */
3157: pc->ops->apply = PCApply_BDDC;
3158: pc->ops->applytranspose = PCApplyTranspose_BDDC;
3159: pc->ops->setup = PCSetUp_BDDC;
3160: pc->ops->destroy = PCDestroy_BDDC;
3161: pc->ops->setfromoptions = PCSetFromOptions_BDDC;
3162: pc->ops->view = PCView_BDDC;
3163: pc->ops->applyrichardson = NULL;
3164: pc->ops->applysymmetricleft = NULL;
3165: pc->ops->applysymmetricright = NULL;
3166: pc->ops->presolve = PCPreSolve_BDDC;
3167: pc->ops->postsolve = PCPostSolve_BDDC;
3168: pc->ops->reset = PCReset_BDDC;
3170: /* composing function */
3171: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDiscreteGradient_C", PCBDDCSetDiscreteGradient_BDDC));
3172: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDivergenceMat_C", PCBDDCSetDivergenceMat_BDDC));
3173: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetChangeOfBasisMat_C", PCBDDCSetChangeOfBasisMat_BDDC));
3174: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesLocalIS_C", PCBDDCSetPrimalVerticesLocalIS_BDDC));
3175: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesIS_C", PCBDDCSetPrimalVerticesIS_BDDC));
3176: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesLocalIS_C", PCBDDCGetPrimalVerticesLocalIS_BDDC));
3177: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesIS_C", PCBDDCGetPrimalVerticesIS_BDDC));
3178: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetCoarseningRatio_C", PCBDDCSetCoarseningRatio_BDDC));
3179: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevel_C", PCBDDCSetLevel_BDDC));
3180: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetUseExactDirichlet_C", PCBDDCSetUseExactDirichlet_BDDC));
3181: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevels_C", PCBDDCSetLevels_BDDC));
3182: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCLoadCustomization_C", PCBDDCLoadCustomization_BDDC));
3183: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSaveCustomization_C", PCBDDCSaveCustomization_BDDC));
3184: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundaries_C", PCBDDCSetDirichletBoundaries_BDDC));
3185: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundariesLocal_C", PCBDDCSetDirichletBoundariesLocal_BDDC));
3186: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundaries_C", PCBDDCSetNeumannBoundaries_BDDC));
3187: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundariesLocal_C", PCBDDCSetNeumannBoundariesLocal_BDDC));
3188: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundaries_C", PCBDDCGetDirichletBoundaries_BDDC));
3189: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundariesLocal_C", PCBDDCGetDirichletBoundariesLocal_BDDC));
3190: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundaries_C", PCBDDCGetNeumannBoundaries_BDDC));
3191: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundariesLocal_C", PCBDDCGetNeumannBoundariesLocal_BDDC));
3192: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplitting_C", PCBDDCSetDofsSplitting_BDDC));
3193: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplittingLocal_C", PCBDDCSetDofsSplittingLocal_BDDC));
3194: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", PCBDDCSetLocalAdjacencyGraph_BDDC));
3195: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCCreateFETIDPOperators_C", PCBDDCCreateFETIDPOperators_BDDC));
3196: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetRHS_C", PCBDDCMatFETIDPGetRHS_BDDC));
3197: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetSolution_C", PCBDDCMatFETIDPGetSolution_BDDC));
3198: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", PCPreSolveChangeRHS_BDDC));
3199: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_BDDC));
3200: PetscFunctionReturn(PETSC_SUCCESS);
3201: }
3203: /*@
3204: PCBDDCInitializePackage - Initializes the `PCBDDC` package
3206: Not Collective
3208: Level: developer
3210: Note:
3211: This routine is called by `PCInitializePackage()`.
3213: .seealso: [](ch_ksp), `PetscInitialize()`, `PCBDDCFinalizePackage()`
3214: @*/
3215: PetscErrorCode PCBDDCInitializePackage(void)
3216: {
3217: int i;
3219: PetscFunctionBegin;
3220: if (PCBDDCPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
3221: PCBDDCPackageInitialized = PETSC_TRUE;
3222: PetscCall(PetscRegisterFinalize(PCBDDCFinalizePackage));
3224: /* general events */
3225: PetscCall(PetscLogEventRegister("PCBDDCTopo", PC_CLASSID, &PC_BDDC_Topology[0]));
3226: PetscCall(PetscLogEventRegister("PCBDDCLKSP", PC_CLASSID, &PC_BDDC_LocalSolvers[0]));
3227: PetscCall(PetscLogEventRegister("PCBDDCLWor", PC_CLASSID, &PC_BDDC_LocalWork[0]));
3228: PetscCall(PetscLogEventRegister("PCBDDCCorr", PC_CLASSID, &PC_BDDC_CorrectionSetUp[0]));
3229: PetscCall(PetscLogEventRegister("PCBDDCASet", PC_CLASSID, &PC_BDDC_ApproxSetUp[0]));
3230: PetscCall(PetscLogEventRegister("PCBDDCAApp", PC_CLASSID, &PC_BDDC_ApproxApply[0]));
3231: PetscCall(PetscLogEventRegister("PCBDDCCSet", PC_CLASSID, &PC_BDDC_CoarseSetUp[0]));
3232: PetscCall(PetscLogEventRegister("PCBDDCCKSP", PC_CLASSID, &PC_BDDC_CoarseSolver[0]));
3233: PetscCall(PetscLogEventRegister("PCBDDCAdap", PC_CLASSID, &PC_BDDC_AdaptiveSetUp[0]));
3234: PetscCall(PetscLogEventRegister("PCBDDCScal", PC_CLASSID, &PC_BDDC_Scaling[0]));
3235: PetscCall(PetscLogEventRegister("PCBDDCSchr", PC_CLASSID, &PC_BDDC_Schurs[0]));
3236: PetscCall(PetscLogEventRegister("PCBDDCDirS", PC_CLASSID, &PC_BDDC_Solves[0][0]));
3237: PetscCall(PetscLogEventRegister("PCBDDCNeuS", PC_CLASSID, &PC_BDDC_Solves[0][1]));
3238: PetscCall(PetscLogEventRegister("PCBDDCCoaS", PC_CLASSID, &PC_BDDC_Solves[0][2]));
3239: for (i = 1; i < PETSC_PCBDDC_MAXLEVELS; i++) {
3240: char ename[32];
3242: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCTopo l%02d", i));
3243: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Topology[i]));
3244: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCLKSP l%02d", i));
3245: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_LocalSolvers[i]));
3246: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCLWor l%02d", i));
3247: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_LocalWork[i]));
3248: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCorr l%02d", i));
3249: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_CorrectionSetUp[i]));
3250: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCASet l%02d", i));
3251: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_ApproxSetUp[i]));
3252: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCAApp l%02d", i));
3253: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_ApproxApply[i]));
3254: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCSet l%02d", i));
3255: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_CoarseSetUp[i]));
3256: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCKSP l%02d", i));
3257: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_CoarseSolver[i]));
3258: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCAdap l%02d", i));
3259: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_AdaptiveSetUp[i]));
3260: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCScal l%02d", i));
3261: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Scaling[i]));
3262: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCSchr l%02d", i));
3263: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Schurs[i]));
3264: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCDirS l%02d", i));
3265: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Solves[i][0]));
3266: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCNeuS l%02d", i));
3267: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Solves[i][1]));
3268: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCoaS l%02d", i));
3269: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Solves[i][2]));
3270: }
3271: PetscFunctionReturn(PETSC_SUCCESS);
3272: }
3274: /*@
3275: PCBDDCFinalizePackage - Finalizes the `PCBDDC` package
3277: Not Collective
3279: Level: developer
3281: Note:
3282: This routine is called automatically by `PetscFinalize()`.
3284: .seealso: [](ch_ksp), `PetscFinalize()`, `PCBDDCInitializePackage()`
3285: @*/
3286: PetscErrorCode PCBDDCFinalizePackage(void)
3287: {
3288: PetscFunctionBegin;
3289: PCBDDCPackageInitialized = PETSC_FALSE;
3290: PetscFunctionReturn(PETSC_SUCCESS);
3291: }