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