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