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