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