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