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: The computed flux constraints supplement the near-nullspace supplied with `MatSetNearNullSpace()`.
367: The near-nullspace attached to the preconditioning matrix is preserved.
369: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDiscreteGradient()`, `MatSetNearNullSpace()`
370: @*/
371: PetscErrorCode PCBDDCSetDivergenceMat(PC pc, Mat divudotp, PetscBool trans, IS vl2l)
372: {
373: PetscBool ismatis;
375: PetscFunctionBegin;
378: PetscCheckSameComm(pc, 1, divudotp, 2);
381: PetscCall(PetscObjectTypeCompare((PetscObject)divudotp, MATIS, &ismatis));
382: PetscCheck(ismatis, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Divergence matrix needs to be of type MATIS");
383: PetscTryMethod(pc, "PCBDDCSetDivergenceMat_C", (PC, Mat, PetscBool, IS), (pc, divudotp, trans, vl2l));
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: static PetscErrorCode PCBDDCSetChangeOfBasisMat_BDDC(PC pc, Mat change, PetscBool interior)
388: {
389: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
391: PetscFunctionBegin;
392: PetscCall(PetscObjectReference((PetscObject)change));
393: PetscCall(MatDestroy(&pcbddc->user_ChangeOfBasisMatrix));
394: pcbddc->user_ChangeOfBasisMatrix = change;
395: pcbddc->change_interior = interior;
396: PetscFunctionReturn(PETSC_SUCCESS);
397: }
399: /*@
400: PCBDDCSetChangeOfBasisMat - Sets a user-defined change of basis for `PCBDDC`
402: Collective
404: Input Parameters:
405: + pc - the preconditioning context
406: . change - the change-of-basis matrix, with the same global and local sizes as the operator
407: - interior - `PETSC_TRUE` if the change of basis modifies interior degrees of freedom
409: Level: intermediate
411: .seealso: [](ch_ksp), `PCBDDC`
412: @*/
413: PetscErrorCode PCBDDCSetChangeOfBasisMat(PC pc, Mat change, PetscBool interior)
414: {
415: PetscFunctionBegin;
418: PetscCheckSameComm(pc, 1, change, 2);
419: if (pc->mat) {
420: PetscInt rows_c, cols_c, rows, cols;
421: PetscCall(MatGetSize(pc->mat, &rows, &cols));
422: PetscCall(MatGetSize(change, &rows_c, &cols_c));
423: 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);
424: 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);
425: PetscCall(MatGetLocalSize(pc->mat, &rows, &cols));
426: PetscCall(MatGetLocalSize(change, &rows_c, &cols_c));
427: 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);
428: 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);
429: }
430: PetscTryMethod(pc, "PCBDDCSetChangeOfBasisMat_C", (PC, Mat, PetscBool), (pc, change, interior));
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: static PetscErrorCode PCBDDCSetPrimalVerticesIS_BDDC(PC pc, IS PrimalVertices)
435: {
436: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
437: PetscBool isequal = PETSC_FALSE;
439: PetscFunctionBegin;
440: PetscCall(PetscObjectReference((PetscObject)PrimalVertices));
441: if (pcbddc->user_primal_vertices) PetscCall(ISEqual(PrimalVertices, pcbddc->user_primal_vertices, &isequal));
442: PetscCall(ISDestroy(&pcbddc->user_primal_vertices));
443: PetscCall(ISDestroy(&pcbddc->user_primal_vertices_local));
444: pcbddc->user_primal_vertices = PrimalVertices;
445: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: /*@
450: PCBDDCSetPrimalVerticesIS - Sets additional user-defined primal vertices in global numbering for `PCBDDC`
452: Collective
454: Input Parameters:
455: + pc - the preconditioning context
456: - PrimalVertices - index set of primal vertices in global numbering (can be empty)
458: Level: intermediate
460: Note:
461: Any process can list any global degree of freedom.
463: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
464: @*/
465: PetscErrorCode PCBDDCSetPrimalVerticesIS(PC pc, IS PrimalVertices)
466: {
467: PetscFunctionBegin;
470: PetscCheckSameComm(pc, 1, PrimalVertices, 2);
471: PetscTryMethod(pc, "PCBDDCSetPrimalVerticesIS_C", (PC, IS), (pc, PrimalVertices));
472: PetscFunctionReturn(PETSC_SUCCESS);
473: }
475: static PetscErrorCode PCBDDCGetPrimalVerticesIS_BDDC(PC pc, IS *is)
476: {
477: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
479: PetscFunctionBegin;
480: *is = pcbddc->user_primal_vertices;
481: PetscFunctionReturn(PETSC_SUCCESS);
482: }
484: /*@
485: PCBDDCGetPrimalVerticesIS - Gets the user-defined primal vertices in global numbering
487: Not Collective
489: Input Parameter:
490: . pc - the preconditioning context
492: Output Parameter:
493: . is - index set of primal vertices in global numbering (`NULL` if not set)
495: Level: intermediate
497: Note:
498: The returned `IS` is owned by `pc`; the caller must not destroy it.
500: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
501: @*/
502: PetscErrorCode PCBDDCGetPrimalVerticesIS(PC pc, IS *is)
503: {
504: PetscFunctionBegin;
506: PetscAssertPointer(is, 2);
507: PetscUseMethod(pc, "PCBDDCGetPrimalVerticesIS_C", (PC, IS *), (pc, is));
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }
511: static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices)
512: {
513: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
514: PetscBool isequal = PETSC_FALSE;
516: PetscFunctionBegin;
517: PetscCall(PetscObjectReference((PetscObject)PrimalVertices));
518: if (pcbddc->user_primal_vertices_local) PetscCall(ISEqual(PrimalVertices, pcbddc->user_primal_vertices_local, &isequal));
519: PetscCall(ISDestroy(&pcbddc->user_primal_vertices));
520: PetscCall(ISDestroy(&pcbddc->user_primal_vertices_local));
521: pcbddc->user_primal_vertices_local = PrimalVertices;
522: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: /*@
527: PCBDDCSetPrimalVerticesLocalIS - Sets additional user-defined primal vertices in local numbering for `PCBDDC`
529: Collective
531: Input Parameters:
532: + pc - the preconditioning context
533: - PrimalVertices - index set of primal vertices in the numbering of the local `MATIS` matrix (can be empty)
535: Level: intermediate
537: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
538: @*/
539: PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices)
540: {
541: PetscFunctionBegin;
544: PetscCheckSameComm(pc, 1, PrimalVertices, 2);
545: PetscTryMethod(pc, "PCBDDCSetPrimalVerticesLocalIS_C", (PC, IS), (pc, PrimalVertices));
546: PetscFunctionReturn(PETSC_SUCCESS);
547: }
549: static PetscErrorCode PCBDDCGetPrimalVerticesLocalIS_BDDC(PC pc, IS *is)
550: {
551: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
553: PetscFunctionBegin;
554: *is = pcbddc->user_primal_vertices_local;
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }
558: /*@
559: PCBDDCGetPrimalVerticesLocalIS - Gets the user-defined primal vertices in local numbering
561: Not Collective
563: Input Parameter:
564: . pc - the preconditioning context
566: Output Parameter:
567: . is - index set of primal vertices in the numbering of the local `MATIS` matrix, or `NULL` if unavailable
569: Level: intermediate
571: Notes:
572: The index set is supplied by `PCBDDCSetPrimalVerticesLocalIS()` or obtained from `PCBDDCSetPrimalVerticesIS()` during `PCSetUp()`.
574: The returned `IS` is owned by `pc`; the caller must not destroy it.
576: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`
577: @*/
578: PetscErrorCode PCBDDCGetPrimalVerticesLocalIS(PC pc, IS *is)
579: {
580: PetscFunctionBegin;
582: PetscAssertPointer(is, 2);
583: PetscUseMethod(pc, "PCBDDCGetPrimalVerticesLocalIS_C", (PC, IS *), (pc, is));
584: PetscFunctionReturn(PETSC_SUCCESS);
585: }
587: static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc, PetscInt k)
588: {
589: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
591: PetscFunctionBegin;
592: pcbddc->coarsening_ratio = k;
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
596: /*@
597: PCBDDCSetCoarseningRatio - Sets the coarsening ratio used by multilevel `PCBDDC`
599: Logically Collective
601: Input Parameters:
602: + pc - the preconditioning context
603: - k - target number of process subdomains or local elements per aggregate
605: Options Database Key:
606: . -pc_bddc_coarsening_ratio k - set the coarsening ratio used in multilevel coarsening
608: Level: intermediate
610: Note:
611: Approximately `k` subdomains at the finer level are aggregated into a single subdomain at the coarser level.
612: When a local `MATIS` matrix stores multiple elements, `k` is the target number of local elements per aggregate.
614: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetLevels()`
615: @*/
616: PetscErrorCode PCBDDCSetCoarseningRatio(PC pc, PetscInt k)
617: {
618: PetscFunctionBegin;
621: PetscTryMethod(pc, "PCBDDCSetCoarseningRatio_C", (PC, PetscInt), (pc, k));
622: PetscFunctionReturn(PETSC_SUCCESS);
623: }
625: /* The following functions (PCBDDCSetUseExactDirichlet PCBDDCSetLevel) are not public */
626: static PetscErrorCode PCBDDCSetUseExactDirichlet_BDDC(PC pc, PetscBool flg)
627: {
628: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
630: PetscFunctionBegin;
631: pcbddc->use_exact_dirichlet_trick = flg;
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }
635: PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc, PetscBool flg)
636: {
637: PetscFunctionBegin;
640: PetscTryMethod(pc, "PCBDDCSetUseExactDirichlet_C", (PC, PetscBool), (pc, flg));
641: PetscFunctionReturn(PETSC_SUCCESS);
642: }
644: static PetscErrorCode PCBDDCSetLevel_BDDC(PC pc, PetscInt level)
645: {
646: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
648: PetscFunctionBegin;
649: pcbddc->current_level = level;
650: PetscFunctionReturn(PETSC_SUCCESS);
651: }
653: PetscErrorCode PCBDDCSetLevel(PC pc, PetscInt level)
654: {
655: PetscFunctionBegin;
658: PetscTryMethod(pc, "PCBDDCSetLevel_C", (PC, PetscInt), (pc, level));
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }
662: static PetscErrorCode PCBDDCSetLevels_BDDC(PC pc, PetscInt levels)
663: {
664: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
666: PetscFunctionBegin;
667: PetscCheck(levels < PETSC_PCBDDC_MAXLEVELS, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Maximum number of additional levels for BDDC is %d", PETSC_PCBDDC_MAXLEVELS - 1);
668: pcbddc->max_levels = levels;
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: /*@
673: PCBDDCSetLevels - Sets the maximum number of additional levels allowed for multilevel `PCBDDC`
675: Logically Collective
677: Input Parameters:
678: + pc - the preconditioning context
679: - levels - the maximum number of additional levels
681: Options Database Key:
682: . -pc_bddc_levels levels - set the maximum number of additional levels for multilevel BDDC
684: Level: intermediate
686: Note:
687: The default value is 0, which gives the classical two-level BDDC algorithm.
689: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetCoarseningRatio()`
690: @*/
691: PetscErrorCode PCBDDCSetLevels(PC pc, PetscInt levels)
692: {
693: PetscFunctionBegin;
696: PetscTryMethod(pc, "PCBDDCSetLevels_C", (PC, PetscInt), (pc, levels));
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: #define PCBDDC_CUSTOMIZATION_VERSION_LEGACY 0
701: #define PCBDDC_CUSTOMIZATION_VERSION_LATEST 1
702: #define PCBDDC_CUSTOMIZATION_HEADER_SIZE_LEGACY 11
703: #define PCBDDC_CUSTOMIZATION_HEADER_SIZE 32
705: static PetscErrorCode PCBDDCLoadOrSaveCustomization_Private(PC pc, PetscBool load, const char *outfile, PetscInt version)
706: {
707: PetscInt header_storage[PCBDDC_CUSTOMIZATION_HEADER_SIZE] = {0};
708: PetscInt *header;
709: PetscInt nheader;
710: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
711: PetscViewer viewer;
712: MPI_Comm comm = PetscObjectComm((PetscObject)pc);
714: PetscFunctionBegin;
715: if (!load && version == PETSC_DECIDE) version = PCBDDC_CUSTOMIZATION_VERSION_LATEST;
716: if (version == PCBDDC_CUSTOMIZATION_VERSION_LEGACY) {
717: header = header_storage;
718: nheader = PCBDDC_CUSTOMIZATION_HEADER_SIZE_LEGACY;
719: } else {
720: header = header_storage + 1;
721: nheader = PCBDDC_CUSTOMIZATION_HEADER_SIZE;
722: }
723: PetscCall(PetscViewerBinaryOpen(comm, outfile ? outfile : "bddc_dump.dat", load ? FILE_MODE_READ : FILE_MODE_WRITE, &viewer));
724: if (load) {
725: IS is;
726: Mat A;
728: PetscCall(PetscViewerBinaryRead(viewer, header_storage, nheader, NULL, PETSC_INT));
729: if (version == PETSC_DECIDE) version = header_storage[0];
730: PetscCheck(header[0] == 0 || header[0] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
731: PetscCheck(header[1] == 0 || header[1] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
732: PetscCheck(header[2] >= 0, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
733: PetscCheck(header[3] == 0 || header[3] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
734: PetscCheck(header[4] == 0 || header[4] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
735: PetscCheck(header[5] >= 0, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
736: PetscCheck(header[7] == 0 || header[7] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
737: PetscCheck(header[8] == 0 || header[8] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
738: PetscCheck(header[9] == 0 || header[9] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
739: PetscCheck(header[10] == 0 || header[10] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
740: if (version >= 1) {
741: PetscCheck(header[11] == 0 || header[11] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
742: PetscCheck(header[12] == 0 || header[12] == 1, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a BDDC dump next in file");
743: }
744: if (header[0]) {
745: PetscCall(ISCreate(comm, &is));
746: PetscCall(ISLoad(is, viewer));
747: PetscCall(PCBDDCSetDirichletBoundaries(pc, is));
748: PetscCall(ISDestroy(&is));
749: }
750: if (header[1]) {
751: PetscCall(ISCreate(comm, &is));
752: PetscCall(ISLoad(is, viewer));
753: PetscCall(PCBDDCSetNeumannBoundaries(pc, is));
754: PetscCall(ISDestroy(&is));
755: }
756: if (header[2]) {
757: IS *isarray;
759: PetscCall(PetscMalloc1(header[2], &isarray));
760: for (PetscInt i = 0; i < header[2]; i++) {
761: PetscCall(ISCreate(comm, &isarray[i]));
762: PetscCall(ISLoad(isarray[i], viewer));
763: }
764: PetscCall(PCBDDCSetDofsSplitting(pc, header[2], isarray));
765: for (PetscInt i = 0; i < header[2]; i++) PetscCall(ISDestroy(&isarray[i]));
766: PetscCall(PetscFree(isarray));
767: }
768: if (header[3]) {
769: PetscCall(ISCreate(comm, &is));
770: PetscCall(ISLoad(is, viewer));
771: PetscCall(PCBDDCSetPrimalVerticesIS(pc, is));
772: PetscCall(ISDestroy(&is));
773: }
774: if (header[4]) {
775: PetscCall(MatCreate(comm, &A));
776: PetscCall(MatSetType(A, MATAIJ));
777: PetscCall(MatLoad(A, viewer));
778: PetscCall(PCBDDCSetDiscreteGradient(pc, A, header[5], header[6], (PetscBool)header[7], (PetscBool)header[8]));
779: PetscCall(MatDestroy(&A));
780: }
781: if (header[9]) {
782: PetscCall(MatCreate(comm, &A));
783: PetscCall(MatSetType(A, MATIS));
784: PetscCall(MatLoad(A, viewer));
785: PetscCall(PCBDDCSetDivergenceMat(pc, A, (PetscBool)header[10], NULL));
786: PetscCall(MatDestroy(&A));
787: }
788: if (header[11]) {
789: PetscCheck(pcbddc->discretegradient, comm, PETSC_ERR_ARG_CORRUPT, "Missing discrete gradient");
790: PetscCall(ISCreate(comm, &is));
791: PetscCall(ISLoad(is, viewer));
792: PetscCall(PetscObjectCompose((PetscObject)pcbddc->discretegradient, "_elements_corners", (PetscObject)is));
793: PetscCall(ISDestroy(&is));
794: }
795: if (header[12]) {
796: MatNullSpace nsp;
798: PetscCheck(pcbddc->discretegradient, comm, PETSC_ERR_ARG_CORRUPT, "Missing discrete gradient");
799: PetscCall(MatNullSpaceLoad(viewer, &nsp));
800: PetscCall(MatSetNullSpace(pcbddc->discretegradient, nsp));
801: PetscCall(MatNullSpaceDestroy(&nsp));
802: }
803: if (header[13]) {
804: PetscReal *coords;
805: PetscInt cdim, nl;
806: PetscCount nc;
808: PetscCheck(pc->pmat, comm, PETSC_ERR_ORDER, "Need to set the matrix first with PCSetOperators()");
809: PetscCall(PetscLayoutGetLocalSize(pc->pmat->rmap, &nl));
810: cdim = header[13];
811: nc = (PetscCount)cdim * nl;
813: PetscCall(PetscMalloc1(nc, &coords));
814: PetscCall(PetscViewerBinaryReadAll(viewer, coords, nc, PETSC_DECIDE, PETSC_DECIDE, PETSC_REAL));
815: PetscCall(PCSetCoordinates(pc, cdim, nl, coords));
816: PetscCall(PetscFree(coords));
817: }
818: } else {
819: if (version != PCBDDC_CUSTOMIZATION_VERSION_LEGACY) header_storage[0] = version;
820: header[0] = (PetscInt)!!pcbddc->DirichletBoundariesLocal;
821: header[1] = (PetscInt)!!pcbddc->NeumannBoundariesLocal;
822: header[2] = pcbddc->n_ISForDofsLocal;
823: header[3] = (PetscInt)!!pcbddc->user_primal_vertices_local;
824: header[4] = (PetscInt)!!pcbddc->discretegradient;
825: header[5] = pcbddc->nedorder;
826: header[6] = pcbddc->nedfield;
827: header[7] = (PetscInt)pcbddc->nedglobal;
828: header[8] = (PetscInt)pcbddc->conforming;
829: header[9] = (PetscInt)!!pcbddc->divudotp;
830: header[10] = (PetscInt)pcbddc->divudotp_trans;
831: if (header[4]) header[3] = 0;
833: if (version >= 1) {
834: if (pcbddc->discretegradient) {
835: IS is;
836: MatNullSpace nsp;
838: PetscCall(PetscObjectQuery((PetscObject)pcbddc->discretegradient, "_elements_corners", (PetscObject *)&is));
839: header[11] = (PetscBool)!!is;
840: PetscCall(MatGetNullSpace(pcbddc->discretegradient, &nsp));
841: header[12] = (PetscBool)!!nsp;
842: }
843: header[13] = pcbddc->mat_graph->cdim;
844: }
846: PetscCall(PetscViewerBinaryWrite(viewer, header_storage, nheader, PETSC_INT));
847: if (header[0]) PetscCall(PCBDDCViewGlobalIS(pc, pcbddc->DirichletBoundariesLocal, viewer));
848: if (header[1]) PetscCall(PCBDDCViewGlobalIS(pc, pcbddc->NeumannBoundariesLocal, viewer));
849: for (PetscInt i = 0; i < header[2]; i++) PetscCall(PCBDDCViewGlobalIS(pc, pcbddc->ISForDofsLocal[i], viewer));
850: if (header[3]) PetscCall(PCBDDCViewGlobalIS(pc, pcbddc->user_primal_vertices_local, viewer));
851: if (header[4]) PetscCall(MatView(pcbddc->discretegradient, viewer));
852: if (header[9]) PetscCall(MatView(pcbddc->divudotp, viewer));
853: if (header[11]) {
854: IS is;
856: PetscCall(PetscObjectQuery((PetscObject)pcbddc->discretegradient, "_elements_corners", (PetscObject *)&is));
857: PetscCall(ISView(is, viewer));
858: }
859: if (header[12]) {
860: MatNullSpace nsp;
862: PetscCall(MatGetNullSpace(pcbddc->discretegradient, &nsp));
863: PetscCall(MatNullSpaceView(nsp, viewer));
864: }
865: if (header[13]) {
866: PetscReal *coords = pcbddc->mat_graph->coords;
867: PetscCount nc = (PetscCount)pcbddc->mat_graph->cdim * pc->pmat->rmap->n;
869: if (pcbddc->mat_graph->cloc) {
870: Mat_IS *matis = (Mat_IS *)pc->pmat->data;
871: PetscMPIInt cdimi;
872: MPI_Datatype dimrealtype;
874: PetscCall(PetscMalloc1(nc, &coords));
875: PetscCall(PetscMPIIntCast(pcbddc->mat_graph->cdim, &cdimi));
876: PetscCallMPI(MPI_Type_contiguous(cdimi, MPIU_REAL, &dimrealtype));
877: PetscCallMPI(MPI_Type_commit(&dimrealtype));
878: PetscCall(PetscSFReduceBegin(matis->sf, dimrealtype, pcbddc->mat_graph->coords, coords, MPI_REPLACE));
879: PetscCall(PetscSFReduceEnd(matis->sf, dimrealtype, pcbddc->mat_graph->coords, coords, MPI_REPLACE));
880: PetscCallMPI(MPI_Type_free(&dimrealtype));
881: }
882: PetscCall(PetscViewerBinaryWriteAll(viewer, coords, nc, PETSC_DECIDE, PETSC_DECIDE, PETSC_REAL));
883: if (pcbddc->mat_graph->cloc) PetscCall(PetscFree(coords));
884: }
885: }
886: PetscCall(PetscViewerDestroy(&viewer));
887: PetscFunctionReturn(PETSC_SUCCESS);
888: }
890: static PetscErrorCode PCBDDCLoadCustomization_BDDC(PC pc, const char filename[], PetscInt version)
891: {
892: PetscFunctionBegin;
893: PetscCall(PCBDDCLoadOrSaveCustomization_Private(pc, PETSC_TRUE, filename, version));
894: PetscFunctionReturn(PETSC_SUCCESS);
895: }
897: /*@
898: PCBDDCLoadCustomization - Loads user-defined customization data for `PCBDDC` from a binary file
900: Collective
902: Input Parameters:
903: + pc - the preconditioning context
904: . filename - path to the binary file, or `NULL` for `bddc_dump.dat`
905: - version - file format version, or `PETSC_DECIDE` to detect the version from the file
907: Level: advanced
909: Note:
910: This routine is normally called before `PCSetUp()`.
912: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSaveCustomization()`, `PCSetUp()`
913: @*/
914: PetscErrorCode PCBDDCLoadCustomization(PC pc, const char filename[], PetscInt version)
915: {
916: PetscFunctionBegin;
918: if (filename) PetscAssertPointer(filename, 2);
920: PetscTryMethod(pc, "PCBDDCLoadCustomization_C", (PC, const char[], PetscInt), (pc, filename, version));
921: PetscFunctionReturn(PETSC_SUCCESS);
922: }
924: static PetscErrorCode PCBDDCSaveCustomization_BDDC(PC pc, const char filename[], PetscInt version)
925: {
926: PetscFunctionBegin;
927: PetscCall(PCBDDCLoadOrSaveCustomization_Private(pc, PETSC_FALSE, filename, version));
928: PetscFunctionReturn(PETSC_SUCCESS);
929: }
931: /*@
932: PCBDDCSaveCustomization - Saves user-defined customization data for `PCBDDC` to a binary file
934: Collective
936: Input Parameters:
937: + pc - the preconditioning context
938: . filename - path to the binary file, or `NULL` for `bddc_dump.dat`
939: - version - file format version, or `PETSC_DECIDE` to use the latest version
941: Level: advanced
943: Note:
944: Call `PCSetUp()` before this routine so that global customization data has been converted to the local representation stored in the file.
946: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCLoadCustomization()`, `PCSetUp()`
947: @*/
948: PetscErrorCode PCBDDCSaveCustomization(PC pc, const char filename[], PetscInt version)
949: {
950: PetscFunctionBegin;
952: if (filename) PetscAssertPointer(filename, 2);
954: PetscTryMethod(pc, "PCBDDCSaveCustomization_C", (PC, const char[], PetscInt), (pc, filename, version));
955: PetscFunctionReturn(PETSC_SUCCESS);
956: }
958: static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc, IS DirichletBoundaries)
959: {
960: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
961: PetscBool isequal = PETSC_FALSE;
963: PetscFunctionBegin;
964: PetscCall(PetscObjectReference((PetscObject)DirichletBoundaries));
965: if (pcbddc->DirichletBoundaries) PetscCall(ISEqual(DirichletBoundaries, pcbddc->DirichletBoundaries, &isequal));
966: /* last user setting takes precedence -> destroy any other customization */
967: PetscCall(ISDestroy(&pcbddc->DirichletBoundariesLocal));
968: PetscCall(ISDestroy(&pcbddc->DirichletBoundaries));
969: pcbddc->DirichletBoundaries = DirichletBoundaries;
970: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
971: PetscFunctionReturn(PETSC_SUCCESS);
972: }
974: /*@
975: PCBDDCSetDirichletBoundaries - Sets the Dirichlet boundary degrees of freedom in global numbering
977: Collective
979: Input Parameters:
980: + pc - the preconditioning context
981: - DirichletBoundaries - index set of Dirichlet boundary degrees of freedom in global numbering
983: Level: intermediate
985: Note:
986: Provide this information when Dirichlet conditions have been imposed with `MatZeroRows()` or `MatZeroRowsColumns()`.
987: Any process can list any global degree of freedom.
989: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCGetDirichletBoundaries()`, `PCBDDCSetDirichletBoundariesLocal()`, `MatZeroRows()`, `MatZeroRowsColumns()`
990: @*/
991: PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc, IS DirichletBoundaries)
992: {
993: PetscFunctionBegin;
996: PetscCheckSameComm(pc, 1, DirichletBoundaries, 2);
997: PetscTryMethod(pc, "PCBDDCSetDirichletBoundaries_C", (PC, IS), (pc, DirichletBoundaries));
998: PetscFunctionReturn(PETSC_SUCCESS);
999: }
1001: static PetscErrorCode PCBDDCSetDirichletBoundariesLocal_BDDC(PC pc, IS DirichletBoundaries)
1002: {
1003: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1004: PetscBool isequal = PETSC_FALSE;
1006: PetscFunctionBegin;
1007: PetscCall(PetscObjectReference((PetscObject)DirichletBoundaries));
1008: if (pcbddc->DirichletBoundariesLocal) PetscCall(ISEqual(DirichletBoundaries, pcbddc->DirichletBoundariesLocal, &isequal));
1009: /* last user setting takes precedence -> destroy any other customization */
1010: PetscCall(ISDestroy(&pcbddc->DirichletBoundariesLocal));
1011: PetscCall(ISDestroy(&pcbddc->DirichletBoundaries));
1012: pcbddc->DirichletBoundariesLocal = DirichletBoundaries;
1013: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1014: PetscFunctionReturn(PETSC_SUCCESS);
1015: }
1017: /*@
1018: PCBDDCSetDirichletBoundariesLocal - Sets the Dirichlet boundary degrees of freedom in local numbering
1020: Collective
1022: Input Parameters:
1023: + pc - the preconditioning context
1024: - DirichletBoundaries - index set of Dirichlet boundary degrees of freedom in the numbering of the local `MATIS` matrix
1026: Level: intermediate
1028: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCGetDirichletBoundariesLocal()`, `PCBDDCSetDirichletBoundaries()`, `MatZeroRows()`, `MatZeroRowsColumns()`
1029: @*/
1030: PetscErrorCode PCBDDCSetDirichletBoundariesLocal(PC pc, IS DirichletBoundaries)
1031: {
1032: PetscFunctionBegin;
1035: PetscCheckSameComm(pc, 1, DirichletBoundaries, 2);
1036: PetscTryMethod(pc, "PCBDDCSetDirichletBoundariesLocal_C", (PC, IS), (pc, DirichletBoundaries));
1037: PetscFunctionReturn(PETSC_SUCCESS);
1038: }
1040: static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc, IS NeumannBoundaries)
1041: {
1042: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1043: PetscBool isequal = PETSC_FALSE;
1045: PetscFunctionBegin;
1046: PetscCall(PetscObjectReference((PetscObject)NeumannBoundaries));
1047: if (pcbddc->NeumannBoundaries) PetscCall(ISEqual(NeumannBoundaries, pcbddc->NeumannBoundaries, &isequal));
1048: /* last user setting takes precedence -> destroy any other customization */
1049: PetscCall(ISDestroy(&pcbddc->NeumannBoundariesLocal));
1050: PetscCall(ISDestroy(&pcbddc->NeumannBoundaries));
1051: pcbddc->NeumannBoundaries = NeumannBoundaries;
1052: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1053: PetscFunctionReturn(PETSC_SUCCESS);
1054: }
1056: /*@
1057: PCBDDCSetNeumannBoundaries - Sets the Neumann boundary degrees of freedom in global numbering
1059: Collective
1061: Input Parameters:
1062: + pc - the preconditioning context
1063: - NeumannBoundaries - index set of Neumann boundary degrees of freedom in global numbering
1065: Level: intermediate
1067: Note:
1068: Any process can list any global degree of freedom.
1070: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCGetNeumannBoundaries()`, `PCBDDCSetNeumannBoundariesLocal()`
1071: @*/
1072: PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc, IS NeumannBoundaries)
1073: {
1074: PetscFunctionBegin;
1077: PetscCheckSameComm(pc, 1, NeumannBoundaries, 2);
1078: PetscTryMethod(pc, "PCBDDCSetNeumannBoundaries_C", (PC, IS), (pc, NeumannBoundaries));
1079: PetscFunctionReturn(PETSC_SUCCESS);
1080: }
1082: static PetscErrorCode PCBDDCSetNeumannBoundariesLocal_BDDC(PC pc, IS NeumannBoundaries)
1083: {
1084: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1085: PetscBool isequal = PETSC_FALSE;
1087: PetscFunctionBegin;
1088: PetscCall(PetscObjectReference((PetscObject)NeumannBoundaries));
1089: if (pcbddc->NeumannBoundariesLocal) PetscCall(ISEqual(NeumannBoundaries, pcbddc->NeumannBoundariesLocal, &isequal));
1090: /* last user setting takes precedence -> destroy any other customization */
1091: PetscCall(ISDestroy(&pcbddc->NeumannBoundariesLocal));
1092: PetscCall(ISDestroy(&pcbddc->NeumannBoundaries));
1093: pcbddc->NeumannBoundariesLocal = NeumannBoundaries;
1094: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1095: PetscFunctionReturn(PETSC_SUCCESS);
1096: }
1098: /*@
1099: PCBDDCSetNeumannBoundariesLocal - Sets the Neumann boundary degrees of freedom in local numbering
1101: Collective
1103: Input Parameters:
1104: + pc - the preconditioning context
1105: - NeumannBoundaries - index set of Neumann boundary degrees of freedom in the numbering of the local `MATIS` matrix
1107: Level: intermediate
1109: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCGetNeumannBoundariesLocal()`, `PCBDDCSetNeumannBoundaries()`
1110: @*/
1111: PetscErrorCode PCBDDCSetNeumannBoundariesLocal(PC pc, IS NeumannBoundaries)
1112: {
1113: PetscFunctionBegin;
1116: PetscCheckSameComm(pc, 1, NeumannBoundaries, 2);
1117: PetscTryMethod(pc, "PCBDDCSetNeumannBoundariesLocal_C", (PC, IS), (pc, NeumannBoundaries));
1118: PetscFunctionReturn(PETSC_SUCCESS);
1119: }
1121: static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc, IS *DirichletBoundaries)
1122: {
1123: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1125: PetscFunctionBegin;
1126: *DirichletBoundaries = pcbddc->DirichletBoundaries;
1127: PetscFunctionReturn(PETSC_SUCCESS);
1128: }
1130: /*@
1131: PCBDDCGetDirichletBoundaries - Gets the Dirichlet boundary degrees of freedom in global numbering
1133: Not Collective
1135: Input Parameter:
1136: . pc - the preconditioning context
1138: Output Parameter:
1139: . DirichletBoundaries - index set of Dirichlet boundary degrees of freedom in global numbering, or `NULL` if not set
1141: Level: intermediate
1143: Notes:
1144: The returned `IS`, if any, is the one supplied to `PCBDDCSetDirichletBoundaries()`.
1146: The returned `IS` is owned by `pc`; the caller must not destroy it.
1148: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDirichletBoundaries()`, `PCBDDCGetDirichletBoundariesLocal()`
1149: @*/
1150: PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc, IS *DirichletBoundaries)
1151: {
1152: PetscFunctionBegin;
1154: PetscUseMethod(pc, "PCBDDCGetDirichletBoundaries_C", (PC, IS *), (pc, DirichletBoundaries));
1155: PetscFunctionReturn(PETSC_SUCCESS);
1156: }
1158: static PetscErrorCode PCBDDCGetDirichletBoundariesLocal_BDDC(PC pc, IS *DirichletBoundaries)
1159: {
1160: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1162: PetscFunctionBegin;
1163: *DirichletBoundaries = pcbddc->DirichletBoundariesLocal;
1164: PetscFunctionReturn(PETSC_SUCCESS);
1165: }
1167: /*@
1168: PCBDDCGetDirichletBoundariesLocal - Gets the Dirichlet boundary degrees of freedom in local numbering
1170: Not Collective
1172: Input Parameter:
1173: . pc - the preconditioning context
1175: Output Parameter:
1176: . DirichletBoundaries - index set of Dirichlet boundary degrees of freedom in the numbering of the local `MATIS` matrix, or `NULL` if unavailable
1178: Level: intermediate
1180: Notes:
1181: The index set is supplied by `PCBDDCSetDirichletBoundariesLocal()` or obtained from `PCBDDCSetDirichletBoundaries()` during `PCSetUp()`.
1182: Setup may update the local index set to make the boundary information consistent across shared degrees of freedom.
1184: The returned `IS` is owned by `pc`; the caller must not destroy it.
1186: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDirichletBoundariesLocal()`, `PCBDDCGetDirichletBoundaries()`, `PCBDDCSetDirichletBoundaries()`
1187: @*/
1188: PetscErrorCode PCBDDCGetDirichletBoundariesLocal(PC pc, IS *DirichletBoundaries)
1189: {
1190: PetscFunctionBegin;
1192: PetscUseMethod(pc, "PCBDDCGetDirichletBoundariesLocal_C", (PC, IS *), (pc, DirichletBoundaries));
1193: PetscFunctionReturn(PETSC_SUCCESS);
1194: }
1196: static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc, IS *NeumannBoundaries)
1197: {
1198: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1200: PetscFunctionBegin;
1201: *NeumannBoundaries = pcbddc->NeumannBoundaries;
1202: PetscFunctionReturn(PETSC_SUCCESS);
1203: }
1205: /*@
1206: PCBDDCGetNeumannBoundaries - Gets the Neumann boundary degrees of freedom in global numbering
1208: Not Collective
1210: Input Parameter:
1211: . pc - the preconditioning context
1213: Output Parameter:
1214: . NeumannBoundaries - index set of Neumann boundary degrees of freedom in global numbering, or `NULL` if not set
1216: Level: intermediate
1218: Notes:
1219: The returned `IS`, if any, is the one supplied to `PCBDDCSetNeumannBoundaries()`.
1221: The returned `IS` is owned by `pc`; the caller must not destroy it.
1223: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCGetNeumannBoundariesLocal()`
1224: @*/
1225: PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc, IS *NeumannBoundaries)
1226: {
1227: PetscFunctionBegin;
1229: PetscUseMethod(pc, "PCBDDCGetNeumannBoundaries_C", (PC, IS *), (pc, NeumannBoundaries));
1230: PetscFunctionReturn(PETSC_SUCCESS);
1231: }
1233: static PetscErrorCode PCBDDCGetNeumannBoundariesLocal_BDDC(PC pc, IS *NeumannBoundaries)
1234: {
1235: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1237: PetscFunctionBegin;
1238: *NeumannBoundaries = pcbddc->NeumannBoundariesLocal;
1239: PetscFunctionReturn(PETSC_SUCCESS);
1240: }
1242: /*@
1243: PCBDDCGetNeumannBoundariesLocal - Gets the Neumann boundary degrees of freedom in local numbering
1245: Not Collective
1247: Input Parameter:
1248: . pc - the preconditioning context
1250: Output Parameter:
1251: . NeumannBoundaries - index set of Neumann boundary degrees of freedom in the numbering of the local `MATIS` matrix, or `NULL` if unavailable
1253: Level: intermediate
1255: Notes:
1256: The index set is supplied by `PCBDDCSetNeumannBoundariesLocal()` or obtained from `PCBDDCSetNeumannBoundaries()` during `PCSetUp()`.
1257: Setup may update the local index set to make the boundary information consistent across shared degrees of freedom.
1259: The returned `IS` is owned by `pc`; the caller must not destroy it.
1261: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCSetNeumannBoundariesLocal()`, `PCBDDCGetNeumannBoundaries()`
1262: @*/
1263: PetscErrorCode PCBDDCGetNeumannBoundariesLocal(PC pc, IS *NeumannBoundaries)
1264: {
1265: PetscFunctionBegin;
1267: PetscUseMethod(pc, "PCBDDCGetNeumannBoundariesLocal_C", (PC, IS *), (pc, NeumannBoundaries));
1268: PetscFunctionReturn(PETSC_SUCCESS);
1269: }
1271: static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs, const PetscInt xadj[], const PetscInt adjncy[], PetscCopyMode copymode)
1272: {
1273: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1274: PCBDDCGraph mat_graph = pcbddc->mat_graph;
1275: PetscBool same_data = PETSC_FALSE;
1277: PetscFunctionBegin;
1278: if (!nvtxs) {
1279: if (copymode == PETSC_OWN_POINTER) {
1280: PetscCall(PetscFree(xadj));
1281: PetscCall(PetscFree(adjncy));
1282: }
1283: PetscCall(PCBDDCGraphResetCSR(mat_graph));
1284: PetscFunctionReturn(PETSC_SUCCESS);
1285: }
1286: if (mat_graph->nvtxs == nvtxs && mat_graph->freecsr) { /* we own the data */
1287: if (mat_graph->xadj == xadj && mat_graph->adjncy == adjncy) same_data = PETSC_TRUE;
1288: if (!same_data && mat_graph->xadj[nvtxs] == xadj[nvtxs]) {
1289: PetscCall(PetscArraycmp(xadj, mat_graph->xadj, nvtxs + 1, &same_data));
1290: if (same_data) PetscCall(PetscArraycmp(adjncy, mat_graph->adjncy, xadj[nvtxs], &same_data));
1291: }
1292: }
1293: if (!same_data) {
1294: /* free old CSR */
1295: PetscCall(PCBDDCGraphResetCSR(mat_graph));
1296: /* get CSR into graph structure */
1297: if (copymode == PETSC_COPY_VALUES) {
1298: PetscCall(PetscMalloc1(nvtxs + 1, &mat_graph->xadj));
1299: PetscCall(PetscMalloc1(xadj[nvtxs], &mat_graph->adjncy));
1300: PetscCall(PetscArraycpy(mat_graph->xadj, xadj, nvtxs + 1));
1301: PetscCall(PetscArraycpy(mat_graph->adjncy, adjncy, xadj[nvtxs]));
1302: mat_graph->freecsr = PETSC_TRUE;
1303: } else if (copymode == PETSC_OWN_POINTER) {
1304: mat_graph->xadj = (PetscInt *)xadj;
1305: mat_graph->adjncy = (PetscInt *)adjncy;
1306: mat_graph->freecsr = PETSC_TRUE;
1307: } else if (copymode == PETSC_USE_POINTER) {
1308: mat_graph->xadj = (PetscInt *)xadj;
1309: mat_graph->adjncy = (PetscInt *)adjncy;
1310: mat_graph->freecsr = PETSC_FALSE;
1311: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported copy mode %d", copymode);
1312: mat_graph->nvtxs_csr = nvtxs;
1313: pcbddc->recompute_topography = PETSC_TRUE;
1314: }
1315: PetscFunctionReturn(PETSC_SUCCESS);
1316: }
1318: /*@
1319: PCBDDCSetLocalAdjacencyGraph - Sets the adjacency graph of the local degrees of freedom in CSR format
1321: Not Collective
1323: Input Parameters:
1324: + pc - the preconditioning context
1325: . nvtxs - number of local graph vertices, equal to the number of local degrees of freedom
1326: . xadj - CSR row offsets, of length `nvtxs` + 1
1327: . adjncy - CSR column indices, of length `xadj[nvtxs]`
1328: - copymode - `PETSC_COPY_VALUES`, `PETSC_USE_POINTER`, or `PETSC_OWN_POINTER`
1330: Level: intermediate
1332: Notes:
1333: Graph vertices use the numbering of the local `MATIS` matrix. The CSR arrays use zero-based indexing.
1334: A degree of freedom `i` is considered connected to all local degrees of freedom if `xadj[i+1] - xadj[i] == 1`
1335: and `adjncy[xadj[i]]` is negative.
1337: With `PETSC_COPY_VALUES`, the arrays are copied. With `PETSC_OWN_POINTER`, ownership of the arrays is transferred to `pc`;
1338: the arrays must have been allocated with `PetscMalloc()`. With `PETSC_USE_POINTER`, the caller must keep the arrays valid
1339: and unchanged until the graph is replaced or the preconditioner is reset or destroyed.
1341: Passing `nvtxs` equal to 0 clears the user-defined graph.
1343: .seealso: [](ch_ksp), `PCBDDC`, `PetscCopyMode`
1344: @*/
1345: PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc, PetscInt nvtxs, const PetscInt xadj[], const PetscInt adjncy[], PetscCopyMode copymode)
1346: {
1347: PetscBool f = PETSC_FALSE;
1349: PetscFunctionBegin;
1351: if (nvtxs) {
1352: PetscAssertPointer(xadj, 3);
1353: if (xadj[nvtxs]) PetscAssertPointer(adjncy, 4);
1354: }
1355: PetscTryMethod(pc, "PCBDDCSetLocalAdjacencyGraph_C", (PC, PetscInt, const PetscInt[], const PetscInt[], PetscCopyMode), (pc, nvtxs, xadj, adjncy, copymode));
1356: /* free arrays if PCBDDC is not the PC type */
1357: PetscCall(PetscObjectHasFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", &f));
1358: if (!f && copymode == PETSC_OWN_POINTER) {
1359: PetscCall(PetscFree(xadj));
1360: PetscCall(PetscFree(adjncy));
1361: }
1362: PetscFunctionReturn(PETSC_SUCCESS);
1363: }
1365: static PetscErrorCode PCBDDCSetDofsSplittingLocal_BDDC(PC pc, PetscInt n_is, IS ISForDofs[])
1366: {
1367: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1368: PetscInt i;
1369: PetscBool isequal = PETSC_FALSE;
1371: PetscFunctionBegin;
1372: if (pcbddc->n_ISForDofsLocal == n_is) {
1373: for (i = 0; i < n_is; i++) {
1374: PetscBool isequalt;
1375: PetscCall(ISEqual(ISForDofs[i], pcbddc->ISForDofsLocal[i], &isequalt));
1376: if (!isequalt) break;
1377: }
1378: if (i == n_is) isequal = PETSC_TRUE;
1379: }
1380: for (i = 0; i < n_is; i++) PetscCall(PetscObjectReference((PetscObject)ISForDofs[i]));
1381: /* Destroy ISes if they were already set */
1382: for (i = 0; i < pcbddc->n_ISForDofsLocal; i++) PetscCall(ISDestroy(&pcbddc->ISForDofsLocal[i]));
1383: PetscCall(PetscFree(pcbddc->ISForDofsLocal));
1384: /* last user setting takes precedence -> destroy any other customization */
1385: for (i = 0; i < pcbddc->n_ISForDofs; i++) PetscCall(ISDestroy(&pcbddc->ISForDofs[i]));
1386: PetscCall(PetscFree(pcbddc->ISForDofs));
1387: pcbddc->n_ISForDofs = 0;
1388: /* allocate space then set */
1389: if (n_is) PetscCall(PetscMalloc1(n_is, &pcbddc->ISForDofsLocal));
1390: for (i = 0; i < n_is; i++) pcbddc->ISForDofsLocal[i] = ISForDofs[i];
1391: pcbddc->n_ISForDofsLocal = n_is;
1392: if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
1393: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1394: PetscFunctionReturn(PETSC_SUCCESS);
1395: }
1397: /*@
1398: PCBDDCSetDofsSplittingLocal - Sets the fields of the local subdomain matrix
1400: Collective
1402: Input Parameters:
1403: + pc - the preconditioning context
1404: . n_is - number of index sets defining the fields, must be the same on all MPI processes
1405: - ISForDofs - array of `n_is` index sets describing the fields in the numbering of the local `MATIS` matrix
1407: Level: intermediate
1409: Note:
1410: Degrees of freedom not listed in any index set belong to the complement field.
1412: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDofsSplitting()`
1413: @*/
1414: PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc, PetscInt n_is, IS ISForDofs[])
1415: {
1416: PetscInt i;
1418: PetscFunctionBegin;
1421: for (i = 0; i < n_is; i++) {
1422: PetscCheckSameComm(pc, 1, ISForDofs[i], 3);
1424: }
1425: PetscTryMethod(pc, "PCBDDCSetDofsSplittingLocal_C", (PC, PetscInt, IS[]), (pc, n_is, ISForDofs));
1426: PetscFunctionReturn(PETSC_SUCCESS);
1427: }
1429: static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc, PetscInt n_is, IS ISForDofs[])
1430: {
1431: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1432: PetscInt i;
1433: PetscBool isequal = PETSC_FALSE;
1435: PetscFunctionBegin;
1436: if (pcbddc->n_ISForDofs == n_is) {
1437: for (i = 0; i < n_is; i++) {
1438: PetscBool isequalt;
1439: PetscCall(ISEqual(ISForDofs[i], pcbddc->ISForDofs[i], &isequalt));
1440: if (!isequalt) break;
1441: }
1442: if (i == n_is) isequal = PETSC_TRUE;
1443: }
1444: for (i = 0; i < n_is; i++) PetscCall(PetscObjectReference((PetscObject)ISForDofs[i]));
1445: /* Destroy ISes if they were already set */
1446: for (i = 0; i < pcbddc->n_ISForDofs; i++) PetscCall(ISDestroy(&pcbddc->ISForDofs[i]));
1447: PetscCall(PetscFree(pcbddc->ISForDofs));
1448: /* last user setting takes precedence -> destroy any other customization */
1449: for (i = 0; i < pcbddc->n_ISForDofsLocal; i++) PetscCall(ISDestroy(&pcbddc->ISForDofsLocal[i]));
1450: PetscCall(PetscFree(pcbddc->ISForDofsLocal));
1451: pcbddc->n_ISForDofsLocal = 0;
1452: /* allocate space then set */
1453: if (n_is) PetscCall(PetscMalloc1(n_is, &pcbddc->ISForDofs));
1454: for (i = 0; i < n_is; i++) pcbddc->ISForDofs[i] = ISForDofs[i];
1455: pcbddc->n_ISForDofs = n_is;
1456: if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
1457: if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1458: PetscFunctionReturn(PETSC_SUCCESS);
1459: }
1461: /*@
1462: PCBDDCSetDofsSplitting - Sets the fields of the global matrix
1464: Collective
1466: Input Parameters:
1467: + pc - the preconditioning context
1468: . n_is - number of index sets defining the fields, must be the same on all MPI processes
1469: - ISForDofs - array of `n_is` index sets describing the fields in global numbering
1471: Level: intermediate
1473: Note:
1474: Any process can list any global degree of freedom. Degrees of freedom not listed in any index set belong to the complement field.
1476: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDofsSplittingLocal()`
1477: @*/
1478: PetscErrorCode PCBDDCSetDofsSplitting(PC pc, PetscInt n_is, IS ISForDofs[])
1479: {
1480: PetscInt i;
1482: PetscFunctionBegin;
1485: for (i = 0; i < n_is; i++) {
1487: PetscCheckSameComm(pc, 1, ISForDofs[i], 3);
1488: }
1489: PetscTryMethod(pc, "PCBDDCSetDofsSplitting_C", (PC, PetscInt, IS[]), (pc, n_is, ISForDofs));
1490: PetscFunctionReturn(PETSC_SUCCESS);
1491: }
1493: static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
1494: {
1495: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1496: PC_IS *pcis = (PC_IS *)pc->data;
1497: Vec used_vec;
1498: PetscBool iscg, save_rhs = PETSC_TRUE, benign_correction_computed;
1500: PetscFunctionBegin;
1501: /* if we are working with CG, one dirichlet solve can be avoided during Krylov iterations */
1502: if (ksp) {
1503: PetscCall(PetscObjectTypeCompareAny((PetscObject)ksp, &iscg, KSPCG, KSPGROPPCG, KSPPIPECG, KSPPIPELCG, KSPPIPECGRR, ""));
1504: if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || !iscg || pc->mat != pc->pmat) PetscCall(PCBDDCSetUseExactDirichlet(pc, PETSC_FALSE));
1505: }
1506: if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || pc->mat != pc->pmat) PetscCall(PCBDDCSetUseExactDirichlet(pc, PETSC_FALSE));
1508: /* Creates parallel work vectors used in presolve */
1509: if (!pcbddc->original_rhs) PetscCall(VecDuplicate(pcis->vec1_global, &pcbddc->original_rhs));
1510: if (!pcbddc->temp_solution) PetscCall(VecDuplicate(pcis->vec1_global, &pcbddc->temp_solution));
1512: pcbddc->temp_solution_used = PETSC_FALSE;
1513: if (x) {
1514: PetscCall(PetscObjectReference((PetscObject)x));
1515: used_vec = x;
1516: } else { /* it can only happen when calling PCBDDCMatFETIDPGetRHS */
1517: PetscCall(PetscObjectReference((PetscObject)pcbddc->temp_solution));
1518: used_vec = pcbddc->temp_solution;
1519: PetscCall(VecSet(used_vec, 0.0));
1520: pcbddc->temp_solution_used = PETSC_TRUE;
1521: PetscCall(VecCopy(rhs, pcbddc->original_rhs));
1522: save_rhs = PETSC_FALSE;
1523: pcbddc->eliminate_dirdofs = PETSC_TRUE;
1524: }
1526: /* hack into ksp data structure since PCPreSolve comes earlier than setting to zero the guess in src/ksp/ksp/interface/itfunc.c */
1527: if (ksp) {
1528: /* store the flag for the initial guess since it will be restored back during PCPostSolve_BDDC */
1529: PetscCall(KSPGetInitialGuessNonzero(ksp, &pcbddc->ksp_guess_nonzero));
1530: if (!pcbddc->ksp_guess_nonzero) PetscCall(VecSet(used_vec, 0.0));
1531: }
1533: pcbddc->rhs_change = PETSC_FALSE;
1534: /* Take into account zeroed rows -> change rhs and store solution removed */
1535: if (rhs && pcbddc->eliminate_dirdofs) {
1536: IS dirIS = NULL;
1538: /* DirichletBoundariesLocal may not be consistent among neighbours; gets a dirichlet dofs IS from graph (may be cached) */
1539: PetscCall(PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph, &dirIS));
1540: if (dirIS) {
1541: Mat_IS *matis = (Mat_IS *)pc->pmat->data;
1542: PetscInt dirsize, i, *is_indices;
1543: PetscScalar *array_x;
1544: const PetscScalar *array_diagonal;
1546: PetscCall(MatGetDiagonal(pc->pmat, pcis->vec1_global));
1547: PetscCall(VecPointwiseDivide(pcis->vec1_global, rhs, pcis->vec1_global));
1548: PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_global, pcis->vec2_N, INSERT_VALUES, SCATTER_FORWARD));
1549: PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_global, pcis->vec2_N, INSERT_VALUES, SCATTER_FORWARD));
1550: PetscCall(VecScatterBegin(matis->rctx, used_vec, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD));
1551: PetscCall(VecScatterEnd(matis->rctx, used_vec, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD));
1552: PetscCall(ISGetLocalSize(dirIS, &dirsize));
1553: PetscCall(VecGetArray(pcis->vec1_N, &array_x));
1554: PetscCall(VecGetArrayRead(pcis->vec2_N, &array_diagonal));
1555: PetscCall(ISGetIndices(dirIS, (const PetscInt **)&is_indices));
1556: for (i = 0; i < dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]];
1557: PetscCall(ISRestoreIndices(dirIS, (const PetscInt **)&is_indices));
1558: PetscCall(VecRestoreArrayRead(pcis->vec2_N, &array_diagonal));
1559: PetscCall(VecRestoreArray(pcis->vec1_N, &array_x));
1560: PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_N, used_vec, INSERT_VALUES, SCATTER_REVERSE));
1561: PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_N, used_vec, INSERT_VALUES, SCATTER_REVERSE));
1562: pcbddc->rhs_change = PETSC_TRUE;
1563: PetscCall(ISDestroy(&dirIS));
1564: }
1565: }
1567: /* remove the computed solution or the initial guess from the rhs */
1568: if (pcbddc->rhs_change || (ksp && pcbddc->ksp_guess_nonzero)) {
1569: /* save the original rhs */
1570: if (save_rhs) {
1571: PetscCall(VecSwap(rhs, pcbddc->original_rhs));
1572: save_rhs = PETSC_FALSE;
1573: }
1574: pcbddc->rhs_change = PETSC_TRUE;
1575: PetscCall(VecScale(used_vec, -1.0));
1576: PetscCall(MatMultAdd(pc->mat, used_vec, pcbddc->original_rhs, rhs));
1577: PetscCall(VecScale(used_vec, -1.0));
1578: PetscCall(VecCopy(used_vec, pcbddc->temp_solution));
1579: pcbddc->temp_solution_used = PETSC_TRUE;
1580: if (ksp) PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_FALSE));
1581: }
1582: PetscCall(VecDestroy(&used_vec));
1584: /* compute initial vector in benign space if needed
1585: and remove non-benign solution from the rhs */
1586: benign_correction_computed = PETSC_FALSE;
1587: if (rhs && pcbddc->benign_compute_correction && (pcbddc->benign_have_null || pcbddc->benign_apply_coarse_only)) {
1588: /* compute u^*_h using ideas similar to those in Xuemin Tu's PhD thesis (see Section 4.8.1)
1589: Recursively apply BDDC in the multilevel case */
1590: if (!pcbddc->benign_vec) PetscCall(VecDuplicate(rhs, &pcbddc->benign_vec));
1591: /* keep applying coarse solver unless we no longer have benign subdomains */
1592: pcbddc->benign_apply_coarse_only = pcbddc->benign_have_null ? PETSC_TRUE : PETSC_FALSE;
1593: if (!pcbddc->benign_skip_correction) {
1594: PetscCall(PCApply_BDDC(pc, rhs, pcbddc->benign_vec));
1595: benign_correction_computed = PETSC_TRUE;
1596: if (pcbddc->temp_solution_used) PetscCall(VecAXPY(pcbddc->temp_solution, 1.0, pcbddc->benign_vec));
1597: PetscCall(VecScale(pcbddc->benign_vec, -1.0));
1598: /* store the original rhs if not done earlier */
1599: if (save_rhs) PetscCall(VecSwap(rhs, pcbddc->original_rhs));
1600: if (pcbddc->rhs_change) {
1601: PetscCall(MatMultAdd(pc->mat, pcbddc->benign_vec, rhs, rhs));
1602: } else {
1603: PetscCall(MatMultAdd(pc->mat, pcbddc->benign_vec, pcbddc->original_rhs, rhs));
1604: }
1605: pcbddc->rhs_change = PETSC_TRUE;
1606: }
1607: pcbddc->benign_apply_coarse_only = PETSC_FALSE;
1608: } else {
1609: PetscCall(VecDestroy(&pcbddc->benign_vec));
1610: }
1612: /* dbg output */
1613: if (pcbddc->dbg_flag && benign_correction_computed) {
1614: Vec v;
1616: PetscCall(VecDuplicate(pcis->vec1_global, &v));
1617: if (pcbddc->ChangeOfBasisMatrix) {
1618: PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, rhs, v));
1619: } else {
1620: PetscCall(VecCopy(rhs, v));
1621: }
1622: PetscCall(PCBDDCBenignGetOrSetP0(pc, v, PETSC_TRUE));
1623: PetscCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "LEVEL %" PetscInt_FMT ": is the correction benign?\n", pcbddc->current_level));
1624: PetscCall(PetscScalarView(pcbddc->benign_n, pcbddc->benign_p0, pcbddc->dbg_viewer));
1625: PetscCall(PetscViewerFlush(pcbddc->dbg_viewer));
1626: PetscCall(VecDestroy(&v));
1627: }
1629: /* set initial guess if using PCG */
1630: pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
1631: if (x && pcbddc->use_exact_dirichlet_trick) {
1632: PetscCall(VecSet(x, 0.0));
1633: if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) {
1634: if (benign_correction_computed) { /* we have already saved the changed rhs */
1635: PetscCall(VecLockReadPop(pcis->vec1_global));
1636: } else {
1637: PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, rhs, pcis->vec1_global));
1638: }
1639: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec1_global, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1640: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec1_global, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1641: } else {
1642: PetscCall(VecScatterBegin(pcis->global_to_D, rhs, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1643: PetscCall(VecScatterEnd(pcis->global_to_D, rhs, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1644: }
1645: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1646: PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D));
1647: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1648: PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D));
1649: if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) {
1650: PetscCall(VecSet(pcis->vec1_global, 0.));
1651: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, pcis->vec1_global, INSERT_VALUES, SCATTER_REVERSE));
1652: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, pcis->vec1_global, INSERT_VALUES, SCATTER_REVERSE));
1653: PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcis->vec1_global, x));
1654: } else {
1655: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, x, INSERT_VALUES, SCATTER_REVERSE));
1656: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, x, INSERT_VALUES, SCATTER_REVERSE));
1657: }
1658: if (ksp) PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
1659: pcbddc->exact_dirichlet_trick_app = PETSC_TRUE;
1660: } else if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior && benign_correction_computed && pcbddc->use_exact_dirichlet_trick) {
1661: PetscCall(VecLockReadPop(pcis->vec1_global));
1662: }
1663: PetscFunctionReturn(PETSC_SUCCESS);
1664: }
1666: static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
1667: {
1668: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1670: PetscFunctionBegin;
1671: /* add solution removed in presolve */
1672: if (x && pcbddc->rhs_change) {
1673: if (pcbddc->temp_solution_used) PetscCall(VecAXPY(x, 1.0, pcbddc->temp_solution));
1674: else if (pcbddc->benign_compute_correction && pcbddc->benign_vec) PetscCall(VecAXPY(x, -1.0, pcbddc->benign_vec));
1675: /* restore to original state (not for FETI-DP) */
1676: if (ksp) pcbddc->temp_solution_used = PETSC_FALSE;
1677: }
1679: /* restore rhs to its original state (not needed for FETI-DP) */
1680: if (rhs && pcbddc->rhs_change) {
1681: PetscCall(VecSwap(rhs, pcbddc->original_rhs));
1682: pcbddc->rhs_change = PETSC_FALSE;
1683: }
1684: /* restore ksp guess state */
1685: if (ksp) {
1686: PetscCall(KSPSetInitialGuessNonzero(ksp, pcbddc->ksp_guess_nonzero));
1687: /* reset flag for exact dirichlet trick */
1688: pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
1689: }
1690: PetscFunctionReturn(PETSC_SUCCESS);
1691: }
1693: static PetscErrorCode PCSetUp_BDDC(PC pc)
1694: {
1695: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1696: PCBDDCSubSchurs sub_schurs;
1697: Mat_IS *matis;
1698: MatNullSpace nearnullspace;
1699: Mat lA;
1700: IS lP, zerodiag = NULL;
1701: PetscInt nrows, ncols;
1702: PetscMPIInt size;
1703: PetscBool computesubschurs;
1704: PetscBool computeconstraintsmatrix;
1705: PetscBool new_nearnullspace_provided, ismatis;
1706: PetscBool isset, issym, isspd;
1708: PetscFunctionBegin;
1709: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATIS, &ismatis));
1710: PetscCheck(ismatis, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PCBDDC preconditioner requires matrix of type MATIS");
1711: PetscCall(MatGetSize(pc->pmat, &nrows, &ncols));
1712: PetscCheck(nrows == ncols, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCBDDC preconditioner requires a square matrix for constructing the preconditioner");
1713: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
1715: matis = (Mat_IS *)pc->pmat->data;
1716: /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */
1717: /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetUp
1718: Also, BDDC builds its own KSP for the Dirichlet problem */
1719: if (!pc->setupcalled || pc->flag == DIFFERENT_NONZERO_PATTERN) pcbddc->recompute_topography = PETSC_TRUE;
1720: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &pcbddc->recompute_topography, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)pc)));
1721: if (pcbddc->recompute_topography) {
1722: pcbddc->graphanalyzed = PETSC_FALSE;
1723: computeconstraintsmatrix = PETSC_TRUE;
1724: } else {
1725: computeconstraintsmatrix = PETSC_FALSE;
1726: }
1728: /* check parameters' compatibility */
1729: if (!pcbddc->use_deluxe_scaling) pcbddc->deluxe_zerorows = PETSC_FALSE;
1730: pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_threshold[0] != 0.0 || pcbddc->adaptive_threshold[1] != 0.0);
1731: pcbddc->use_deluxe_scaling = (PetscBool)(pcbddc->use_deluxe_scaling && (size > 1 || matis->allow_repeated));
1732: pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_selection && (size > 1 || matis->allow_repeated));
1733: pcbddc->adaptive_userdefined = (PetscBool)(pcbddc->adaptive_selection && pcbddc->adaptive_userdefined);
1734: if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE;
1736: computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling);
1738: /* activate all connected components if the netflux has been requested */
1739: if (pcbddc->compute_nonetflux) {
1740: pcbddc->use_vertices = PETSC_TRUE;
1741: pcbddc->use_edges = PETSC_TRUE;
1742: pcbddc->use_faces = PETSC_TRUE;
1743: }
1745: /* Get stdout for dbg */
1746: if (pcbddc->dbg_flag) {
1747: if (!pcbddc->dbg_viewer) pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc));
1748: PetscCall(PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer));
1749: PetscCall(PetscViewerASCIIAddTab(pcbddc->dbg_viewer, 2 * pcbddc->current_level));
1750: }
1752: /* process topology information */
1753: PetscCall(PetscLogEventBegin(PC_BDDC_Topology[pcbddc->current_level], pc, 0, 0, 0));
1754: if (pcbddc->recompute_topography) {
1755: PetscCall(PCBDDCComputeLocalTopologyInfo(pc));
1756: if (pcbddc->discretegradient) PetscCall(PCBDDCNedelecSupport(pc));
1757: }
1758: if (pcbddc->corner_selected) pcbddc->use_vertices = PETSC_TRUE;
1760: /* change basis if requested by the user */
1761: if (pcbddc->user_ChangeOfBasisMatrix) {
1762: /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */
1763: pcbddc->use_change_of_basis = PETSC_FALSE;
1764: PetscCall(PCBDDCComputeLocalMatrix(pc, pcbddc->user_ChangeOfBasisMatrix));
1765: } else {
1766: PetscCall(MatDestroy(&pcbddc->local_mat));
1767: PetscCall(PetscObjectReference((PetscObject)matis->A));
1768: pcbddc->local_mat = matis->A;
1769: }
1771: /*
1772: Compute change of basis on local pressures (aka zerodiag dofs) with the benign trick
1773: This should come earlier than PCISSetUp for extracting the correct subdomain matrices
1774: */
1775: PetscCall(PCBDDCBenignShellMat(pc, PETSC_TRUE));
1776: if (pcbddc->benign_saddle_point) {
1777: PC_IS *pcis = (PC_IS *)pc->data;
1779: if (pcbddc->user_ChangeOfBasisMatrix || pcbddc->use_change_of_basis || !computesubschurs) pcbddc->benign_change_explicit = PETSC_TRUE;
1780: /* detect local saddle point and change the basis in pcbddc->local_mat */
1781: PetscCall(PCBDDCBenignDetectSaddlePoint(pc, (PetscBool)(!pcbddc->recompute_topography), &zerodiag));
1782: /* pop B0 mat from local mat */
1783: PetscCall(PCBDDCBenignPopOrPushB0(pc, PETSC_TRUE));
1784: /* give pcis a hint to not reuse submatrices during PCISCreate */
1785: if (pc->flag == SAME_NONZERO_PATTERN && pcis->reusesubmatrices == PETSC_TRUE) {
1786: if (pcbddc->benign_n && (pcbddc->benign_change_explicit || pcbddc->dbg_flag)) {
1787: pcis->reusesubmatrices = PETSC_FALSE;
1788: } else {
1789: pcis->reusesubmatrices = PETSC_TRUE;
1790: }
1791: } else {
1792: pcis->reusesubmatrices = PETSC_FALSE;
1793: }
1794: }
1796: /* propagate relevant information */
1797: PetscCall(MatIsSymmetricKnown(matis->A, &isset, &issym));
1798: if (isset) PetscCall(MatSetOption(pcbddc->local_mat, MAT_SYMMETRIC, issym));
1799: PetscCall(MatIsSPDKnown(matis->A, &isset, &isspd));
1800: if (isset) PetscCall(MatSetOption(pcbddc->local_mat, MAT_SPD, isspd));
1802: /* Set up all the "iterative substructuring" common block without computing solvers */
1803: {
1804: Mat temp_mat;
1806: temp_mat = matis->A;
1807: matis->A = pcbddc->local_mat;
1808: PetscCall(PCISSetUp(pc, PETSC_TRUE, PETSC_FALSE));
1809: pcbddc->local_mat = matis->A;
1810: matis->A = temp_mat;
1811: }
1813: /* Analyze interface */
1814: if (!pcbddc->graphanalyzed) {
1815: PetscCall(PCBDDCAnalyzeInterface(pc));
1816: computeconstraintsmatrix = PETSC_TRUE;
1817: 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");
1818: PetscCall(MatNullSpaceDestroy(&pcbddc->nonetflux));
1819: if (pcbddc->compute_nonetflux) {
1820: PetscCheck(pcbddc->divudotp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Missing divudotp operator");
1821: PetscCall(PCBDDCComputeNoNetFlux(pc->pmat, pcbddc->divudotp, pcbddc->divudotp_trans, pcbddc->divudotp_vl2l, pcbddc->mat_graph, &pcbddc->nonetflux));
1822: }
1823: }
1824: PetscCall(PetscLogEventEnd(PC_BDDC_Topology[pcbddc->current_level], pc, 0, 0, 0));
1826: /* check existence of a divergence free extension, i.e.
1827: b(v_I,p_0) = 0 for all v_I (raise error if not).
1828: Also, check that PCBDDCBenignGetOrSetP0 works */
1829: if (pcbddc->benign_saddle_point && pcbddc->dbg_flag > 1) PetscCall(PCBDDCBenignCheck(pc, zerodiag));
1830: PetscCall(ISDestroy(&zerodiag));
1832: /* Setup local dirichlet solver ksp_D and sub_schurs solvers */
1833: if (computesubschurs && pcbddc->recompute_topography) PetscCall(PCBDDCInitSubSchurs(pc));
1834: /* SetUp Scaling operator (scaling matrices could be needed in SubSchursSetUp)*/
1835: if (!pcbddc->use_deluxe_scaling) PetscCall(PCBDDCScalingSetUp(pc));
1837: /* finish setup solvers and do adaptive selection of constraints */
1838: sub_schurs = pcbddc->sub_schurs;
1839: if (sub_schurs && sub_schurs->schur_explicit) {
1840: if (computesubschurs) PetscCall(PCBDDCSetUpSubSchurs(pc));
1841: PetscCall(PCBDDCSetUpLocalSolvers(pc, PETSC_TRUE, PETSC_FALSE));
1842: } else {
1843: PetscCall(PCBDDCSetUpLocalSolvers(pc, PETSC_TRUE, PETSC_FALSE));
1844: if (computesubschurs) PetscCall(PCBDDCSetUpSubSchurs(pc));
1845: }
1846: if (pcbddc->adaptive_selection) {
1847: PetscCall(PCBDDCAdaptiveSelection(pc));
1848: computeconstraintsmatrix = PETSC_TRUE;
1849: }
1851: /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */
1852: new_nearnullspace_provided = PETSC_FALSE;
1853: PetscCall(MatGetNearNullSpace(pc->pmat, &nearnullspace));
1854: if (pcbddc->onearnullspace) { /* already used nearnullspace */
1855: if (!nearnullspace) { /* near null space attached to mat has been destroyed */
1856: new_nearnullspace_provided = PETSC_TRUE;
1857: } else {
1858: /* determine if the two nullspaces are different (should be lightweight) */
1859: if (nearnullspace != pcbddc->onearnullspace) {
1860: new_nearnullspace_provided = PETSC_TRUE;
1861: } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */
1862: const Vec *nearnullvecs;
1863: PetscObjectState state;
1864: PetscInt nnsp_size;
1865: PetscCall(MatNullSpaceGetVecs(nearnullspace, NULL, &nnsp_size, &nearnullvecs));
1866: for (PetscInt i = 0; i < nnsp_size; i++) {
1867: PetscCall(PetscObjectStateGet((PetscObject)nearnullvecs[i], &state));
1868: if (pcbddc->onearnullvecs_state[i] != state) {
1869: new_nearnullspace_provided = PETSC_TRUE;
1870: break;
1871: }
1872: }
1873: }
1874: }
1875: } else {
1876: if (!nearnullspace) { /* both nearnullspaces are null */
1877: new_nearnullspace_provided = PETSC_FALSE;
1878: } else { /* nearnullspace attached later */
1879: new_nearnullspace_provided = PETSC_TRUE;
1880: }
1881: }
1883: /* Setup constraints and related work vectors */
1884: /* reset primal space flags */
1885: PetscCall(PetscLogEventBegin(PC_BDDC_LocalWork[pcbddc->current_level], pc, 0, 0, 0));
1886: pcbddc->new_primal_space = PETSC_FALSE;
1887: pcbddc->new_primal_space_local = PETSC_FALSE;
1888: if (computeconstraintsmatrix || new_nearnullspace_provided) {
1889: /* It also sets the primal space flags */
1890: PetscCall(PCBDDCConstraintsSetUp(pc));
1891: }
1892: /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */
1893: PetscCall(PCBDDCSetUpLocalWorkVectors(pc));
1895: if (pcbddc->use_change_of_basis) {
1896: PC_IS *pcis = (PC_IS *)pc->data;
1898: PetscCall(PCBDDCComputeLocalMatrix(pc, pcbddc->ChangeOfBasisMatrix));
1899: if (pcbddc->benign_change) {
1900: PetscCall(MatDestroy(&pcbddc->benign_B0));
1901: /* pop B0 from pcbddc->local_mat */
1902: PetscCall(PCBDDCBenignPopOrPushB0(pc, PETSC_TRUE));
1903: }
1904: /* get submatrices */
1905: PetscCall(MatDestroy(&pcis->A_IB));
1906: PetscCall(MatDestroy(&pcis->A_BI));
1907: PetscCall(MatDestroy(&pcis->A_BB));
1908: PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcis->is_B_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &pcis->A_BB));
1909: PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcis->is_I_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &pcis->A_IB));
1910: PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcis->is_B_local, pcis->is_I_local, MAT_INITIAL_MATRIX, &pcis->A_BI));
1911: /* set flag in pcis to not reuse submatrices during PCISCreate */
1912: pcis->reusesubmatrices = PETSC_FALSE;
1913: } else if (!pcbddc->user_ChangeOfBasisMatrix && !pcbddc->benign_change) {
1914: PetscCall(MatDestroy(&pcbddc->local_mat));
1915: PetscCall(PetscObjectReference((PetscObject)matis->A));
1916: pcbddc->local_mat = matis->A;
1917: }
1919: /* interface pressure block row for B_C */
1920: PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_lP", (PetscObject *)&lP));
1921: PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_lA", (PetscObject *)&lA));
1922: if (lA && lP) {
1923: PC_IS *pcis = (PC_IS *)pc->data;
1924: Mat B_BI, B_BB, Bt_BI, Bt_BB;
1925: PetscBool issym;
1927: PetscCall(MatIsSymmetric(lA, PETSC_SMALL, &issym));
1928: if (issym) {
1929: PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_I_local, MAT_INITIAL_MATRIX, &B_BI));
1930: PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_B_local, MAT_INITIAL_MATRIX, &B_BB));
1931: PetscCall(MatCreateTranspose(B_BI, &Bt_BI));
1932: PetscCall(MatCreateTranspose(B_BB, &Bt_BB));
1933: } else {
1934: PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_I_local, MAT_INITIAL_MATRIX, &B_BI));
1935: PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_B_local, MAT_INITIAL_MATRIX, &B_BB));
1936: PetscCall(MatCreateSubMatrix(lA, pcis->is_I_local, lP, MAT_INITIAL_MATRIX, &Bt_BI));
1937: PetscCall(MatCreateSubMatrix(lA, pcis->is_B_local, lP, MAT_INITIAL_MATRIX, &Bt_BB));
1938: }
1939: PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_B_BI", (PetscObject)B_BI));
1940: PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_B_BB", (PetscObject)B_BB));
1941: PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_Bt_BI", (PetscObject)Bt_BI));
1942: PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_Bt_BB", (PetscObject)Bt_BB));
1943: PetscCall(MatDestroy(&B_BI));
1944: PetscCall(MatDestroy(&B_BB));
1945: PetscCall(MatDestroy(&Bt_BI));
1946: PetscCall(MatDestroy(&Bt_BB));
1947: }
1948: PetscCall(PetscLogEventEnd(PC_BDDC_LocalWork[pcbddc->current_level], pc, 0, 0, 0));
1950: /* SetUp coarse and local Neumann solvers */
1951: PetscCall(PCBDDCSetUpSolvers(pc));
1952: /* SetUp Scaling operator */
1953: if (pcbddc->use_deluxe_scaling) PetscCall(PCBDDCScalingSetUp(pc));
1955: /* mark topography as done */
1956: pcbddc->recompute_topography = PETSC_FALSE;
1958: /* wrap pcis->A_IB and pcis->A_BI if we did not change explicitly the variables on the pressures */
1959: PetscCall(PCBDDCBenignShellMat(pc, PETSC_FALSE));
1961: if (pcbddc->dbg_flag) {
1962: PetscCall(PetscViewerASCIISubtractTab(pcbddc->dbg_viewer, 2 * pcbddc->current_level));
1963: PetscCall(PetscViewerASCIIPopSynchronized(pcbddc->dbg_viewer));
1964: }
1966: { /* Dump customization */
1967: PetscInt save_version = PETSC_DECIDE;
1968: PetscBool flg;
1969: char save[PETSC_MAX_PATH_LEN] = {'\0'};
1971: PetscCall(PetscOptionsGetString(NULL, ((PetscObject)pc)->prefix, "-pc_bddc_save", save, sizeof(save), &flg));
1972: PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)pc)->prefix, "-pc_bddc_save_version", &save_version, NULL));
1973: if (flg) {
1974: size_t len;
1976: PetscCall(PetscStrlen(save, &len));
1977: PetscCall(PCBDDCSaveCustomization(pc, len ? save : NULL, save_version));
1978: }
1979: }
1980: PetscFunctionReturn(PETSC_SUCCESS);
1981: }
1983: static PetscErrorCode PCApply_BDDC(PC pc, Vec r, Vec z)
1984: {
1985: PC_IS *pcis = (PC_IS *)pc->data;
1986: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
1987: Mat lA = NULL;
1988: PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B;
1989: const PetscScalar one = 1.0;
1990: const PetscScalar m_one = -1.0;
1991: const PetscScalar zero = 0.0;
1992: /* This code is similar to that provided in nn.c for PCNN
1993: NN interface preconditioner changed to BDDC
1994: Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */
1996: PetscFunctionBegin;
1997: PetscCall(PetscCitationsRegister(citation, &cited));
1998: if (pcbddc->switch_static) PetscCall(MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat, &lA));
2000: if (pcbddc->ChangeOfBasisMatrix) {
2001: Vec swap;
2003: PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, r, pcbddc->work_change));
2004: swap = pcbddc->work_change;
2005: pcbddc->work_change = r;
2006: r = swap;
2007: /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */
2008: if (pcbddc->benign_apply_coarse_only && pcbddc->use_exact_dirichlet_trick && pcbddc->change_interior) {
2009: PetscCall(VecCopy(r, pcis->vec1_global));
2010: PetscCall(VecLockReadPush(pcis->vec1_global));
2011: }
2012: }
2013: if (pcbddc->benign_have_null) { /* get p0 from r */
2014: PetscCall(PCBDDCBenignGetOrSetP0(pc, r, PETSC_TRUE));
2015: }
2016: if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_DIRICHLET && !pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
2017: PetscCall(VecCopy(r, z));
2018: /* First Dirichlet solve */
2019: PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
2020: PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
2021: /*
2022: Assembling right-hand side for BDDC operator
2023: - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE)
2024: - pcis->vec1_B the interface part of the global vector z
2025: */
2026: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2027: if (n_D) {
2028: PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D));
2029: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2030: PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D));
2031: PetscCall(VecScale(pcis->vec2_D, m_one));
2032: if (pcbddc->switch_static) {
2033: PetscCall(VecSet(pcis->vec1_N, 0.));
2034: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2035: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2036: if (!pcbddc->switch_static_change) PetscCall(MatMult(lA, pcis->vec1_N, pcis->vec2_N));
2037: else {
2038: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2039: PetscCall(MatMult(lA, pcis->vec2_N, pcis->vec1_N));
2040: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2041: }
2042: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
2043: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
2044: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2045: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2046: } else {
2047: PetscCall(MatMult(pcis->A_BI, pcis->vec2_D, pcis->vec1_B));
2048: }
2049: } else {
2050: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2051: PetscCall(VecSet(pcis->vec1_B, zero));
2052: }
2053: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
2054: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
2055: PetscCall(PCBDDCScalingRestriction(pc, z, pcis->vec1_B));
2056: } else {
2057: if (!pcbddc->benign_apply_coarse_only) PetscCall(PCBDDCScalingRestriction(pc, r, pcis->vec1_B));
2058: }
2059: if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_LUMP) {
2060: PetscCheck(pcbddc->switch_static, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "You forgot to pass -pc_bddc_switch_static");
2061: PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
2062: PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
2063: }
2065: /* Apply interface preconditioner
2066: input/output vecs: pcis->vec1_B and pcis->vec1_D */
2067: PetscCall(PCBDDCApplyInterfacePreconditioner(pc, PETSC_FALSE));
2069: /* Apply transpose of partition of unity operator */
2070: PetscCall(PCBDDCScalingExtension(pc, pcis->vec1_B, z));
2071: if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_LUMP) {
2072: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec1_D, z, INSERT_VALUES, SCATTER_REVERSE));
2073: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec1_D, z, INSERT_VALUES, SCATTER_REVERSE));
2074: PetscFunctionReturn(PETSC_SUCCESS);
2075: }
2076: /* Second Dirichlet solve and assembling of output */
2077: PetscCall(VecScatterBegin(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2078: PetscCall(VecScatterEnd(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2079: if (n_B) {
2080: if (pcbddc->switch_static) {
2081: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2082: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2083: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2084: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2085: if (!pcbddc->switch_static_change) PetscCall(MatMult(lA, pcis->vec1_N, pcis->vec2_N));
2086: else {
2087: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2088: PetscCall(MatMult(lA, pcis->vec2_N, pcis->vec1_N));
2089: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2090: }
2091: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
2092: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
2093: } else {
2094: PetscCall(MatMult(pcis->A_IB, pcis->vec1_B, pcis->vec3_D));
2095: }
2096: } else if (pcbddc->switch_static) { /* n_B is zero */
2097: if (!pcbddc->switch_static_change) PetscCall(MatMult(lA, pcis->vec1_D, pcis->vec3_D));
2098: else {
2099: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_D, pcis->vec1_N));
2100: PetscCall(MatMult(lA, pcis->vec1_N, pcis->vec2_N));
2101: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec2_N, pcis->vec3_D));
2102: }
2103: }
2104: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2105: PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec3_D, pcis->vec4_D));
2106: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2107: PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec4_D));
2109: if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
2110: if (pcbddc->switch_static) PetscCall(VecAXPBYPCZ(pcis->vec2_D, m_one, one, m_one, pcis->vec4_D, pcis->vec1_D));
2111: else PetscCall(VecAXPBY(pcis->vec2_D, m_one, m_one, pcis->vec4_D));
2112: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
2113: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
2114: } else {
2115: if (pcbddc->switch_static) PetscCall(VecAXPBY(pcis->vec4_D, one, m_one, pcis->vec1_D));
2116: else PetscCall(VecScale(pcis->vec4_D, m_one));
2117: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
2118: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
2119: }
2120: if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */
2121: if (pcbddc->benign_apply_coarse_only) PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
2122: PetscCall(PCBDDCBenignGetOrSetP0(pc, z, PETSC_FALSE));
2123: }
2125: if (pcbddc->ChangeOfBasisMatrix) {
2126: pcbddc->work_change = r;
2127: PetscCall(VecCopy(z, pcbddc->work_change));
2128: PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcbddc->work_change, z));
2129: }
2130: PetscFunctionReturn(PETSC_SUCCESS);
2131: }
2133: static PetscErrorCode PCApplyTranspose_BDDC(PC pc, Vec r, Vec z)
2134: {
2135: PC_IS *pcis = (PC_IS *)pc->data;
2136: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
2137: Mat lA = NULL;
2138: PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B;
2139: const PetscScalar one = 1.0;
2140: const PetscScalar m_one = -1.0;
2141: const PetscScalar zero = 0.0;
2143: PetscFunctionBegin;
2144: PetscCall(PetscCitationsRegister(citation, &cited));
2145: if (pcbddc->switch_static) PetscCall(MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat, &lA));
2146: if (pcbddc->ChangeOfBasisMatrix) {
2147: Vec swap;
2149: PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, r, pcbddc->work_change));
2150: swap = pcbddc->work_change;
2151: pcbddc->work_change = r;
2152: r = swap;
2153: /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */
2154: if (pcbddc->benign_apply_coarse_only && pcbddc->exact_dirichlet_trick_app && pcbddc->change_interior) {
2155: PetscCall(VecCopy(r, pcis->vec1_global));
2156: PetscCall(VecLockReadPush(pcis->vec1_global));
2157: }
2158: }
2159: if (pcbddc->benign_have_null) { /* get p0 from r */
2160: PetscCall(PCBDDCBenignGetOrSetP0(pc, r, PETSC_TRUE));
2161: }
2162: if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
2163: PetscCall(VecCopy(r, z));
2164: /* First Dirichlet solve */
2165: PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
2166: PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
2167: /*
2168: Assembling right-hand side for BDDC operator
2169: - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE)
2170: - pcis->vec1_B the interface part of the global vector z
2171: */
2172: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2173: if (n_D) {
2174: PetscCall(KSPSolveTranspose(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D));
2175: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2176: PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D));
2177: PetscCall(VecScale(pcis->vec2_D, m_one));
2178: if (pcbddc->switch_static) {
2179: PetscCall(VecSet(pcis->vec1_N, 0.));
2180: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2181: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2182: if (!pcbddc->switch_static_change) {
2183: PetscCall(MatMultTranspose(lA, pcis->vec1_N, pcis->vec2_N));
2184: } else {
2185: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2186: PetscCall(MatMultTranspose(lA, pcis->vec2_N, pcis->vec1_N));
2187: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2188: }
2189: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
2190: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
2191: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2192: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2193: } else {
2194: PetscCall(MatMultTranspose(pcis->A_IB, pcis->vec2_D, pcis->vec1_B));
2195: }
2196: } else {
2197: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2198: PetscCall(VecSet(pcis->vec1_B, zero));
2199: }
2200: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
2201: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
2202: PetscCall(PCBDDCScalingRestriction(pc, z, pcis->vec1_B));
2203: } else {
2204: PetscCall(PCBDDCScalingRestriction(pc, r, pcis->vec1_B));
2205: }
2207: /* Apply interface preconditioner
2208: input/output vecs: pcis->vec1_B and pcis->vec1_D */
2209: PetscCall(PCBDDCApplyInterfacePreconditioner(pc, PETSC_TRUE));
2211: /* Apply transpose of partition of unity operator */
2212: PetscCall(PCBDDCScalingExtension(pc, pcis->vec1_B, z));
2214: /* Second Dirichlet solve and assembling of output */
2215: PetscCall(VecScatterBegin(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2216: PetscCall(VecScatterEnd(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2217: if (n_B) {
2218: if (pcbddc->switch_static) {
2219: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2220: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2221: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2222: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2223: if (!pcbddc->switch_static_change) {
2224: PetscCall(MatMultTranspose(lA, pcis->vec1_N, pcis->vec2_N));
2225: } else {
2226: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2227: PetscCall(MatMultTranspose(lA, pcis->vec2_N, pcis->vec1_N));
2228: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2229: }
2230: PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
2231: PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
2232: } else {
2233: PetscCall(MatMultTranspose(pcis->A_BI, pcis->vec1_B, pcis->vec3_D));
2234: }
2235: } else if (pcbddc->switch_static) { /* n_B is zero */
2236: if (!pcbddc->switch_static_change) {
2237: PetscCall(MatMultTranspose(lA, pcis->vec1_D, pcis->vec3_D));
2238: } else {
2239: PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_D, pcis->vec1_N));
2240: PetscCall(MatMultTranspose(lA, pcis->vec1_N, pcis->vec2_N));
2241: PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec2_N, pcis->vec3_D));
2242: }
2243: }
2244: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2245: PetscCall(KSPSolveTranspose(pcbddc->ksp_D, pcis->vec3_D, pcis->vec4_D));
2246: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2247: PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec4_D));
2248: if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
2249: if (pcbddc->switch_static) PetscCall(VecAXPBYPCZ(pcis->vec2_D, m_one, one, m_one, pcis->vec4_D, pcis->vec1_D));
2250: else PetscCall(VecAXPBY(pcis->vec2_D, m_one, m_one, pcis->vec4_D));
2251: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
2252: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
2253: } else {
2254: if (pcbddc->switch_static) PetscCall(VecAXPBY(pcis->vec4_D, one, m_one, pcis->vec1_D));
2255: else PetscCall(VecScale(pcis->vec4_D, m_one));
2256: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
2257: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
2258: }
2259: if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */
2260: PetscCall(PCBDDCBenignGetOrSetP0(pc, z, PETSC_FALSE));
2261: }
2262: if (pcbddc->ChangeOfBasisMatrix) {
2263: pcbddc->work_change = r;
2264: PetscCall(VecCopy(z, pcbddc->work_change));
2265: PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcbddc->work_change, z));
2266: }
2267: PetscFunctionReturn(PETSC_SUCCESS);
2268: }
2270: static PetscErrorCode PCReset_BDDC(PC pc)
2271: {
2272: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
2273: PC_IS *pcis = (PC_IS *)pc->data;
2274: PetscScalar scaling_factor = pcis->scaling_factor;
2276: PetscFunctionBegin;
2277: /* free BDDC custom data */
2278: PetscCall(PCBDDCResetCustomization(pc));
2279: /* destroy objects related to topography */
2280: PetscCall(PCBDDCResetTopography(pc));
2281: /* destroy objects for scaling operator */
2282: PetscCall(PCBDDCScalingDestroy(pc));
2283: /* free solvers stuff */
2284: PetscCall(PCBDDCResetSolvers(pc));
2285: /* free global vectors needed in presolve */
2286: PetscCall(VecDestroy(&pcbddc->temp_solution));
2287: PetscCall(VecDestroy(&pcbddc->original_rhs));
2288: /* free data created by PCIS */
2289: PetscCall(PCISReset(pc));
2291: // Restore the PCIS callbacks and clear derived state while preserving configured options.
2292: PetscCall(PCISInitialize(pc));
2293: pcis->scaling_factor = scaling_factor;
2294: pcis->n = 0;
2295: pcis->n_B = 0;
2296: pcis->pure_neumann = PETSC_FALSE;
2297: pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
2298: pcbddc->ksp_guess_nonzero = PETSC_FALSE;
2299: pcbddc->rhs_change = PETSC_FALSE;
2300: pcbddc->temp_solution_used = PETSC_FALSE;
2301: PetscFunctionReturn(PETSC_SUCCESS);
2302: }
2304: static PetscErrorCode PCDestroy_BDDC(PC pc)
2305: {
2306: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
2308: PetscFunctionBegin;
2309: PetscCall(PCReset_BDDC(pc));
2310: PetscCall(PCBDDCGraphDestroy(&pcbddc->mat_graph));
2311: PetscCall(PCISReset(pc));
2312: PetscCall(KSPDestroy(&pcbddc->ksp_D));
2313: PetscCall(KSPDestroy(&pcbddc->ksp_R));
2314: PetscCall(KSPDestroy(&pcbddc->coarse_ksp));
2315: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDiscreteGradient_C", NULL));
2316: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDivergenceMat_C", NULL));
2317: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetChangeOfBasisMat_C", NULL));
2318: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesLocalIS_C", NULL));
2319: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesIS_C", NULL));
2320: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesLocalIS_C", NULL));
2321: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesIS_C", NULL));
2322: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetCoarseningRatio_C", NULL));
2323: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevel_C", NULL));
2324: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetUseExactDirichlet_C", NULL));
2325: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevels_C", NULL));
2326: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCLoadCustomization_C", NULL));
2327: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSaveCustomization_C", NULL));
2328: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundaries_C", NULL));
2329: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundariesLocal_C", NULL));
2330: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundaries_C", NULL));
2331: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundariesLocal_C", NULL));
2332: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundaries_C", NULL));
2333: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundariesLocal_C", NULL));
2334: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundaries_C", NULL));
2335: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundariesLocal_C", NULL));
2336: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplitting_C", NULL));
2337: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplittingLocal_C", NULL));
2338: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", NULL));
2339: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCCreateFETIDPOperators_C", NULL));
2340: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetRHS_C", NULL));
2341: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetSolution_C", NULL));
2342: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", NULL));
2343: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
2344: PetscCall(PetscFree(pc->data));
2345: PetscFunctionReturn(PETSC_SUCCESS);
2346: }
2348: static PetscErrorCode PCSetCoordinates_BDDC(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
2349: {
2350: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
2351: PCBDDCGraph mat_graph = pcbddc->mat_graph;
2353: PetscFunctionBegin;
2354: PetscCall(PetscFree(mat_graph->coords));
2355: PetscCall(PetscMalloc1(nloc * dim, &mat_graph->coords));
2356: PetscCall(PetscArraycpy(mat_graph->coords, coords, nloc * dim));
2357: mat_graph->cnloc = nloc;
2358: mat_graph->cdim = dim;
2359: mat_graph->cloc = PETSC_FALSE;
2360: /* flg setup */
2361: pcbddc->recompute_topography = PETSC_TRUE;
2362: pcbddc->corner_selected = PETSC_FALSE;
2363: PetscFunctionReturn(PETSC_SUCCESS);
2364: }
2366: static PetscErrorCode PCPreSolveChangeRHS_BDDC(PC pc, PetscBool *change)
2367: {
2368: PetscFunctionBegin;
2369: *change = PETSC_TRUE;
2370: PetscFunctionReturn(PETSC_SUCCESS);
2371: }
2373: static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
2374: {
2375: FETIDPMat_ctx mat_ctx;
2376: Vec work;
2377: PC_IS *pcis;
2378: PC_BDDC *pcbddc;
2380: PetscFunctionBegin;
2381: PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2382: pcis = (PC_IS *)mat_ctx->pc->data;
2383: pcbddc = (PC_BDDC *)mat_ctx->pc->data;
2385: PetscCall(VecSet(fetidp_flux_rhs, 0.0));
2386: /* copy rhs since we may change it during PCPreSolve_BDDC */
2387: if (!pcbddc->original_rhs) PetscCall(VecDuplicate(pcis->vec1_global, &pcbddc->original_rhs));
2388: if (mat_ctx->rhs_flip) {
2389: PetscCall(VecPointwiseMult(pcbddc->original_rhs, standard_rhs, mat_ctx->rhs_flip));
2390: } else {
2391: PetscCall(VecCopy(standard_rhs, pcbddc->original_rhs));
2392: }
2393: if (mat_ctx->g2g_p) {
2394: /* interface pressure rhs */
2395: PetscCall(VecScatterBegin(mat_ctx->g2g_p, fetidp_flux_rhs, pcbddc->original_rhs, INSERT_VALUES, SCATTER_REVERSE));
2396: PetscCall(VecScatterEnd(mat_ctx->g2g_p, fetidp_flux_rhs, pcbddc->original_rhs, INSERT_VALUES, SCATTER_REVERSE));
2397: PetscCall(VecScatterBegin(mat_ctx->g2g_p, standard_rhs, fetidp_flux_rhs, INSERT_VALUES, SCATTER_FORWARD));
2398: PetscCall(VecScatterEnd(mat_ctx->g2g_p, standard_rhs, fetidp_flux_rhs, INSERT_VALUES, SCATTER_FORWARD));
2399: if (!mat_ctx->rhs_flip) PetscCall(VecScale(fetidp_flux_rhs, -1.));
2400: }
2401: /*
2402: change of basis for physical rhs if needed
2403: It also changes the rhs in case of dirichlet boundaries
2404: */
2405: PetscCall(PCPreSolve_BDDC(mat_ctx->pc, NULL, pcbddc->original_rhs, NULL));
2406: if (pcbddc->ChangeOfBasisMatrix) {
2407: PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, pcbddc->original_rhs, pcbddc->work_change));
2408: work = pcbddc->work_change;
2409: } else {
2410: work = pcbddc->original_rhs;
2411: }
2412: /* store vectors for computation of fetidp final solution */
2413: PetscCall(VecScatterBegin(pcis->global_to_D, work, mat_ctx->temp_solution_D, INSERT_VALUES, SCATTER_FORWARD));
2414: PetscCall(VecScatterEnd(pcis->global_to_D, work, mat_ctx->temp_solution_D, INSERT_VALUES, SCATTER_FORWARD));
2415: /* scale rhs since it should be unassembled */
2416: /* TODO use counter scaling? (also below) */
2417: PetscCall(VecScatterBegin(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2418: PetscCall(VecScatterEnd(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2419: /* Apply partition of unity */
2420: PetscCall(VecPointwiseMult(mat_ctx->temp_solution_B, pcis->D, mat_ctx->temp_solution_B));
2421: /* PetscCall(PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B)); */
2422: if (!pcbddc->switch_static) {
2423: /* compute partially subassembled Schur complement right-hand side */
2424: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2425: PetscCall(KSPSolve(pcbddc->ksp_D, mat_ctx->temp_solution_D, pcis->vec1_D));
2426: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2427: /* Cannot propagate up error in KSPSolve() because there is no access to the PC */
2428: PetscCall(MatMult(pcis->A_BI, pcis->vec1_D, pcis->vec1_B));
2429: PetscCall(VecAXPY(mat_ctx->temp_solution_B, -1.0, pcis->vec1_B));
2430: PetscCall(VecSet(work, 0.0));
2431: PetscCall(VecScatterBegin(pcis->global_to_B, mat_ctx->temp_solution_B, work, ADD_VALUES, SCATTER_REVERSE));
2432: PetscCall(VecScatterEnd(pcis->global_to_B, mat_ctx->temp_solution_B, work, ADD_VALUES, SCATTER_REVERSE));
2433: /* PetscCall(PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B)); */
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: PetscCall(VecPointwiseMult(mat_ctx->temp_solution_B, pcis->D, mat_ctx->temp_solution_B));
2437: }
2438: /* BDDC rhs */
2439: PetscCall(VecCopy(mat_ctx->temp_solution_B, pcis->vec1_B));
2440: if (pcbddc->switch_static) PetscCall(VecCopy(mat_ctx->temp_solution_D, pcis->vec1_D));
2441: /* apply BDDC */
2442: PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
2443: PetscCall(PCBDDCApplyInterfacePreconditioner(mat_ctx->pc, PETSC_FALSE));
2444: PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
2446: /* Application of B_delta and assembling of rhs for fetidp fluxes */
2447: PetscCall(MatMult(mat_ctx->B_delta, pcis->vec1_B, mat_ctx->lambda_local));
2448: PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, mat_ctx->lambda_local, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2449: PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, mat_ctx->lambda_local, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2450: /* Add contribution to interface pressures */
2451: if (mat_ctx->l2g_p) {
2452: PetscCall(VecISSet(pcis->vec1_B, mat_ctx->lP_B, 0));
2453: PetscCall(MatMult(mat_ctx->B_BB, pcis->vec1_B, mat_ctx->vP));
2454: if (pcbddc->switch_static) {
2455: PetscCall(VecISSet(pcis->vec1_D, mat_ctx->lP_I, 0));
2456: PetscCall(MatMultAdd(mat_ctx->B_BI, pcis->vec1_D, mat_ctx->vP, mat_ctx->vP));
2457: }
2458: PetscCall(VecScatterBegin(mat_ctx->l2g_p, mat_ctx->vP, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2459: PetscCall(VecScatterEnd(mat_ctx->l2g_p, mat_ctx->vP, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2460: }
2461: PetscFunctionReturn(PETSC_SUCCESS);
2462: }
2464: /*@
2465: PCBDDCMatFETIDPGetRHS - Computes the right-hand side of a FETI-DP system from the physical right-hand side
2467: Collective
2469: Input Parameters:
2470: + fetidp_mat - the FETI-DP matrix obtained from `PCBDDCCreateFETIDPOperators()`
2471: - standard_rhs - the right-hand side of the original linear system
2473: Output Parameter:
2474: . fetidp_flux_rhs - vector in which to store the right-hand side of the FETI-DP system
2476: Level: developer
2478: Note:
2479: Most users should employ the `KSP` interface for linear solvers and create a solver of type `KSPFETIDP`.
2481: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCCreateFETIDPOperators()`, `PCBDDCMatFETIDPGetSolution()`
2482: @*/
2483: PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
2484: {
2485: FETIDPMat_ctx mat_ctx;
2487: PetscFunctionBegin;
2491: PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2492: PetscUseMethod(mat_ctx->pc, "PCBDDCMatFETIDPGetRHS_C", (Mat, Vec, Vec), (fetidp_mat, standard_rhs, fetidp_flux_rhs));
2493: PetscFunctionReturn(PETSC_SUCCESS);
2494: }
2496: static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
2497: {
2498: FETIDPMat_ctx mat_ctx;
2499: PC_IS *pcis;
2500: PC_BDDC *pcbddc;
2501: Vec work;
2503: PetscFunctionBegin;
2504: PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2505: pcis = (PC_IS *)mat_ctx->pc->data;
2506: pcbddc = (PC_BDDC *)mat_ctx->pc->data;
2508: /* apply B_delta^T */
2509: PetscCall(VecSet(pcis->vec1_B, 0.));
2510: PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, fetidp_flux_sol, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
2511: PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, fetidp_flux_sol, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
2512: PetscCall(MatMultTranspose(mat_ctx->B_delta, mat_ctx->lambda_local, pcis->vec1_B));
2513: if (mat_ctx->l2g_p) {
2514: PetscCall(VecScatterBegin(mat_ctx->l2g_p, fetidp_flux_sol, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
2515: PetscCall(VecScatterEnd(mat_ctx->l2g_p, fetidp_flux_sol, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
2516: PetscCall(MatMultAdd(mat_ctx->Bt_BB, mat_ctx->vP, pcis->vec1_B, pcis->vec1_B));
2517: }
2519: /* compute rhs for BDDC application */
2520: PetscCall(VecAYPX(pcis->vec1_B, -1.0, mat_ctx->temp_solution_B));
2521: if (pcbddc->switch_static) {
2522: PetscCall(VecCopy(mat_ctx->temp_solution_D, pcis->vec1_D));
2523: if (mat_ctx->l2g_p) {
2524: PetscCall(VecScale(mat_ctx->vP, -1.));
2525: PetscCall(MatMultAdd(mat_ctx->Bt_BI, mat_ctx->vP, pcis->vec1_D, pcis->vec1_D));
2526: }
2527: }
2529: /* apply BDDC */
2530: PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
2531: PetscCall(PCBDDCApplyInterfacePreconditioner(mat_ctx->pc, PETSC_FALSE));
2533: /* put values into global vector */
2534: if (pcbddc->ChangeOfBasisMatrix) work = pcbddc->work_change;
2535: else work = standard_sol;
2536: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, work, INSERT_VALUES, SCATTER_REVERSE));
2537: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, work, INSERT_VALUES, SCATTER_REVERSE));
2538: if (!pcbddc->switch_static) {
2539: /* compute values into the interior if solved for the partially subassembled Schur complement */
2540: PetscCall(MatMult(pcis->A_IB, pcis->vec1_B, pcis->vec1_D));
2541: PetscCall(VecAYPX(pcis->vec1_D, -1.0, mat_ctx->temp_solution_D));
2542: PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2543: PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec1_D));
2544: PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2545: /* Cannot propagate up error in KSPSolve() because there is no access to the PC */
2546: }
2548: PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec1_D, work, INSERT_VALUES, SCATTER_REVERSE));
2549: PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec1_D, work, INSERT_VALUES, SCATTER_REVERSE));
2550: /* add p0 solution to final solution */
2551: PetscCall(PCBDDCBenignGetOrSetP0(mat_ctx->pc, work, PETSC_FALSE));
2552: if (pcbddc->ChangeOfBasisMatrix) PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, work, standard_sol));
2553: PetscCall(PCPostSolve_BDDC(mat_ctx->pc, NULL, NULL, standard_sol));
2554: if (mat_ctx->g2g_p) {
2555: PetscCall(VecScatterBegin(mat_ctx->g2g_p, fetidp_flux_sol, standard_sol, INSERT_VALUES, SCATTER_REVERSE));
2556: PetscCall(VecScatterEnd(mat_ctx->g2g_p, fetidp_flux_sol, standard_sol, INSERT_VALUES, SCATTER_REVERSE));
2557: }
2558: PetscFunctionReturn(PETSC_SUCCESS);
2559: }
2561: static PetscErrorCode PCView_BDDCIPC(PC pc, PetscViewer viewer)
2562: {
2563: BDDCIPC_ctx bddcipc_ctx;
2564: PetscBool isascii;
2566: PetscFunctionBegin;
2567: PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2568: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2569: if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "BDDC interface preconditioner\n"));
2570: PetscCall(PetscViewerASCIIPushTab(viewer));
2571: PetscCall(PCView(bddcipc_ctx->bddc, viewer));
2572: PetscCall(PetscViewerASCIIPopTab(viewer));
2573: PetscFunctionReturn(PETSC_SUCCESS);
2574: }
2576: static PetscErrorCode PCSetUp_BDDCIPC(PC pc)
2577: {
2578: BDDCIPC_ctx bddcipc_ctx;
2579: PetscBool isbddc;
2580: Vec vv;
2581: IS is;
2582: PC_IS *pcis;
2584: PetscFunctionBegin;
2585: PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2586: PetscCall(PetscObjectTypeCompare((PetscObject)bddcipc_ctx->bddc, PCBDDC, &isbddc));
2587: PetscCheck(isbddc, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Invalid type %s. Must be of type bddc", ((PetscObject)bddcipc_ctx->bddc)->type_name);
2588: PetscCall(PCSetUp(bddcipc_ctx->bddc));
2590: /* create interface scatter */
2591: pcis = (PC_IS *)bddcipc_ctx->bddc->data;
2592: PetscCall(VecScatterDestroy(&bddcipc_ctx->g2l));
2593: PetscCall(MatCreateVecs(pc->pmat, &vv, NULL));
2594: PetscCall(ISRenumber(pcis->is_B_global, NULL, NULL, &is));
2595: PetscCall(VecScatterCreate(vv, is, pcis->vec1_B, NULL, &bddcipc_ctx->g2l));
2596: PetscCall(ISDestroy(&is));
2597: PetscCall(VecDestroy(&vv));
2598: PetscFunctionReturn(PETSC_SUCCESS);
2599: }
2601: static PetscErrorCode PCApply_BDDCIPC(PC pc, Vec r, Vec x)
2602: {
2603: BDDCIPC_ctx bddcipc_ctx;
2604: PC_IS *pcis;
2605: VecScatter tmps;
2607: PetscFunctionBegin;
2608: PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2609: pcis = (PC_IS *)bddcipc_ctx->bddc->data;
2610: tmps = pcis->global_to_B;
2611: pcis->global_to_B = bddcipc_ctx->g2l;
2612: PetscCall(PCBDDCScalingRestriction(bddcipc_ctx->bddc, r, pcis->vec1_B));
2613: PetscCall(PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc, PETSC_FALSE));
2614: PetscCall(PCBDDCScalingExtension(bddcipc_ctx->bddc, pcis->vec1_B, x));
2615: pcis->global_to_B = tmps;
2616: PetscFunctionReturn(PETSC_SUCCESS);
2617: }
2619: static PetscErrorCode PCApplyTranspose_BDDCIPC(PC pc, Vec r, Vec x)
2620: {
2621: BDDCIPC_ctx bddcipc_ctx;
2622: PC_IS *pcis;
2623: VecScatter tmps;
2625: PetscFunctionBegin;
2626: PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2627: pcis = (PC_IS *)bddcipc_ctx->bddc->data;
2628: tmps = pcis->global_to_B;
2629: pcis->global_to_B = bddcipc_ctx->g2l;
2630: PetscCall(PCBDDCScalingRestriction(bddcipc_ctx->bddc, r, pcis->vec1_B));
2631: PetscCall(PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc, PETSC_TRUE));
2632: PetscCall(PCBDDCScalingExtension(bddcipc_ctx->bddc, pcis->vec1_B, x));
2633: pcis->global_to_B = tmps;
2634: PetscFunctionReturn(PETSC_SUCCESS);
2635: }
2637: static PetscErrorCode PCDestroy_BDDCIPC(PC pc)
2638: {
2639: BDDCIPC_ctx bddcipc_ctx;
2641: PetscFunctionBegin;
2642: PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2643: PetscCall(PCDestroy(&bddcipc_ctx->bddc));
2644: PetscCall(VecScatterDestroy(&bddcipc_ctx->g2l));
2645: PetscCall(PetscFree(bddcipc_ctx));
2646: PetscFunctionReturn(PETSC_SUCCESS);
2647: }
2649: /*@
2650: PCBDDCMatFETIDPGetSolution - Computes the physical solution from the solution of a FETI-DP system
2652: Collective
2654: Input Parameters:
2655: + fetidp_mat - the FETI-DP matrix obtained from `PCBDDCCreateFETIDPOperators()`
2656: - fetidp_flux_sol - the solution of the FETI-DP linear system
2658: Output Parameter:
2659: . standard_sol - vector in which to store the solution on the physical domain
2661: Level: developer
2663: Note:
2664: Most users should employ the `KSP` interface for linear solvers and create a solver of type `KSPFETIDP`.
2666: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCCreateFETIDPOperators()`, `PCBDDCMatFETIDPGetRHS()`
2667: @*/
2668: PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
2669: {
2670: FETIDPMat_ctx mat_ctx;
2672: PetscFunctionBegin;
2676: PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2677: PetscUseMethod(mat_ctx->pc, "PCBDDCMatFETIDPGetSolution_C", (Mat, Vec, Vec), (fetidp_mat, fetidp_flux_sol, standard_sol));
2678: PetscFunctionReturn(PETSC_SUCCESS);
2679: }
2681: static PetscErrorCode MatISSubMatrixEmbedLocalIS(Mat A, IS oldis, IS *newis)
2682: {
2683: Mat_IS *matis = (Mat_IS *)A->data;
2684: ISLocalToGlobalMapping ltog;
2685: IS is;
2687: PetscFunctionBegin;
2688: PetscCheck(matis->getsub_ris, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing getsub IS");
2689: PetscCall(ISLocalToGlobalMappingCreateIS(matis->getsub_ris, <og));
2690: PetscCall(ISGlobalToLocalMappingApplyIS(ltog, IS_GTOLM_DROP, oldis, &is));
2691: PetscCall(ISOnComm(is, PetscObjectComm((PetscObject)A), PETSC_COPY_VALUES, newis));
2692: PetscCall(ISLocalToGlobalMappingDestroy(<og));
2693: PetscCall(ISDestroy(&is));
2694: PetscFunctionReturn(PETSC_SUCCESS);
2695: }
2697: static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, PetscBool fully_redundant, const char *prefix, Mat *fetidp_mat, PC *fetidp_pc)
2698: {
2699: FETIDPMat_ctx fetidpmat_ctx;
2700: Mat newmat;
2701: FETIDPPC_ctx fetidppc_ctx;
2702: PC newpc;
2703: MPI_Comm comm;
2705: PetscFunctionBegin;
2706: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
2707: /* FETI-DP matrix */
2708: PetscCall(PCBDDCCreateFETIDPMatContext(pc, &fetidpmat_ctx));
2709: fetidpmat_ctx->fully_redundant = fully_redundant;
2710: PetscCall(PCBDDCSetupFETIDPMatContext(fetidpmat_ctx));
2711: PetscCall(MatCreateShell(comm, fetidpmat_ctx->n, fetidpmat_ctx->n, fetidpmat_ctx->N, fetidpmat_ctx->N, fetidpmat_ctx, &newmat));
2712: PetscCall(PetscObjectSetName((PetscObject)newmat, !fetidpmat_ctx->l2g_lambda_only ? "F" : "G"));
2713: PetscCall(MatShellSetOperation(newmat, MATOP_MULT, (PetscErrorCodeFn *)FETIDPMatMult));
2714: PetscCall(MatShellSetOperation(newmat, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)FETIDPMatMultTranspose));
2715: PetscCall(MatShellSetOperation(newmat, MATOP_DESTROY, (PetscErrorCodeFn *)PCBDDCDestroyFETIDPMat));
2716: /* propagate MatOptions */
2717: {
2718: PC_BDDC *pcbddc = (PC_BDDC *)fetidpmat_ctx->pc->data;
2719: PetscBool isset, issym;
2721: PetscCall(MatIsSymmetricKnown(pc->mat, &isset, &issym));
2722: if ((isset && issym) || pcbddc->symmetric_primal) PetscCall(MatSetOption(newmat, MAT_SYMMETRIC, PETSC_TRUE));
2723: }
2724: PetscCall(MatSetOptionsPrefix(newmat, prefix));
2725: PetscCall(MatAppendOptionsPrefix(newmat, "fetidp_"));
2726: PetscCall(MatSetUp(newmat));
2727: /* FETI-DP preconditioner */
2728: PetscCall(PCBDDCCreateFETIDPPCContext(pc, &fetidppc_ctx));
2729: PetscCall(PCBDDCSetupFETIDPPCContext(newmat, fetidppc_ctx));
2730: PetscCall(PCCreate(comm, &newpc));
2731: PetscCall(PCSetOperators(newpc, newmat, newmat));
2732: PetscCall(PCSetOptionsPrefix(newpc, prefix));
2733: PetscCall(PCAppendOptionsPrefix(newpc, "fetidp_"));
2734: PetscCall(PCSetErrorIfFailure(newpc, pc->erroriffailure));
2735: if (!fetidpmat_ctx->l2g_lambda_only) { /* standard FETI-DP */
2736: PetscCall(PCSetType(newpc, PCSHELL));
2737: PetscCall(PCShellSetName(newpc, "FETI-DP multipliers"));
2738: PetscCall(PCShellSetContext(newpc, fetidppc_ctx));
2739: PetscCall(PCShellSetApply(newpc, FETIDPPCApply));
2740: PetscCall(PCShellSetApplyTranspose(newpc, FETIDPPCApplyTranspose));
2741: PetscCall(PCShellSetView(newpc, FETIDPPCView));
2742: PetscCall(PCShellSetDestroy(newpc, PCBDDCDestroyFETIDPPC));
2743: } else { /* saddle-point FETI-DP */
2744: Mat M;
2745: PetscInt psize;
2746: PetscBool fake = PETSC_FALSE, isfieldsplit;
2748: PetscCall(ISViewFromOptions(fetidpmat_ctx->lagrange, NULL, "-lag_view"));
2749: PetscCall(ISViewFromOptions(fetidpmat_ctx->pressure, NULL, "-press_view"));
2750: PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_PPmat", (PetscObject *)&M));
2751: PetscCall(PCSetType(newpc, PCFIELDSPLIT));
2752: PetscCall(PCFieldSplitSetIS(newpc, "lag", fetidpmat_ctx->lagrange));
2753: PetscCall(PCFieldSplitSetIS(newpc, "p", fetidpmat_ctx->pressure));
2754: PetscCall(PCFieldSplitSetType(newpc, PC_COMPOSITE_SCHUR));
2755: PetscCall(PCFieldSplitSetSchurFactType(newpc, PC_FIELDSPLIT_SCHUR_FACT_DIAG));
2756: PetscCall(ISGetSize(fetidpmat_ctx->pressure, &psize));
2757: if (psize != M->rmap->N) {
2758: Mat M2;
2759: PetscInt lpsize;
2761: fake = PETSC_TRUE;
2762: PetscCall(ISGetLocalSize(fetidpmat_ctx->pressure, &lpsize));
2763: PetscCall(MatCreate(comm, &M2));
2764: PetscCall(MatSetType(M2, MATAIJ));
2765: PetscCall(MatSetSizes(M2, lpsize, lpsize, psize, psize));
2766: PetscCall(MatSetUp(M2));
2767: PetscCall(MatAssemblyBegin(M2, MAT_FINAL_ASSEMBLY));
2768: PetscCall(MatAssemblyEnd(M2, MAT_FINAL_ASSEMBLY));
2769: PetscCall(PCFieldSplitSetSchurPre(newpc, PC_FIELDSPLIT_SCHUR_PRE_USER, M2));
2770: PetscCall(MatDestroy(&M2));
2771: } else {
2772: PetscCall(PCFieldSplitSetSchurPre(newpc, PC_FIELDSPLIT_SCHUR_PRE_USER, M));
2773: }
2774: PetscCall(PCFieldSplitSetSchurScale(newpc, 1.0));
2776: /* we need to setfromoptions and setup here to access the blocks */
2777: PetscCall(PCSetFromOptions(newpc));
2778: PetscCall(PCSetUp(newpc));
2780: /* user may have changed the type (e.g. -fetidp_pc_type none) */
2781: PetscCall(PetscObjectTypeCompare((PetscObject)newpc, PCFIELDSPLIT, &isfieldsplit));
2782: if (isfieldsplit) {
2783: KSP *ksps;
2784: PC ppc, lagpc;
2785: PetscInt nn;
2786: PetscBool ismatis, matisok = PETSC_FALSE, check = PETSC_FALSE;
2788: /* set the solver for the (0,0) block */
2789: PetscCall(PCFieldSplitSchurGetSubKSP(newpc, &nn, &ksps));
2790: if (!nn) { /* not of type PC_COMPOSITE_SCHUR */
2791: PetscCall(PCFieldSplitGetSubKSP(newpc, &nn, &ksps));
2792: if (!fake) { /* pass pmat to the pressure solver */
2793: Mat F;
2795: PetscCall(KSPGetOperators(ksps[1], &F, NULL));
2796: PetscCall(KSPSetOperators(ksps[1], F, M));
2797: }
2798: } else {
2799: PetscBool issym, isset;
2800: Mat S;
2802: PetscCall(PCFieldSplitSchurGetS(newpc, &S));
2803: PetscCall(MatIsSymmetricKnown(newmat, &isset, &issym));
2804: if (isset) PetscCall(MatSetOption(S, MAT_SYMMETRIC, issym));
2805: }
2806: PetscCall(KSPGetPC(ksps[0], &lagpc));
2807: PetscCall(PCSetType(lagpc, PCSHELL));
2808: PetscCall(PCShellSetName(lagpc, "FETI-DP multipliers"));
2809: PetscCall(PCShellSetContext(lagpc, fetidppc_ctx));
2810: PetscCall(PCShellSetApply(lagpc, FETIDPPCApply));
2811: PetscCall(PCShellSetApplyTranspose(lagpc, FETIDPPCApplyTranspose));
2812: PetscCall(PCShellSetView(lagpc, FETIDPPCView));
2813: PetscCall(PCShellSetDestroy(lagpc, PCBDDCDestroyFETIDPPC));
2815: /* Olof's idea: interface Schur complement preconditioner for the mass matrix */
2816: PetscCall(KSPGetPC(ksps[1], &ppc));
2817: if (fake) {
2818: PC_BDDC *pcbddc = (PC_BDDC *)fetidpmat_ctx->pc->data;
2819: BDDCIPC_ctx bddcipc_ctx;
2820: PetscContainer c;
2822: matisok = PETSC_TRUE;
2824: /* create inner BDDC solver */
2825: PetscCall(PetscNew(&bddcipc_ctx));
2826: PetscCall(PCCreate(comm, &bddcipc_ctx->bddc));
2827: PetscCall(PCSetType(bddcipc_ctx->bddc, PCBDDC));
2828: PetscCall(PCSetOperators(bddcipc_ctx->bddc, M, M));
2829: PetscCall(PetscObjectTypeCompare((PetscObject)M, MATIS, &ismatis));
2830: PetscCheck(ismatis, comm, PETSC_ERR_PLIB, "Matrix type %s not of type MATIS", ((PetscObject)M)->type_name);
2831: /* the inner bddc for FETI-DP is already setup, we have local info available */
2832: if (pcbddc->user_primal_vertices_local || pcbddc->n_ISForDofsLocal > 2) {
2833: if (pcbddc->user_primal_vertices_local) {
2834: IS primals;
2836: PetscCall(MatISSubMatrixEmbedLocalIS(M, pcbddc->user_primal_vertices_local, &primals));
2837: PetscCall(PCBDDCSetPrimalVerticesLocalIS(bddcipc_ctx->bddc, primals));
2838: PetscCall(ISDestroy(&primals));
2839: }
2840: if (pcbddc->n_ISForDofsLocal > 2) { /* no need to propagate info if nfields < 3 */
2841: IS *split;
2842: PetscInt i, nf;
2844: PetscCall(PetscCalloc1(pcbddc->n_ISForDofsLocal, &split));
2845: for (i = 0, nf = 0; i < pcbddc->n_ISForDofsLocal; i++) {
2846: PetscInt ns;
2848: PetscCall(MatISSubMatrixEmbedLocalIS(M, pcbddc->ISForDofsLocal[i], &split[nf]));
2849: PetscCall(ISGetSize(split[nf], &ns));
2850: if (!ns) PetscCall(ISDestroy(&split[nf]));
2851: else nf++;
2852: }
2853: PetscCall(PCBDDCSetDofsSplittingLocal(bddcipc_ctx->bddc, nf, split));
2854: for (i = 0; i < nf; i++) PetscCall(ISDestroy(&split[i]));
2855: PetscCall(PetscFree(split));
2856: }
2857: }
2858: PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_pCSR", (PetscObject *)&c));
2859: PetscCall(PetscObjectTypeCompare((PetscObject)M, MATIS, &ismatis));
2860: if (c && ismatis) {
2861: Mat lM;
2862: PetscInt *csr, n;
2864: PetscCall(MatISGetLocalMat(M, &lM));
2865: PetscCall(MatGetSize(lM, &n, NULL));
2866: PetscCall(PetscContainerGetPointer(c, &csr));
2867: PetscCall(PCBDDCSetLocalAdjacencyGraph(bddcipc_ctx->bddc, n, csr, csr + (n + 1), PETSC_COPY_VALUES));
2868: PetscCall(MatISRestoreLocalMat(M, &lM));
2869: }
2870: PetscCall(PCSetOptionsPrefix(bddcipc_ctx->bddc, ((PetscObject)ksps[1])->prefix));
2871: PetscCall(PCSetErrorIfFailure(bddcipc_ctx->bddc, pc->erroriffailure));
2872: PetscCall(PCSetFromOptions(bddcipc_ctx->bddc));
2874: /* wrap the interface application */
2875: PetscCall(PCSetType(ppc, PCSHELL));
2876: PetscCall(PCShellSetName(ppc, "FETI-DP pressure"));
2877: PetscCall(PCShellSetContext(ppc, bddcipc_ctx));
2878: PetscCall(PCShellSetSetUp(ppc, PCSetUp_BDDCIPC));
2879: PetscCall(PCShellSetApply(ppc, PCApply_BDDCIPC));
2880: PetscCall(PCShellSetApplyTranspose(ppc, PCApplyTranspose_BDDCIPC));
2881: PetscCall(PCShellSetView(ppc, PCView_BDDCIPC));
2882: PetscCall(PCShellSetDestroy(ppc, PCDestroy_BDDCIPC));
2883: }
2885: /* determine if we need to assemble M to construct a preconditioner */
2886: if (!matisok) {
2887: PetscCall(PetscObjectTypeCompare((PetscObject)M, MATIS, &ismatis));
2888: PetscCall(PetscObjectTypeCompareAny((PetscObject)ppc, &matisok, PCBDDC, PCJACOBI, PCNONE, PCMG, ""));
2889: if (ismatis && !matisok) PetscCall(MatConvert(M, MATAIJ, MAT_INPLACE_MATRIX, &M));
2890: }
2892: /* run the subproblems to check convergence */
2893: PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)newmat)->prefix, "-check_saddlepoint", &check, NULL));
2894: if (check) {
2895: for (PetscInt i = 0; i < nn; i++) {
2896: KSP kspC;
2897: PC npc;
2898: Mat F, pF;
2899: Vec x, y;
2900: PetscBool isschur, prec = PETSC_TRUE;
2902: PetscCall(KSPCreate(PetscObjectComm((PetscObject)ksps[i]), &kspC));
2903: PetscCall(KSPSetNestLevel(kspC, pc->kspnestlevel));
2904: PetscCall(KSPSetOptionsPrefix(kspC, ((PetscObject)ksps[i])->prefix));
2905: PetscCall(KSPAppendOptionsPrefix(kspC, "check_"));
2906: PetscCall(KSPGetOperators(ksps[i], &F, &pF));
2907: PetscCall(PetscObjectTypeCompare((PetscObject)F, MATSCHURCOMPLEMENT, &isschur));
2908: if (isschur) {
2909: KSP kspS, kspS2;
2910: Mat A00, pA00, A10, A01, A11;
2911: char prefix[256];
2913: PetscCall(MatSchurComplementGetKSP(F, &kspS));
2914: PetscCall(MatSchurComplementGetSubMatrices(F, &A00, &pA00, &A01, &A10, &A11));
2915: PetscCall(MatCreateSchurComplement(A00, pA00, A01, A10, A11, &F));
2916: PetscCall(MatSchurComplementGetKSP(F, &kspS2));
2917: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sschur_", ((PetscObject)kspC)->prefix));
2918: PetscCall(KSPSetOptionsPrefix(kspS2, prefix));
2919: PetscCall(KSPGetPC(kspS2, &npc));
2920: PetscCall(PCSetType(npc, PCKSP));
2921: PetscCall(PCKSPSetKSP(npc, kspS));
2922: PetscCall(KSPSetFromOptions(kspS2));
2923: PetscCall(KSPGetPC(kspS2, &npc));
2924: PetscCall(PCSetUseAmat(npc, PETSC_TRUE));
2925: } else {
2926: PetscCall(PetscObjectReference((PetscObject)F));
2927: }
2928: PetscCall(KSPSetFromOptions(kspC));
2929: PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)kspC)->prefix, "-preconditioned", &prec, NULL));
2930: if (prec) {
2931: PetscCall(KSPGetPC(ksps[i], &npc));
2932: PetscCall(KSPSetPC(kspC, npc));
2933: }
2934: PetscCall(KSPSetOperators(kspC, F, pF));
2935: PetscCall(MatCreateVecs(F, &x, &y));
2936: PetscCall(VecSetRandom(x, NULL));
2937: PetscCall(MatMult(F, x, y));
2938: PetscCall(KSPSolve(kspC, y, x));
2939: PetscCall(KSPCheckSolve(kspC, npc, x));
2940: PetscCall(KSPDestroy(&kspC));
2941: PetscCall(MatDestroy(&F));
2942: PetscCall(VecDestroy(&x));
2943: PetscCall(VecDestroy(&y));
2944: }
2945: }
2946: PetscCall(PetscFree(ksps));
2947: }
2948: }
2949: /* return pointers for objects created */
2950: *fetidp_mat = newmat;
2951: *fetidp_pc = newpc;
2952: PetscFunctionReturn(PETSC_SUCCESS);
2953: }
2955: /*@
2956: PCBDDCCreateFETIDPOperators - Creates the FETI-DP matrix and its Dirichlet preconditioner
2958: Collective
2960: Input Parameters:
2961: + pc - the `PCBDDC` preconditioning context, after `PCSetUp()` has been called
2962: . fully_redundant - `PETSC_TRUE` for a fully redundant set of Lagrange multipliers
2963: - prefix - options database prefix for the objects to be created, or `NULL`
2965: Output Parameters:
2966: + fetidp_mat - the FETI-DP shell matrix
2967: - fetidp_pc - the shell Dirichlet preconditioner for the FETI-DP matrix
2969: Level: developer
2971: Notes:
2972: Most users should employ the `KSP` interface for linear solvers and create a solver of type `KSPFETIDP`.
2973: The FETI-DP matrix supports `MatMult()` and `MatMultTranspose()`.
2975: The caller must destroy the returned objects with `MatDestroy()` and `PCDestroy()`.
2977: .seealso: [](ch_ksp), `KSPFETIDP`, `PCBDDC`, `PCBDDCMatFETIDPGetRHS()`, `PCBDDCMatFETIDPGetSolution()`
2978: @*/
2979: PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, PetscBool fully_redundant, const char *prefix, Mat *fetidp_mat, PC *fetidp_pc)
2980: {
2981: PetscFunctionBegin;
2983: PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "You must call PCSetup_BDDC() first");
2984: PetscUseMethod(pc, "PCBDDCCreateFETIDPOperators_C", (PC, PetscBool, const char *, Mat *, PC *), (pc, fully_redundant, prefix, fetidp_mat, fetidp_pc));
2985: PetscFunctionReturn(PETSC_SUCCESS);
2986: }
2988: /*MC
2989: PCBDDC - Balancing Domain Decomposition by Constraints preconditioner
2991: Options Database Keys:
2992: + -pc_bddc_use_vertices (true|false) - include vertices in the primal space
2993: . -pc_bddc_use_edges (true|false) - include edge constraints in the primal space
2994: . -pc_bddc_use_faces (true|false) - include face constraints in the primal space
2995: . -pc_bddc_vertex_size size - classify connected components of at most this size as primal vertices
2996: . -pc_bddc_corner_selection (true|false) - select corners using subdomain faces and coordinates
2997: . -pc_bddc_use_local_mat_graph (true|false) - use the local matrix adjacency graph for interface analysis
2998: . -pc_bddc_local_mat_graph_square count - number of times to square the local matrix graph before interface analysis
2999: . -pc_bddc_graph_maxcount count - classify components shared by more than this many neighboring subdomains as primal vertices
3000: . -pc_bddc_detect_disconnected (true|false) - detect disconnected local subdomains
3001: . -pc_bddc_detect_disconnected_filter (true|false) - filter small local matrix entries when detecting disconnected subdomains
3002: . -pc_bddc_monolithic (true|false) - discard information about splitting degrees of freedom into fields
3003: . -pc_bddc_use_nnsp (true|false) - use the matrix near nullspace to construct constraints
3004: . -pc_bddc_use_nnsp_true (true|false) - use the supplied near-nullspace vectors directly, without orthonormalization
3005: . -pc_bddc_constraint_near_null_space_tol tol - discard near-nullspace vectors whose restriction to a connected component has norm at most tol
3006: . -pc_bddc_constraint_singular_tol tol - relative tolerance for retaining independent constraint modes
3007: . -pc_bddc_symmetric (true|false) - compute primal basis functions assuming symmetry; use false for nonsymmetric problems
3008: . -pc_bddc_use_change_of_basis (true|false) - construct a change of basis on edges
3009: . -pc_bddc_use_change_on_faces (true|false) - construct the requested change of basis on faces
3010: . -pc_bddc_interface_ext_type (dirichlet|lump) - select how interface corrections are extended to subdomain interiors
3011: . -pc_bddc_dirichlet_approximate (true|false) - enable nullspace corrections for approximate Dirichlet solvers
3012: . -pc_bddc_dirichlet_approximate_scale (true|false) - scale the approximate Dirichlet solver when applying nullspace corrections
3013: . -pc_bddc_neumann_approximate (true|false) - enable nullspace corrections for approximate Neumann solvers
3014: . -pc_bddc_neumann_approximate_scale (true|false) - scale the approximate Neumann solver when applying nullspace corrections
3015: . -pc_bddc_switch_static (true|false) - switch from the default $M_2$ operator to $M_3$ in {cite}`dohrmann2007approximate`
3016: . -pc_bddc_levels levels - maximum number of additional levels (default 0)
3017: . -pc_bddc_coarsening_ratio ratio - target number of process subdomains or local elements per aggregate (default 8)
3018: . -pc_bddc_coarse_eqs_per_proc neq - target number of equations per process at the coarsest level; a negative value uses one process
3019: . -pc_bddc_coarse_eqs_limit neq - stop adding coarse levels when the coarse problem has at most this many equations
3020: . -pc_bddc_coarse_adj nprocs - number of processes used to hold the coarse adjacency graph for partitioning
3021: . -pc_bddc_use_coarse_estimates (true|false) - use estimated eigenvalues to configure an iterative coarse solver
3022: . -pc_bddc_use_deluxe_scaling (true|false) - use deluxe scaling
3023: . -pc_bddc_deluxe_zerorows (true|false) - zero rows and columns of deluxe operators associated with primal degrees of freedom
3024: . -pc_bddc_deluxe_singlemat (true|false) - combine the deluxe scaling operations into one matrix per interface component
3025: . -pc_bddc_schur_rebuild (true|false) - rebuild the interface graph without adjacency information for computing Schur complement principal minors
3026: . -pc_bddc_schur_layers layers - number of layers used for economic deluxe scaling; -1 uses all degrees of freedom
3027: . -pc_bddc_schur_use_useradj (true|false) - use the graph supplied with `PCBDDCSetLocalAdjacencyGraph()` to select Schur complement layers
3028: . -pc_bddc_schur_exact (true|false) - use the full Schur complement, including components of size one, for adaptive constraint selection
3029: . -pc_bddc_adaptive_threshold thresholds - one or two comma-separated thresholds for adaptive constraint selection; one value sets both thresholds
3030: . -pc_bddc_adaptive_nmin count - minimum number of adaptive constraints per connected component
3031: . -pc_bddc_adaptive_nmax count - maximum number of adaptive constraints per connected component
3032: . -pc_bddc_adaptive_userdefined (true|false) - retain constraints from `MatSetNearNullSpace()` in addition to adaptive constraints
3033: . -pc_bddc_benign_trick (true|false) - use the benign subspace approach for saddle-point problems with discontinuous spaces
3034: . -pc_bddc_nonetflux (true|false) - compute quadrature weights for no-net-flux constraints automatically
3035: . -pc_bddc_nedelec_field_primal (true|false) - make Nedelec degrees of freedom shared by more than two subdomains primal
3036: . -pc_bddc_nedelec_order order - override the Nedelec order for testing; 0 selects variable order
3037: . -pc_bddc_nedelec_print (true|false) - print Nedelec setup diagnostics
3038: . -pc_bddc_load filename - load BDDC customization from a binary file for debugging
3039: . -pc_bddc_load_version version - version of the customization file to load
3040: . -pc_bddc_save filename - save BDDC customization to a binary file after setup for debugging
3041: . -pc_bddc_save_version version - version of the customization file to write
3042: . -pc_bddc_check_level level - verbosity level of debugging output
3044: Level: intermediate
3046: Notes:
3047: `PCBDDC` requires `MATIS` matrices and supports nonsymmetric and indefinite problems.
3048: The implementation and its customization are described in {cite}`zampini2016pcbddc`; see also [](sec_bddc).
3049: See {cite}`dohrmann2007approximate`, {cite}`klawonn2006dual`, and {cite}`mandel2008multispace` for the underlying methods.
3051: `PCBDDC` acts on all degrees of freedom, including subdomain interiors. This allows approximate subdomain solvers.
3052: Approximate local solvers are automatically adapted as described in {cite}`dohrmann2007approximate` when a nullspace object
3053: is attached to the subdomain matrices and approximate solvers are selected through the options database.
3055: Interface nodes are classified as vertices, edges, or faces using the local-to-global mapping of degrees of freedom
3056: and the local connectivity graph. The graph can be customized with `PCBDDCSetLocalAdjacencyGraph()`.
3057: Additional information about degrees of freedom can be supplied with `PCBDDCSetDofsSplitting()`, `PCBDDCSetDirichletBoundaries()`,
3058: `PCBDDCSetNeumannBoundaries()`, `PCBDDCSetPrimalVerticesIS()`, and their local counterparts.
3060: Support for $H(\mathrm{div})$ and $H(\mathrm{curl})$ problems is provided through `PCBDDCSetDivergenceMat()` and `PCBDDCSetDiscreteGradient()`.
3062: Constraints can be customized by attaching a `MatNullSpace` object to the `MATIS` matrix with `MatSetNearNullSpace()`.
3063: Linearly independent modes are retained using a singular value decomposition.
3065: When requested, a change of basis is performed as in {cite}`klawonn2006dual`. Local QR factorizations are used when more than
3066: one constraint is present on a connected component, such as an edge or a face. A user-defined change of basis can be supplied
3067: with `PCBDDCSetChangeOfBasisMat()`.
3069: Multilevel `PCBDDC` is supported as described in {cite}`mandel2008multispace`. Process subdomains are partitioned using a `MatPartitioning` object.
3070: When a local `MATIS` matrix stores multiple elements, their local coarse contributions are first aggregated using a `PetscPartitioner` object.
3071: In this case, the coarsening ratio is the target number of local elements per aggregate.
3073: Adaptive selection of primal constraints is supported for symmetric positive definite systems with high contrast in the coefficients
3074: when MUMPS or MKL_PARDISO is available. See {cite}`ohwidlundzampinidohrmann2017` for adaptive deluxe methods for Raviart-Thomas fields.
3075: The benign subspace approach for saddle-point problems with discontinuous spaces is described in {cite}`zampinitu2017`.
3077: Options for the Dirichlet, Neumann, coarse solver, and aggregation objects use the following prefixes, preceded by any user prefix:
3078: .vb
3079: -pc_bddc_dirichlet_
3080: -pc_bddc_neumann_
3081: -pc_bddc_coarse_
3082: -pc_bddc_aggregator_n_
3083: .ve
3084: For example, `-pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg` selects an approximate Dirichlet solver.
3085: By default, local solvers use `KSPPREONLY` with a direct factorization.
3086: At level `n`, the aggregator prefix is `pc_bddc_aggregator_n_`. Numeric-prefix fallback allows
3087: `-pc_bddc_aggregator_mat_partitioning_type type` and `-pc_bddc_aggregator_petscpartitioner_type type` to configure all levels.
3089: For BDDC level `N` > 0, the solver prefixes are:
3090: .vb
3091: -pc_bddc_dirichlet_lN_
3092: -pc_bddc_neumann_lN_
3093: -pc_bddc_coarse_lN_
3094: .ve
3095: Level 0 is the finest level. A coarse-level `PCBDDC` inherits the corresponding coarse-solver prefix. For example,
3096: .vb
3097: -pc_bddc_coarse_pc_bddc_adaptive_threshold 5
3098: .ve
3099: sets the adaptive constraint threshold to 5 for the first coarse-level `PCBDDC`.
3101: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `MATIS`, `KSPFETIDP`, `PCLU`, `PCGAMG`, `PCBDDCSetLocalAdjacencyGraph()`, `PCBDDCSetDofsSplitting()`,
3102: `PCBDDCSetDirichletBoundaries()`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCSetPrimalVerticesIS()`, `MatNullSpace`, `MatSetNearNullSpace()`,
3103: `PCBDDCSetChangeOfBasisMat()`, `PCBDDCSetDivergenceMat()`, `PCBDDCSetDiscreteGradient()`
3104: M*/
3106: PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc)
3107: {
3108: PC_BDDC *pcbddc;
3110: PetscFunctionBegin;
3111: PetscCall(PetscNew(&pcbddc));
3112: pc->data = pcbddc;
3114: PetscCall(PCISInitialize(pc));
3116: /* create local graph structure */
3117: PetscCall(PCBDDCGraphCreate(&pcbddc->mat_graph));
3119: /* BDDC nonzero defaults */
3120: pcbddc->use_nnsp = PETSC_TRUE;
3121: pcbddc->use_local_adj = PETSC_TRUE;
3122: pcbddc->use_vertices = PETSC_TRUE;
3123: pcbddc->use_edges = PETSC_TRUE;
3124: pcbddc->symmetric_primal = PETSC_TRUE;
3125: pcbddc->vertex_size = 1;
3126: pcbddc->recompute_topography = PETSC_TRUE;
3127: pcbddc->coarse_size = -1;
3128: pcbddc->use_exact_dirichlet_trick = PETSC_TRUE;
3129: pcbddc->coarsening_ratio = 8;
3130: pcbddc->coarse_eqs_per_proc = 1;
3131: pcbddc->benign_compute_correction = PETSC_TRUE;
3132: pcbddc->nedfield = -1;
3133: pcbddc->nedglobal = PETSC_TRUE;
3134: pcbddc->graphmaxcount = PETSC_INT_MAX;
3135: pcbddc->sub_schurs_layers = -1;
3136: pcbddc->adaptive_threshold[0] = 0.0;
3137: pcbddc->adaptive_threshold[1] = 0.0;
3139: /* function pointers */
3140: pc->ops->apply = PCApply_BDDC;
3141: pc->ops->applytranspose = PCApplyTranspose_BDDC;
3142: pc->ops->setup = PCSetUp_BDDC;
3143: pc->ops->destroy = PCDestroy_BDDC;
3144: pc->ops->setfromoptions = PCSetFromOptions_BDDC;
3145: pc->ops->view = PCView_BDDC;
3146: pc->ops->applyrichardson = NULL;
3147: pc->ops->applysymmetricleft = NULL;
3148: pc->ops->applysymmetricright = NULL;
3149: pc->ops->presolve = PCPreSolve_BDDC;
3150: pc->ops->postsolve = PCPostSolve_BDDC;
3151: pc->ops->reset = PCReset_BDDC;
3153: /* composing function */
3154: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDiscreteGradient_C", PCBDDCSetDiscreteGradient_BDDC));
3155: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDivergenceMat_C", PCBDDCSetDivergenceMat_BDDC));
3156: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetChangeOfBasisMat_C", PCBDDCSetChangeOfBasisMat_BDDC));
3157: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesLocalIS_C", PCBDDCSetPrimalVerticesLocalIS_BDDC));
3158: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesIS_C", PCBDDCSetPrimalVerticesIS_BDDC));
3159: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesLocalIS_C", PCBDDCGetPrimalVerticesLocalIS_BDDC));
3160: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesIS_C", PCBDDCGetPrimalVerticesIS_BDDC));
3161: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetCoarseningRatio_C", PCBDDCSetCoarseningRatio_BDDC));
3162: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevel_C", PCBDDCSetLevel_BDDC));
3163: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetUseExactDirichlet_C", PCBDDCSetUseExactDirichlet_BDDC));
3164: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevels_C", PCBDDCSetLevels_BDDC));
3165: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCLoadCustomization_C", PCBDDCLoadCustomization_BDDC));
3166: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSaveCustomization_C", PCBDDCSaveCustomization_BDDC));
3167: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundaries_C", PCBDDCSetDirichletBoundaries_BDDC));
3168: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundariesLocal_C", PCBDDCSetDirichletBoundariesLocal_BDDC));
3169: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundaries_C", PCBDDCSetNeumannBoundaries_BDDC));
3170: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundariesLocal_C", PCBDDCSetNeumannBoundariesLocal_BDDC));
3171: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundaries_C", PCBDDCGetDirichletBoundaries_BDDC));
3172: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundariesLocal_C", PCBDDCGetDirichletBoundariesLocal_BDDC));
3173: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundaries_C", PCBDDCGetNeumannBoundaries_BDDC));
3174: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundariesLocal_C", PCBDDCGetNeumannBoundariesLocal_BDDC));
3175: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplitting_C", PCBDDCSetDofsSplitting_BDDC));
3176: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplittingLocal_C", PCBDDCSetDofsSplittingLocal_BDDC));
3177: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", PCBDDCSetLocalAdjacencyGraph_BDDC));
3178: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCCreateFETIDPOperators_C", PCBDDCCreateFETIDPOperators_BDDC));
3179: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetRHS_C", PCBDDCMatFETIDPGetRHS_BDDC));
3180: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetSolution_C", PCBDDCMatFETIDPGetSolution_BDDC));
3181: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", PCPreSolveChangeRHS_BDDC));
3182: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_BDDC));
3183: PetscFunctionReturn(PETSC_SUCCESS);
3184: }
3186: /*@
3187: PCBDDCInitializePackage - Initializes the `PCBDDC` package
3189: Not Collective
3191: Level: developer
3193: Note:
3194: This routine is called by `PCInitializePackage()`.
3196: .seealso: [](ch_ksp), `PetscInitialize()`, `PCBDDCFinalizePackage()`
3197: @*/
3198: PetscErrorCode PCBDDCInitializePackage(void)
3199: {
3200: int i;
3202: PetscFunctionBegin;
3203: if (PCBDDCPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
3204: PCBDDCPackageInitialized = PETSC_TRUE;
3205: PetscCall(PetscRegisterFinalize(PCBDDCFinalizePackage));
3207: /* general events */
3208: PetscCall(PetscLogEventRegister("PCBDDCTopo", PC_CLASSID, &PC_BDDC_Topology[0]));
3209: PetscCall(PetscLogEventRegister("PCBDDCLKSP", PC_CLASSID, &PC_BDDC_LocalSolvers[0]));
3210: PetscCall(PetscLogEventRegister("PCBDDCLWor", PC_CLASSID, &PC_BDDC_LocalWork[0]));
3211: PetscCall(PetscLogEventRegister("PCBDDCCorr", PC_CLASSID, &PC_BDDC_CorrectionSetUp[0]));
3212: PetscCall(PetscLogEventRegister("PCBDDCASet", PC_CLASSID, &PC_BDDC_ApproxSetUp[0]));
3213: PetscCall(PetscLogEventRegister("PCBDDCAApp", PC_CLASSID, &PC_BDDC_ApproxApply[0]));
3214: PetscCall(PetscLogEventRegister("PCBDDCCSet", PC_CLASSID, &PC_BDDC_CoarseSetUp[0]));
3215: PetscCall(PetscLogEventRegister("PCBDDCCKSP", PC_CLASSID, &PC_BDDC_CoarseSolver[0]));
3216: PetscCall(PetscLogEventRegister("PCBDDCAdap", PC_CLASSID, &PC_BDDC_AdaptiveSetUp[0]));
3217: PetscCall(PetscLogEventRegister("PCBDDCScal", PC_CLASSID, &PC_BDDC_Scaling[0]));
3218: PetscCall(PetscLogEventRegister("PCBDDCSchr", PC_CLASSID, &PC_BDDC_Schurs[0]));
3219: PetscCall(PetscLogEventRegister("PCBDDCDirS", PC_CLASSID, &PC_BDDC_Solves[0][0]));
3220: PetscCall(PetscLogEventRegister("PCBDDCNeuS", PC_CLASSID, &PC_BDDC_Solves[0][1]));
3221: PetscCall(PetscLogEventRegister("PCBDDCCoaS", PC_CLASSID, &PC_BDDC_Solves[0][2]));
3222: for (i = 1; i < PETSC_PCBDDC_MAXLEVELS; i++) {
3223: char ename[32];
3225: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCTopo l%02d", i));
3226: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Topology[i]));
3227: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCLKSP l%02d", i));
3228: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_LocalSolvers[i]));
3229: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCLWor l%02d", i));
3230: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_LocalWork[i]));
3231: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCorr l%02d", i));
3232: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_CorrectionSetUp[i]));
3233: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCASet l%02d", i));
3234: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_ApproxSetUp[i]));
3235: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCAApp l%02d", i));
3236: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_ApproxApply[i]));
3237: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCSet l%02d", i));
3238: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_CoarseSetUp[i]));
3239: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCKSP l%02d", i));
3240: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_CoarseSolver[i]));
3241: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCAdap l%02d", i));
3242: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_AdaptiveSetUp[i]));
3243: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCScal l%02d", i));
3244: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Scaling[i]));
3245: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCSchr l%02d", i));
3246: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Schurs[i]));
3247: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCDirS l%02d", i));
3248: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Solves[i][0]));
3249: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCNeuS l%02d", i));
3250: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Solves[i][1]));
3251: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCoaS l%02d", i));
3252: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Solves[i][2]));
3253: }
3254: PetscFunctionReturn(PETSC_SUCCESS);
3255: }
3257: /*@
3258: PCBDDCFinalizePackage - Finalizes the `PCBDDC` package
3260: Not Collective
3262: Level: developer
3264: Note:
3265: This routine is called automatically by `PetscFinalize()`.
3267: .seealso: [](ch_ksp), `PetscFinalize()`, `PCBDDCInitializePackage()`
3268: @*/
3269: PetscErrorCode PCBDDCFinalizePackage(void)
3270: {
3271: PetscFunctionBegin;
3272: PCBDDCPackageInitialized = PETSC_FALSE;
3273: PetscFunctionReturn(PETSC_SUCCESS);
3274: }