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