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