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