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