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