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