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