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 IS zerodiag = NULL; 1531 PetscInt nrows,ncols; 1532 PetscBool computesubschurs; 1533 PetscBool computeconstraintsmatrix; 1534 PetscBool new_nearnullspace_provided,ismatis; 1535 PetscErrorCode ierr; 1536 1537 PetscFunctionBegin; 1538 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&ismatis);CHKERRQ(ierr); 1539 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PCBDDC preconditioner requires matrix of type MATIS"); 1540 ierr = MatGetSize(pc->pmat,&nrows,&ncols);CHKERRQ(ierr); 1541 if (nrows != ncols) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCBDDC preconditioner requires a square preconditioning matrix"); 1542 matis = (Mat_IS*)pc->pmat->data; 1543 /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */ 1544 /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup 1545 Also, BDDC builds its own KSP for the Dirichlet problem */ 1546 if (pc->setupcalled && pc->flag != SAME_NONZERO_PATTERN) pcbddc->recompute_topography = PETSC_TRUE; 1547 if (pcbddc->recompute_topography) { 1548 pcbddc->graphanalyzed = PETSC_FALSE; 1549 computeconstraintsmatrix = PETSC_TRUE; 1550 } else { 1551 computeconstraintsmatrix = PETSC_FALSE; 1552 } 1553 1554 /* check parameters' compatibility */ 1555 if (!pcbddc->use_deluxe_scaling) pcbddc->deluxe_zerorows = PETSC_FALSE; 1556 pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_threshold > 0.0); 1557 pcbddc->adaptive_userdefined = (PetscBool)(pcbddc->adaptive_selection && pcbddc->adaptive_userdefined); 1558 if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE; 1559 1560 computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling); 1561 if (pcbddc->switch_static) { 1562 PetscBool ismatis; 1563 ierr = PetscObjectTypeCompare((PetscObject)pc->mat,MATIS,&ismatis);CHKERRQ(ierr); 1564 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"When the static switch is one, the iteration matrix should be of type MATIS"); 1565 } 1566 1567 /* activate all connected components if the netflux has been requested */ 1568 if (pcbddc->compute_nonetflux) { 1569 pcbddc->use_vertices = PETSC_TRUE; 1570 pcbddc->use_edges = PETSC_TRUE; 1571 pcbddc->use_faces = PETSC_TRUE; 1572 } 1573 1574 /* Get stdout for dbg */ 1575 if (pcbddc->dbg_flag) { 1576 if (!pcbddc->dbg_viewer) { 1577 pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)); 1578 ierr = PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer);CHKERRQ(ierr); 1579 } 1580 ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1581 } 1582 1583 /* process topology information */ 1584 if (pcbddc->recompute_topography) { 1585 PetscContainer c; 1586 1587 ierr = PetscObjectQuery((PetscObject)pc->pmat,"_convert_nest_lfields",(PetscObject*)&c);CHKERRQ(ierr); 1588 if (c) { 1589 MatISLocalFields lf; 1590 ierr = PetscContainerGetPointer(c,(void**)&lf);CHKERRQ(ierr); 1591 ierr = PCBDDCSetDofsSplittingLocal(pc,lf->nr,lf->rf);CHKERRQ(ierr); 1592 } 1593 ierr = PCBDDCComputeLocalTopologyInfo(pc);CHKERRQ(ierr); 1594 /* detect local disconnected subdomains if requested (use matis->A) */ 1595 if (pcbddc->detect_disconnected) { 1596 PetscInt i; 1597 for (i=0;i<pcbddc->n_local_subs;i++) { 1598 ierr = ISDestroy(&pcbddc->local_subs[i]);CHKERRQ(ierr); 1599 } 1600 ierr = PetscFree(pcbddc->local_subs);CHKERRQ(ierr); 1601 ierr = MatDetectDisconnectedComponents(matis->A,PETSC_FALSE,&pcbddc->n_local_subs,&pcbddc->local_subs);CHKERRQ(ierr); 1602 } 1603 if (pcbddc->discretegradient) { 1604 ierr = PCBDDCNedelecSupport(pc);CHKERRQ(ierr); 1605 } 1606 } 1607 1608 /* change basis if requested by the user */ 1609 if (pcbddc->user_ChangeOfBasisMatrix) { 1610 /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */ 1611 pcbddc->use_change_of_basis = PETSC_FALSE; 1612 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 1613 } else { 1614 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1615 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1616 pcbddc->local_mat = matis->A; 1617 } 1618 1619 /* 1620 Compute change of basis on local pressures (aka zerodiag dofs) with the benign trick 1621 This should come earlier then PCISSetUp for extracting the correct subdomain matrices 1622 */ 1623 ierr = PCBDDCBenignShellMat(pc,PETSC_TRUE);CHKERRQ(ierr); 1624 if (pcbddc->benign_saddle_point) { 1625 PC_IS* pcis = (PC_IS*)pc->data; 1626 1627 if (pcbddc->user_ChangeOfBasisMatrix || pcbddc->use_change_of_basis || !computesubschurs) pcbddc->benign_change_explicit = PETSC_TRUE; 1628 /* detect local saddle point and change the basis in pcbddc->local_mat (TODO: reuse case) */ 1629 ierr = PCBDDCBenignDetectSaddlePoint(pc,&zerodiag);CHKERRQ(ierr); 1630 /* pop B0 mat from local mat */ 1631 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 1632 /* give pcis a hint to not reuse submatrices during PCISCreate */ 1633 if (pc->flag == SAME_NONZERO_PATTERN && pcis->reusesubmatrices == PETSC_TRUE) { 1634 if (pcbddc->benign_n && (pcbddc->benign_change_explicit || pcbddc->dbg_flag)) { 1635 pcis->reusesubmatrices = PETSC_FALSE; 1636 } else { 1637 pcis->reusesubmatrices = PETSC_TRUE; 1638 } 1639 } else { 1640 pcis->reusesubmatrices = PETSC_FALSE; 1641 } 1642 } 1643 1644 /* propagate relevant information -> TODO remove*/ 1645 #if !defined(PETSC_USE_COMPLEX) /* workaround for reals */ 1646 if (matis->A->symmetric_set) { 1647 ierr = MatSetOption(pcbddc->local_mat,MAT_HERMITIAN,matis->A->symmetric);CHKERRQ(ierr); 1648 } 1649 #endif 1650 if (matis->A->symmetric_set) { 1651 ierr = MatSetOption(pcbddc->local_mat,MAT_SYMMETRIC,matis->A->symmetric);CHKERRQ(ierr); 1652 } 1653 if (matis->A->spd_set) { 1654 ierr = MatSetOption(pcbddc->local_mat,MAT_SPD,matis->A->spd);CHKERRQ(ierr); 1655 } 1656 1657 /* Set up all the "iterative substructuring" common block without computing solvers */ 1658 { 1659 Mat temp_mat; 1660 1661 temp_mat = matis->A; 1662 matis->A = pcbddc->local_mat; 1663 ierr = PCISSetUp(pc,PETSC_FALSE);CHKERRQ(ierr); 1664 pcbddc->local_mat = matis->A; 1665 matis->A = temp_mat; 1666 } 1667 1668 /* Analyze interface */ 1669 if (!pcbddc->graphanalyzed) { 1670 ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr); 1671 computeconstraintsmatrix = PETSC_TRUE; 1672 if (pcbddc->adaptive_selection && !pcbddc->use_deluxe_scaling && !pcbddc->mat_graph->twodim) { 1673 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot compute the adaptive primal space for a problem with 3D edges without deluxe scaling"); 1674 } 1675 if (pcbddc->compute_nonetflux) { 1676 MatNullSpace nnfnnsp; 1677 1678 ierr = PCBDDCComputeNoNetFlux(pc->pmat,pcbddc->divudotp,pcbddc->divudotp_trans,pcbddc->divudotp_vl2l,pcbddc->mat_graph,&nnfnnsp);CHKERRQ(ierr); 1679 /* TODO what if a nearnullspace is already attached? */ 1680 ierr = MatSetNearNullSpace(pc->pmat,nnfnnsp);CHKERRQ(ierr); 1681 ierr = MatNullSpaceDestroy(&nnfnnsp);CHKERRQ(ierr); 1682 } 1683 } 1684 1685 /* check existence of a divergence free extension, i.e. 1686 b(v_I,p_0) = 0 for all v_I (raise error if not). 1687 Also, check that PCBDDCBenignGetOrSetP0 works */ 1688 if (pcbddc->benign_saddle_point && pcbddc->dbg_flag > 1) { 1689 ierr = PCBDDCBenignCheck(pc,zerodiag);CHKERRQ(ierr); 1690 } 1691 ierr = ISDestroy(&zerodiag);CHKERRQ(ierr); 1692 1693 /* Setup local dirichlet solver ksp_D and sub_schurs solvers */ 1694 if (computesubschurs && pcbddc->recompute_topography) { 1695 ierr = PCBDDCInitSubSchurs(pc);CHKERRQ(ierr); 1696 } 1697 /* SetUp Scaling operator (scaling matrices could be needed in SubSchursSetUp)*/ 1698 if (!pcbddc->use_deluxe_scaling) { 1699 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1700 } 1701 1702 /* finish setup solvers and do adaptive selection of constraints */ 1703 sub_schurs = pcbddc->sub_schurs; 1704 if (sub_schurs && sub_schurs->schur_explicit) { 1705 if (computesubschurs) { 1706 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1707 } 1708 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1709 } else { 1710 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1711 if (computesubschurs) { 1712 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1713 } 1714 } 1715 if (pcbddc->adaptive_selection) { 1716 ierr = PCBDDCAdaptiveSelection(pc);CHKERRQ(ierr); 1717 computeconstraintsmatrix = PETSC_TRUE; 1718 } 1719 1720 /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */ 1721 new_nearnullspace_provided = PETSC_FALSE; 1722 ierr = MatGetNearNullSpace(pc->pmat,&nearnullspace);CHKERRQ(ierr); 1723 if (pcbddc->onearnullspace) { /* already used nearnullspace */ 1724 if (!nearnullspace) { /* near null space attached to mat has been destroyed */ 1725 new_nearnullspace_provided = PETSC_TRUE; 1726 } else { 1727 /* determine if the two nullspaces are different (should be lightweight) */ 1728 if (nearnullspace != pcbddc->onearnullspace) { 1729 new_nearnullspace_provided = PETSC_TRUE; 1730 } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */ 1731 PetscInt i; 1732 const Vec *nearnullvecs; 1733 PetscObjectState state; 1734 PetscInt nnsp_size; 1735 ierr = MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 1736 for (i=0;i<nnsp_size;i++) { 1737 ierr = PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);CHKERRQ(ierr); 1738 if (pcbddc->onearnullvecs_state[i] != state) { 1739 new_nearnullspace_provided = PETSC_TRUE; 1740 break; 1741 } 1742 } 1743 } 1744 } 1745 } else { 1746 if (!nearnullspace) { /* both nearnullspaces are null */ 1747 new_nearnullspace_provided = PETSC_FALSE; 1748 } else { /* nearnullspace attached later */ 1749 new_nearnullspace_provided = PETSC_TRUE; 1750 } 1751 } 1752 1753 /* Setup constraints and related work vectors */ 1754 /* reset primal space flags */ 1755 pcbddc->new_primal_space = PETSC_FALSE; 1756 pcbddc->new_primal_space_local = PETSC_FALSE; 1757 if (computeconstraintsmatrix || new_nearnullspace_provided) { 1758 /* It also sets the primal space flags */ 1759 ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr); 1760 } 1761 /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */ 1762 ierr = PCBDDCSetUpLocalWorkVectors(pc);CHKERRQ(ierr); 1763 1764 if (pcbddc->use_change_of_basis) { 1765 PC_IS *pcis = (PC_IS*)(pc->data); 1766 1767 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 1768 if (pcbddc->benign_change) { 1769 ierr = MatDestroy(&pcbddc->benign_B0);CHKERRQ(ierr); 1770 /* pop B0 from pcbddc->local_mat */ 1771 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 1772 } 1773 /* get submatrices */ 1774 ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr); 1775 ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr); 1776 ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr); 1777 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);CHKERRQ(ierr); 1778 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);CHKERRQ(ierr); 1779 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);CHKERRQ(ierr); 1780 /* set flag in pcis to not reuse submatrices during PCISCreate */ 1781 pcis->reusesubmatrices = PETSC_FALSE; 1782 } else if (!pcbddc->user_ChangeOfBasisMatrix && !pcbddc->benign_change) { 1783 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1784 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1785 pcbddc->local_mat = matis->A; 1786 } 1787 /* SetUp coarse and local Neumann solvers */ 1788 ierr = PCBDDCSetUpSolvers(pc);CHKERRQ(ierr); 1789 /* SetUp Scaling operator */ 1790 if (pcbddc->use_deluxe_scaling) { 1791 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1792 } 1793 1794 /* mark topography as done */ 1795 pcbddc->recompute_topography = PETSC_FALSE; 1796 1797 /* wrap pcis->A_IB and pcis->A_BI if we did not change explicitly the variables on the pressures */ 1798 ierr = PCBDDCBenignShellMat(pc,PETSC_FALSE);CHKERRQ(ierr); 1799 1800 if (pcbddc->dbg_flag) { 1801 ierr = PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1802 } 1803 PetscFunctionReturn(0); 1804 } 1805 1806 /* 1807 PCApply_BDDC - Applies the BDDC operator to a vector. 1808 1809 Input Parameters: 1810 + pc - the preconditioner context 1811 - r - input vector (global) 1812 1813 Output Parameter: 1814 . z - output vector (global) 1815 1816 Application Interface Routine: PCApply() 1817 */ 1818 #undef __FUNCT__ 1819 #define __FUNCT__ "PCApply_BDDC" 1820 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z) 1821 { 1822 PC_IS *pcis = (PC_IS*)(pc->data); 1823 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1824 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1825 PetscErrorCode ierr; 1826 const PetscScalar one = 1.0; 1827 const PetscScalar m_one = -1.0; 1828 const PetscScalar zero = 0.0; 1829 1830 /* This code is similar to that provided in nn.c for PCNN 1831 NN interface preconditioner changed to BDDC 1832 Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */ 1833 1834 PetscFunctionBegin; 1835 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 1836 if (pcbddc->ChangeOfBasisMatrix) { 1837 Vec swap; 1838 1839 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);CHKERRQ(ierr); 1840 swap = pcbddc->work_change; 1841 pcbddc->work_change = r; 1842 r = swap; 1843 /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */ 1844 if (pcbddc->benign_apply_coarse_only && pcbddc->use_exact_dirichlet_trick && pcbddc->change_interior) { 1845 ierr = VecCopy(r,pcis->vec1_global);CHKERRQ(ierr); 1846 ierr = VecLockPush(pcis->vec1_global);CHKERRQ(ierr); 1847 } 1848 } 1849 if (pcbddc->benign_have_null) { /* get p0 from r */ 1850 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 1851 } 1852 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 1853 ierr = VecCopy(r,z);CHKERRQ(ierr); 1854 /* First Dirichlet solve */ 1855 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1856 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1857 /* 1858 Assembling right hand side for BDDC operator 1859 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1860 - pcis->vec1_B the interface part of the global vector z 1861 */ 1862 if (n_D) { 1863 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1864 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1865 if (pcbddc->switch_static) { 1866 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1867 1868 ierr = VecSet(pcis->vec1_N,0.);CHKERRQ(ierr); 1869 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1870 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1871 if (!pcbddc->switch_static_change) { 1872 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1873 } else { 1874 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1875 ierr = MatMult(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1876 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1877 } 1878 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1879 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1880 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1881 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1882 } else { 1883 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1884 } 1885 } else { 1886 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1887 } 1888 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1889 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1890 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1891 } else { 1892 if (!pcbddc->benign_apply_coarse_only) { 1893 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1894 } 1895 } 1896 1897 /* Apply interface preconditioner 1898 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1899 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);CHKERRQ(ierr); 1900 1901 /* Apply transpose of partition of unity operator */ 1902 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1903 1904 /* Second Dirichlet solve and assembling of output */ 1905 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1906 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1907 if (n_B) { 1908 if (pcbddc->switch_static) { 1909 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1910 1911 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1912 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1913 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1914 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1915 if (!pcbddc->switch_static_change) { 1916 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1917 } else { 1918 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1919 ierr = MatMult(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1920 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1921 } 1922 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1923 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1924 } else { 1925 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1926 } 1927 } else if (pcbddc->switch_static) { /* n_B is zero */ 1928 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1929 1930 if (!pcbddc->switch_static_change) { 1931 ierr = MatMult(matis->A,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 1932 } else { 1933 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);CHKERRQ(ierr); 1934 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1935 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);CHKERRQ(ierr); 1936 } 1937 } 1938 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1939 1940 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 1941 if (pcbddc->switch_static) { 1942 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 1943 } else { 1944 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 1945 } 1946 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1947 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1948 } else { 1949 if (pcbddc->switch_static) { 1950 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 1951 } else { 1952 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1953 } 1954 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1955 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1956 } 1957 if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */ 1958 if (pcbddc->benign_apply_coarse_only) { 1959 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 1960 } 1961 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 1962 } 1963 1964 if (pcbddc->ChangeOfBasisMatrix) { 1965 pcbddc->work_change = r; 1966 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 1967 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 1968 } 1969 PetscFunctionReturn(0); 1970 } 1971 1972 /* 1973 PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector. 1974 1975 Input Parameters: 1976 + pc - the preconditioner context 1977 - r - input vector (global) 1978 1979 Output Parameter: 1980 . z - output vector (global) 1981 1982 Application Interface Routine: PCApplyTranspose() 1983 */ 1984 #undef __FUNCT__ 1985 #define __FUNCT__ "PCApplyTranspose_BDDC" 1986 PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z) 1987 { 1988 PC_IS *pcis = (PC_IS*)(pc->data); 1989 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1990 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1991 PetscErrorCode ierr; 1992 const PetscScalar one = 1.0; 1993 const PetscScalar m_one = -1.0; 1994 const PetscScalar zero = 0.0; 1995 1996 PetscFunctionBegin; 1997 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 1998 if (pcbddc->ChangeOfBasisMatrix) { 1999 Vec swap; 2000 2001 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);CHKERRQ(ierr); 2002 swap = pcbddc->work_change; 2003 pcbddc->work_change = r; 2004 r = swap; 2005 /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */ 2006 if (pcbddc->benign_apply_coarse_only && pcbddc->exact_dirichlet_trick_app && pcbddc->change_interior) { 2007 ierr = VecCopy(r,pcis->vec1_global);CHKERRQ(ierr); 2008 ierr = VecLockPush(pcis->vec1_global);CHKERRQ(ierr); 2009 } 2010 } 2011 if (pcbddc->benign_have_null) { /* get p0 from r */ 2012 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 2013 } 2014 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 2015 ierr = VecCopy(r,z);CHKERRQ(ierr); 2016 /* First Dirichlet solve */ 2017 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2018 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2019 /* 2020 Assembling right hand side for BDDC operator 2021 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 2022 - pcis->vec1_B the interface part of the global vector z 2023 */ 2024 if (n_D) { 2025 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 2026 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 2027 if (pcbddc->switch_static) { 2028 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 2029 2030 ierr = VecSet(pcis->vec1_N,0.);CHKERRQ(ierr); 2031 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2032 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2033 if (!pcbddc->switch_static_change) { 2034 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2035 } else { 2036 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2037 ierr = MatMultTranspose(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 2038 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2039 } 2040 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2041 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2042 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2043 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2044 } else { 2045 ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 2046 } 2047 } else { 2048 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 2049 } 2050 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2051 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2052 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 2053 } else { 2054 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 2055 } 2056 2057 /* Apply interface preconditioner 2058 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 2059 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr); 2060 2061 /* Apply transpose of partition of unity operator */ 2062 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 2063 2064 /* Second Dirichlet solve and assembling of output */ 2065 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2066 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2067 if (n_B) { 2068 if (pcbddc->switch_static) { 2069 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 2070 2071 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2072 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2073 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2074 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2075 if (!pcbddc->switch_static_change) { 2076 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2077 } else { 2078 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2079 ierr = MatMultTranspose(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 2080 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2081 } 2082 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2083 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2084 } else { 2085 ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 2086 } 2087 } else if (pcbddc->switch_static) { /* n_B is zero */ 2088 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 2089 2090 if (!pcbddc->switch_static_change) { 2091 ierr = MatMultTranspose(matis->A,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 2092 } else { 2093 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);CHKERRQ(ierr); 2094 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2095 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);CHKERRQ(ierr); 2096 } 2097 } 2098 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 2099 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 2100 if (pcbddc->switch_static) { 2101 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 2102 } else { 2103 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 2104 } 2105 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2106 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2107 } else { 2108 if (pcbddc->switch_static) { 2109 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 2110 } else { 2111 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 2112 } 2113 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2114 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2115 } 2116 if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */ 2117 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 2118 } 2119 if (pcbddc->ChangeOfBasisMatrix) { 2120 pcbddc->work_change = r; 2121 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 2122 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 2123 } 2124 PetscFunctionReturn(0); 2125 } 2126 2127 #undef __FUNCT__ 2128 #define __FUNCT__ "PCDestroy_BDDC" 2129 PetscErrorCode PCDestroy_BDDC(PC pc) 2130 { 2131 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 2132 PetscErrorCode ierr; 2133 2134 PetscFunctionBegin; 2135 /* free BDDC custom data */ 2136 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 2137 /* destroy objects related to topography */ 2138 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 2139 /* free allocated graph structure */ 2140 ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr); 2141 /* free allocated sub schurs structure */ 2142 ierr = PetscFree(pcbddc->sub_schurs);CHKERRQ(ierr); 2143 /* destroy objects for scaling operator */ 2144 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 2145 ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr); 2146 /* free solvers stuff */ 2147 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 2148 /* free global vectors needed in presolve */ 2149 ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr); 2150 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 2151 /* free data created by PCIS */ 2152 ierr = PCISDestroy(pc);CHKERRQ(ierr); 2153 /* remove functions */ 2154 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDiscreteGradient_C",NULL);CHKERRQ(ierr); 2155 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",NULL);CHKERRQ(ierr); 2156 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);CHKERRQ(ierr); 2157 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 2158 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",NULL);CHKERRQ(ierr); 2159 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 2160 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr); 2161 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr); 2162 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr); 2163 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2164 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2165 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2166 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2167 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2168 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2169 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2170 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2171 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 2172 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr); 2173 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 2174 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 2175 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 2176 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 2177 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",NULL);CHKERRQ(ierr); 2178 /* Free the private data structure */ 2179 ierr = PetscFree(pc->data);CHKERRQ(ierr); 2180 PetscFunctionReturn(0); 2181 } 2182 2183 #undef __FUNCT__ 2184 #define __FUNCT__ "PCPreSolveChangeRHS_BDDC" 2185 static PetscErrorCode PCPreSolveChangeRHS_BDDC(PC pc, PetscBool* change) 2186 { 2187 PetscFunctionBegin; 2188 *change = PETSC_TRUE; 2189 PetscFunctionReturn(0); 2190 } 2191 2192 #undef __FUNCT__ 2193 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC" 2194 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 2195 { 2196 FETIDPMat_ctx mat_ctx; 2197 Vec work; 2198 PC_IS* pcis; 2199 PC_BDDC* pcbddc; 2200 PetscErrorCode ierr; 2201 2202 PetscFunctionBegin; 2203 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2204 pcis = (PC_IS*)mat_ctx->pc->data; 2205 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2206 2207 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 2208 /* copy rhs since we may change it during PCPreSolve_BDDC */ 2209 if (!pcbddc->original_rhs) { 2210 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 2211 } 2212 ierr = VecCopy(standard_rhs,pcbddc->original_rhs);CHKERRQ(ierr); 2213 if (mat_ctx->g2g_p) { 2214 /* interface pressure rhs */ 2215 ierr = VecScatterBegin(mat_ctx->g2g_p,fetidp_flux_rhs,pcbddc->original_rhs,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2216 ierr = VecScatterEnd(mat_ctx->g2g_p,fetidp_flux_rhs,pcbddc->original_rhs,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2217 ierr = VecScatterBegin(mat_ctx->g2g_p,standard_rhs,fetidp_flux_rhs,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2218 ierr = VecScatterEnd(mat_ctx->g2g_p,standard_rhs,fetidp_flux_rhs,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2219 ierr = VecScale(fetidp_flux_rhs,-1.);CHKERRQ(ierr); 2220 } 2221 /* 2222 change of basis for physical rhs if needed 2223 It also changes the rhs in case of dirichlet boundaries 2224 */ 2225 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,pcbddc->original_rhs,NULL);CHKERRQ(ierr); 2226 if (pcbddc->ChangeOfBasisMatrix) { 2227 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,pcbddc->original_rhs,pcbddc->work_change);CHKERRQ(ierr); 2228 work = pcbddc->work_change; 2229 } else { 2230 work = pcbddc->original_rhs; 2231 } 2232 /* store vectors for computation of fetidp final solution */ 2233 ierr = VecScatterBegin(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2234 ierr = VecScatterEnd(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2235 /* scale rhs since it should be unassembled */ 2236 /* TODO use counter scaling? (also below) */ 2237 ierr = VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2238 ierr = VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2239 /* Apply partition of unity */ 2240 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2241 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2242 if (!pcbddc->switch_static) { 2243 /* compute partially subassembled Schur complement right-hand side */ 2244 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2245 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 2246 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 2247 ierr = VecSet(work,0.0);CHKERRQ(ierr); 2248 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2249 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2250 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2251 ierr = VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2252 ierr = VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2253 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2254 } 2255 /* BDDC rhs */ 2256 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 2257 if (pcbddc->switch_static) { 2258 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2259 } 2260 /* apply BDDC */ 2261 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2262 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2263 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2264 2265 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 2266 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 2267 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2268 ierr = VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2269 /* Add contribution to interface pressures */ 2270 if (mat_ctx->l2g_p) { 2271 ierr = MatMult(mat_ctx->B_BB,pcis->vec1_B,mat_ctx->vP);CHKERRQ(ierr); 2272 if (pcbddc->switch_static) { 2273 ierr = MatMultAdd(mat_ctx->B_BI,pcis->vec1_D,mat_ctx->vP,mat_ctx->vP);CHKERRQ(ierr); 2274 } 2275 ierr = VecScatterBegin(mat_ctx->l2g_p,mat_ctx->vP,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2276 ierr = VecScatterEnd(mat_ctx->l2g_p,mat_ctx->vP,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2277 } 2278 PetscFunctionReturn(0); 2279 } 2280 2281 #undef __FUNCT__ 2282 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS" 2283 /*@ 2284 PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETI-DP linear system using the physical right-hand side 2285 2286 Collective 2287 2288 Input Parameters: 2289 + fetidp_mat - the FETI-DP matrix object obtained by a call to PCBDDCCreateFETIDPOperators 2290 - standard_rhs - the right-hand side of the original linear system 2291 2292 Output Parameters: 2293 . fetidp_flux_rhs - the right-hand side for the FETI-DP linear system 2294 2295 Level: developer 2296 2297 Notes: 2298 2299 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetSolution 2300 @*/ 2301 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 2302 { 2303 FETIDPMat_ctx mat_ctx; 2304 PetscErrorCode ierr; 2305 2306 PetscFunctionBegin; 2307 PetscValidHeaderSpecific(fetidp_mat,MAT_CLASSID,1); 2308 PetscValidHeaderSpecific(standard_rhs,VEC_CLASSID,2); 2309 PetscValidHeaderSpecific(fetidp_flux_rhs,VEC_CLASSID,3); 2310 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2311 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 2312 PetscFunctionReturn(0); 2313 } 2314 2315 #undef __FUNCT__ 2316 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC" 2317 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2318 { 2319 FETIDPMat_ctx mat_ctx; 2320 PC_IS* pcis; 2321 PC_BDDC* pcbddc; 2322 PetscErrorCode ierr; 2323 Vec work; 2324 2325 PetscFunctionBegin; 2326 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2327 pcis = (PC_IS*)mat_ctx->pc->data; 2328 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2329 2330 /* apply B_delta^T */ 2331 ierr = VecSet(pcis->vec1_B,0.);CHKERRQ(ierr); 2332 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2333 ierr = VecScatterEnd(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2334 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 2335 if (mat_ctx->l2g_p) { 2336 ierr = VecScatterBegin(mat_ctx->l2g_p,fetidp_flux_sol,mat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2337 ierr = VecScatterEnd(mat_ctx->l2g_p,fetidp_flux_sol,mat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2338 ierr = MatMultAdd(mat_ctx->Bt_BB,mat_ctx->vP,pcis->vec1_B,pcis->vec1_B);CHKERRQ(ierr); 2339 } 2340 2341 /* compute rhs for BDDC application */ 2342 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2343 if (pcbddc->switch_static) { 2344 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2345 if (mat_ctx->l2g_p) { 2346 ierr = VecScale(mat_ctx->vP,-1.);CHKERRQ(ierr); 2347 ierr = MatMultAdd(mat_ctx->Bt_BI,mat_ctx->vP,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); 2348 } 2349 } 2350 2351 /* apply BDDC */ 2352 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2353 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2354 2355 /* put values into global vector */ 2356 if (pcbddc->ChangeOfBasisMatrix) work = pcbddc->work_change; 2357 else work = standard_sol; 2358 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2359 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2360 if (!pcbddc->switch_static) { 2361 /* compute values into the interior if solved for the partially subassembled Schur complement */ 2362 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 2363 ierr = VecAYPX(pcis->vec1_D,-1.0,mat_ctx->temp_solution_D);CHKERRQ(ierr); 2364 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); 2365 } 2366 2367 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2368 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec1_D,work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2369 /* add p0 solution to final solution */ 2370 ierr = PCBDDCBenignGetOrSetP0(mat_ctx->pc,work,PETSC_FALSE);CHKERRQ(ierr); 2371 if (pcbddc->ChangeOfBasisMatrix) { 2372 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,work,standard_sol);CHKERRQ(ierr); 2373 } 2374 ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr); 2375 if (mat_ctx->g2g_p) { 2376 ierr = VecScale(fetidp_flux_sol,-1.);CHKERRQ(ierr); 2377 ierr = VecScatterBegin(mat_ctx->g2g_p,fetidp_flux_sol,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2378 ierr = VecScatterEnd(mat_ctx->g2g_p,fetidp_flux_sol,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2379 ierr = VecScale(fetidp_flux_sol,-1.);CHKERRQ(ierr); 2380 } 2381 PetscFunctionReturn(0); 2382 } 2383 2384 #undef __FUNCT__ 2385 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution" 2386 /*@ 2387 PCBDDCMatFETIDPGetSolution - Compute the physical solution using the solution of the FETI-DP linear system 2388 2389 Collective 2390 2391 Input Parameters: 2392 + fetidp_mat - the FETI-DP matrix obtained by a call to PCBDDCCreateFETIDPOperators 2393 - fetidp_flux_sol - the solution of the FETI-DP linear system 2394 2395 Output Parameters: 2396 . standard_sol - the solution defined on the physical domain 2397 2398 Level: developer 2399 2400 Notes: 2401 2402 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetRHS 2403 @*/ 2404 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2405 { 2406 FETIDPMat_ctx mat_ctx; 2407 PetscErrorCode ierr; 2408 2409 PetscFunctionBegin; 2410 PetscValidHeaderSpecific(fetidp_mat,MAT_CLASSID,1); 2411 PetscValidHeaderSpecific(fetidp_flux_sol,VEC_CLASSID,2); 2412 PetscValidHeaderSpecific(standard_sol,VEC_CLASSID,3); 2413 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2414 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 2415 PetscFunctionReturn(0); 2416 } 2417 2418 PETSC_EXTERN PetscErrorCode FETIDPMatMult(Mat,Vec,Vec); 2419 PETSC_EXTERN PetscErrorCode FETIDPMatMultTranspose(Mat,Vec,Vec); 2420 PETSC_EXTERN PetscErrorCode PCBDDCDestroyFETIDPMat(Mat); 2421 PETSC_EXTERN PetscErrorCode FETIDPPCApply(PC,Vec,Vec); 2422 PETSC_EXTERN PetscErrorCode FETIDPPCApplyTranspose(PC,Vec,Vec); 2423 PETSC_EXTERN PetscErrorCode PCBDDCDestroyFETIDPPC(PC); 2424 2425 #undef __FUNCT__ 2426 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC" 2427 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, PetscBool fully_redundant, Mat *fetidp_mat, PC *fetidp_pc) 2428 { 2429 2430 FETIDPMat_ctx fetidpmat_ctx; 2431 Mat newmat; 2432 FETIDPPC_ctx fetidppc_ctx; 2433 PC newpc; 2434 MPI_Comm comm; 2435 PetscErrorCode ierr; 2436 2437 PetscFunctionBegin; 2438 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 2439 /* FETIDP linear matrix */ 2440 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 2441 fetidpmat_ctx->fully_redundant = fully_redundant; 2442 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 2443 ierr = MatCreateShell(comm,fetidpmat_ctx->n,fetidpmat_ctx->n,fetidpmat_ctx->N,fetidpmat_ctx->N,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 2444 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 2445 ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr); 2446 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 2447 ierr = MatSetUp(newmat);CHKERRQ(ierr); 2448 /* FETIDP preconditioner */ 2449 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 2450 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 2451 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 2452 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 2453 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 2454 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 2455 ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 2456 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 2457 ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr); 2458 ierr = PCSetUp(newpc);CHKERRQ(ierr); 2459 /* return pointers for objects created */ 2460 *fetidp_mat=newmat; 2461 *fetidp_pc=newpc; 2462 PetscFunctionReturn(0); 2463 } 2464 2465 #undef __FUNCT__ 2466 #define __FUNCT__ "PCBDDCCreateFETIDPOperators" 2467 /*@ 2468 PCBDDCCreateFETIDPOperators - Create FETI-DP operators 2469 2470 Collective 2471 2472 Input Parameters: 2473 + pc - the BDDC preconditioning context (setup should have been called before) 2474 - fully_redundant - true for a fully redundant set of Lagrange multipliers 2475 2476 Output Parameters: 2477 + fetidp_mat - shell FETI-DP matrix object 2478 - fetidp_pc - shell Dirichlet preconditioner for FETI-DP matrix 2479 2480 Level: developer 2481 2482 Notes: 2483 Currently the only operations provided for FETI-DP matrix are MatMult and MatMultTranspose 2484 2485 .seealso: PCBDDC, PCBDDCMatFETIDPGetRHS, PCBDDCMatFETIDPGetSolution 2486 @*/ 2487 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, PetscBool fully_redundant, Mat *fetidp_mat, PC *fetidp_pc) 2488 { 2489 PetscErrorCode ierr; 2490 2491 PetscFunctionBegin; 2492 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2493 if (pc->setupcalled) { 2494 ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,PetscBool,Mat*,PC*),(pc,fully_redundant,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 2495 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n"); 2496 PetscFunctionReturn(0); 2497 } 2498 /* -------------------------------------------------------------------------- */ 2499 /*MC 2500 PCBDDC - Balancing Domain Decomposition by Constraints. 2501 2502 An implementation of the BDDC preconditioner based on 2503 2504 .vb 2505 [1] C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 2506 [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 2507 [3] J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", http://arxiv.org/abs/0712.3977 2508 [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 2509 .ve 2510 2511 The matrix to be preconditioned (Pmat) must be of type MATIS. 2512 2513 Currently works with MATIS matrices with local matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers. 2514 2515 It also works with unsymmetric and indefinite problems. 2516 2517 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. 2518 2519 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). 2520 2521 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() 2522 Additional information on dofs can be provided by using PCBDDCSetDofsSplitting(), PCBDDCSetDirichletBoundaries(), PCBDDCSetNeumannBoundaries(), and PCBDDCSetPrimalVerticesIS() and their local counterparts. 2523 2524 Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace(). Non-singular modes are retained via SVD. 2525 2526 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. 2527 User defined change of basis can be passed to PCBDDC by using PCBDDCSetChangeOfBasisMat() 2528 2529 The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using a MatPartitioning object. 2530 2531 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. 2532 2533 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. 2534 Deluxe scaling is not supported yet for FETI-DP. 2535 2536 Options Database Keys (some of them, run with -h for a complete list): 2537 2538 . -pc_bddc_use_vertices <true> - use or not vertices in primal space 2539 . -pc_bddc_use_edges <true> - use or not edges in primal space 2540 . -pc_bddc_use_faces <false> - use or not faces in primal space 2541 . -pc_bddc_symmetric <true> - symmetric computation of primal basis functions. Specify false for unsymmetric problems 2542 . -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only) 2543 . -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested 2544 . -pc_bddc_switch_static <false> - switches from M_2 (default) to M_3 operator (see reference article [1]) 2545 . -pc_bddc_levels <0> - maximum number of levels for multilevel 2546 . -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) 2547 . -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) 2548 . -pc_bddc_use_deluxe_scaling <false> - use deluxe scaling 2549 . -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) 2550 . -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) 2551 - -pc_bddc_check_level <0> - set verbosity level of debugging output 2552 2553 Options for Dirichlet, Neumann or coarse solver can be set with 2554 .vb 2555 -pc_bddc_dirichlet_ 2556 -pc_bddc_neumann_ 2557 -pc_bddc_coarse_ 2558 .ve 2559 e.g -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. PCBDDC uses by default KPSPREONLY and PCLU. 2560 2561 When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified as 2562 .vb 2563 -pc_bddc_dirichlet_lN_ 2564 -pc_bddc_neumann_lN_ 2565 -pc_bddc_coarse_lN_ 2566 .ve 2567 Note that level number ranges from the finest (0) to the coarsest (N). 2568 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. 2569 .vb 2570 -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3 2571 .ve 2572 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 2573 2574 Level: intermediate 2575 2576 Developer notes: 2577 2578 Contributed by Stefano Zampini 2579 2580 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 2581 M*/ 2582 2583 #undef __FUNCT__ 2584 #define __FUNCT__ "PCCreate_BDDC" 2585 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 2586 { 2587 PetscErrorCode ierr; 2588 PC_BDDC *pcbddc; 2589 2590 PetscFunctionBegin; 2591 /* Creates the private data structure for this preconditioner and attach it to the PC object. */ 2592 ierr = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr); 2593 pc->data = (void*)pcbddc; 2594 2595 /* create PCIS data structure */ 2596 ierr = PCISCreate(pc);CHKERRQ(ierr); 2597 2598 /* BDDC customization */ 2599 pcbddc->use_local_adj = PETSC_TRUE; 2600 pcbddc->use_vertices = PETSC_TRUE; 2601 pcbddc->use_edges = PETSC_TRUE; 2602 pcbddc->use_faces = PETSC_FALSE; 2603 pcbddc->use_change_of_basis = PETSC_FALSE; 2604 pcbddc->use_change_on_faces = PETSC_FALSE; 2605 pcbddc->switch_static = PETSC_FALSE; 2606 pcbddc->use_nnsp_true = PETSC_FALSE; 2607 pcbddc->use_qr_single = PETSC_FALSE; 2608 pcbddc->symmetric_primal = PETSC_TRUE; 2609 pcbddc->benign_saddle_point = PETSC_FALSE; 2610 pcbddc->benign_have_null = PETSC_FALSE; 2611 pcbddc->vertex_size = 1; 2612 pcbddc->dbg_flag = 0; 2613 /* private */ 2614 pcbddc->local_primal_size = 0; 2615 pcbddc->local_primal_size_cc = 0; 2616 pcbddc->local_primal_ref_node = 0; 2617 pcbddc->local_primal_ref_mult = 0; 2618 pcbddc->n_vertices = 0; 2619 pcbddc->primal_indices_local_idxs = 0; 2620 pcbddc->recompute_topography = PETSC_TRUE; 2621 pcbddc->coarse_size = -1; 2622 pcbddc->new_primal_space = PETSC_FALSE; 2623 pcbddc->new_primal_space_local = PETSC_FALSE; 2624 pcbddc->global_primal_indices = 0; 2625 pcbddc->onearnullspace = 0; 2626 pcbddc->onearnullvecs_state = 0; 2627 pcbddc->user_primal_vertices = 0; 2628 pcbddc->user_primal_vertices_local = 0; 2629 pcbddc->temp_solution = 0; 2630 pcbddc->original_rhs = 0; 2631 pcbddc->local_mat = 0; 2632 pcbddc->ChangeOfBasisMatrix = 0; 2633 pcbddc->user_ChangeOfBasisMatrix = 0; 2634 pcbddc->coarse_vec = 0; 2635 pcbddc->coarse_ksp = 0; 2636 pcbddc->coarse_phi_B = 0; 2637 pcbddc->coarse_phi_D = 0; 2638 pcbddc->coarse_psi_B = 0; 2639 pcbddc->coarse_psi_D = 0; 2640 pcbddc->vec1_P = 0; 2641 pcbddc->vec1_R = 0; 2642 pcbddc->vec2_R = 0; 2643 pcbddc->local_auxmat1 = 0; 2644 pcbddc->local_auxmat2 = 0; 2645 pcbddc->R_to_B = 0; 2646 pcbddc->R_to_D = 0; 2647 pcbddc->ksp_D = 0; 2648 pcbddc->ksp_R = 0; 2649 pcbddc->NeumannBoundaries = 0; 2650 pcbddc->NeumannBoundariesLocal = 0; 2651 pcbddc->DirichletBoundaries = 0; 2652 pcbddc->DirichletBoundariesLocal = 0; 2653 pcbddc->user_provided_isfordofs = PETSC_FALSE; 2654 pcbddc->n_ISForDofs = 0; 2655 pcbddc->n_ISForDofsLocal = 0; 2656 pcbddc->ISForDofs = 0; 2657 pcbddc->ISForDofsLocal = 0; 2658 pcbddc->ConstraintMatrix = 0; 2659 pcbddc->use_exact_dirichlet_trick = PETSC_TRUE; 2660 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 2661 pcbddc->coarse_loc_to_glob = 0; 2662 pcbddc->coarsening_ratio = 8; 2663 pcbddc->coarse_adj_red = 0; 2664 pcbddc->current_level = 0; 2665 pcbddc->max_levels = 0; 2666 pcbddc->use_coarse_estimates = PETSC_FALSE; 2667 pcbddc->coarse_eqs_per_proc = 1; 2668 pcbddc->coarse_subassembling = 0; 2669 pcbddc->detect_disconnected = PETSC_FALSE; 2670 pcbddc->n_local_subs = 0; 2671 pcbddc->local_subs = NULL; 2672 2673 /* benign subspace trick */ 2674 pcbddc->benign_change = 0; 2675 pcbddc->benign_compute_correction = PETSC_TRUE; 2676 pcbddc->benign_vec = 0; 2677 pcbddc->benign_original_mat = 0; 2678 pcbddc->benign_sf = 0; 2679 pcbddc->benign_B0 = 0; 2680 pcbddc->benign_n = 0; 2681 pcbddc->benign_p0 = NULL; 2682 pcbddc->benign_p0_lidx = NULL; 2683 pcbddc->benign_p0_gidx = NULL; 2684 pcbddc->benign_null = PETSC_FALSE; 2685 2686 /* Nedelec */ 2687 pcbddc->nedfield = -1; 2688 pcbddc->nedglobal = PETSC_TRUE; 2689 2690 /* create local graph structure */ 2691 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 2692 pcbddc->graphmaxcount = PETSC_MAX_INT; 2693 2694 /* scaling */ 2695 pcbddc->work_scaling = 0; 2696 pcbddc->use_deluxe_scaling = PETSC_FALSE; 2697 2698 /* sub schurs options */ 2699 pcbddc->sub_schurs_rebuild = PETSC_FALSE; 2700 pcbddc->sub_schurs_layers = -1; 2701 pcbddc->sub_schurs_use_useradj = PETSC_FALSE; 2702 2703 pcbddc->computed_rowadj = PETSC_FALSE; 2704 2705 /* adaptivity */ 2706 pcbddc->adaptive_threshold = 0.0; 2707 pcbddc->adaptive_nmax = 0; 2708 pcbddc->adaptive_nmin = 0; 2709 2710 /* function pointers */ 2711 pc->ops->apply = PCApply_BDDC; 2712 pc->ops->applytranspose = PCApplyTranspose_BDDC; 2713 pc->ops->setup = PCSetUp_BDDC; 2714 pc->ops->destroy = PCDestroy_BDDC; 2715 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 2716 pc->ops->view = PCView_BDDC; 2717 pc->ops->applyrichardson = 0; 2718 pc->ops->applysymmetricleft = 0; 2719 pc->ops->applysymmetricright = 0; 2720 pc->ops->presolve = PCPreSolve_BDDC; 2721 pc->ops->postsolve = PCPostSolve_BDDC; 2722 2723 /* composing function */ 2724 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDiscreteGradient_C",PCBDDCSetDiscreteGradient_BDDC);CHKERRQ(ierr); 2725 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",PCBDDCSetDivergenceMat_BDDC);CHKERRQ(ierr); 2726 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);CHKERRQ(ierr); 2727 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 2728 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",PCBDDCSetPrimalVerticesIS_BDDC);CHKERRQ(ierr); 2729 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 2730 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr); 2731 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr); 2732 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr); 2733 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2734 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2735 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2736 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2737 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2738 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2739 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2740 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2741 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 2742 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr); 2743 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 2744 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 2745 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 2746 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 2747 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",PCPreSolveChangeRHS_BDDC);CHKERRQ(ierr); 2748 PetscFunctionReturn(0); 2749 } 2750 2751