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