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