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