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