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