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