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