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