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