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