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