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