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