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