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