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: Graph max count: %d\n",pcbddc->graphmaxcount);CHKERRQ(ierr); 120 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use vertices: %d (vertex size %d)\n",pcbddc->use_vertices,pcbddc->vertex_size);CHKERRQ(ierr); 121 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use edges: %d\n",pcbddc->use_edges);CHKERRQ(ierr); 122 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use faces: %d\n",pcbddc->use_faces);CHKERRQ(ierr); 123 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use true near null space: %d\n",pcbddc->use_nnsp_true);CHKERRQ(ierr); 124 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use QR for single constraints on cc: %d\n",pcbddc->use_qr_single);CHKERRQ(ierr); 125 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use change of basis on local edge nodes: %d\n",pcbddc->use_change_of_basis);CHKERRQ(ierr); 126 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use change of basis on local face nodes: %d\n",pcbddc->use_change_on_faces);CHKERRQ(ierr); 127 ierr = PetscViewerASCIIPrintf(viewer," BDDC: User defined change of basis matrix: %d\n",!!pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 128 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Has change of basis matrix: %d\n",!!pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 129 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Eliminate dirichlet boundary dofs: %d\n",pcbddc->eliminate_dirdofs);CHKERRQ(ierr); 130 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Switch on static condensation ops around the interface preconditioner: %d\n",pcbddc->switch_static);CHKERRQ(ierr); 131 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use exact dirichlet trick: %d\n",pcbddc->use_exact_dirichlet_trick);CHKERRQ(ierr); 132 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Multilevel max levels: %d\n",pcbddc->max_levels);CHKERRQ(ierr); 133 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Multilevel coarsening ratio: %d\n",pcbddc->coarsening_ratio);CHKERRQ(ierr); 134 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use estimated eigs for coarse problem: %d\n",pcbddc->use_coarse_estimates);CHKERRQ(ierr); 135 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use deluxe scaling: %d\n",pcbddc->use_deluxe_scaling);CHKERRQ(ierr); 136 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use deluxe zerorows: %d\n",pcbddc->deluxe_zerorows);CHKERRQ(ierr); 137 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Rebuild interface graph for Schur principal minors: %d\n",pcbddc->sub_schurs_rebuild);CHKERRQ(ierr); 138 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Number of dofs' layers for the computation of principal minors: %d\n",pcbddc->sub_schurs_layers);CHKERRQ(ierr); 139 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use user CSR graph to compute successive layers: %d\n",pcbddc->sub_schurs_use_useradj);CHKERRQ(ierr); 140 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); 141 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Min constraints / connected component: %d\n",pcbddc->adaptive_nmin);CHKERRQ(ierr); 142 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Max constraints / connected component: %d\n",pcbddc->adaptive_nmax);CHKERRQ(ierr); 143 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Invert exact Schur complement for adaptive selection: %d\n",pcbddc->sub_schurs_exact_schur);CHKERRQ(ierr); 144 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Symmetric computation of primal basis functions: %d\n",pcbddc->symmetric_primal);CHKERRQ(ierr); 145 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Num. Procs. to map coarse adjacency list: %d\n",pcbddc->coarse_adj_red);CHKERRQ(ierr); 146 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Coarse eqs per proc (significant at the coarsest level): %d\n",pcbddc->coarse_eqs_per_proc);CHKERRQ(ierr); 147 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Detect disconnected: %d\n",pcbddc->detect_disconnected);CHKERRQ(ierr); 148 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Benign subspace trick: %d (change explicit %d)\n",pcbddc->benign_saddle_point,pcbddc->benign_change_explicit);CHKERRQ(ierr); 149 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Benign subspace trick is active: %d\n",pcbddc->benign_have_null);CHKERRQ(ierr); 150 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Algebraic computation of no-net-flux %d\n",pcbddc->compute_nonetflux);CHKERRQ(ierr); 151 152 /* compute some numbers on the domain decomposition */ 153 ierr = VecSet(pcis->vec1_B,1.0);CHKERRQ(ierr); 154 ierr = MatCreateVecs(pc->pmat,&counter,0);CHKERRQ(ierr); 155 ierr = VecSet(counter,0.0);CHKERRQ(ierr); 156 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,counter,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 157 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,counter,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 158 ierr = VecSum(counter,&interface_size);CHKERRQ(ierr); 159 ierr = VecDestroy(&counter);CHKERRQ(ierr); 160 gsum[0] = 1; 161 gsum[1] = gsum[2] = gsum[3] = gsum[4] = 0; 162 loc[0] = !!pcis->n; 163 loc[1] = pcis->n - pcis->n_B; 164 loc[2] = pcis->n_B; 165 loc[3] = pcbddc->local_primal_size; 166 loc[4] = pcis->n; 167 loc[5] = pcbddc->benign_n; 168 ierr = MPI_Reduce(loc,gsum,5,MPIU_INT64,MPI_SUM,0,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 169 if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = -1; 170 ierr = MPI_Reduce(loc,gmax,5,MPIU_INT64,MPI_MAX,0,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 171 if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = PETSC_MAX_INT; 172 ierr = MPI_Reduce(loc,gmin,5,MPIU_INT64,MPI_MIN,0,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 173 ierr = MPI_Reduce(&loc[5],&totbenign,1,MPIU_INT64,MPI_SUM,0,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 174 if (pcbddc->coarse_size) { 175 ratio1 = pc->pmat->rmap->N/(1.*pcbddc->coarse_size); 176 ratio2 = PetscRealPart(interface_size)/pcbddc->coarse_size; 177 } 178 ierr = PetscViewerASCIIPrintf(viewer," BDDC: ********************************** STATISTICS AT LEVEL %d **********************************\n",pcbddc->current_level);CHKERRQ(ierr); 179 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); 180 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Coarsening ratios: all/coarse %d interface/coarse %d\n",(PetscInt)ratio1,(PetscInt)ratio2);CHKERRQ(ierr); 181 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Active processes : %d\n",(PetscInt)gsum[0]);CHKERRQ(ierr); 182 if (pcbddc->benign_have_null) { 183 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Benign subs : %d\n",(PetscInt)totbenign);CHKERRQ(ierr); 184 } 185 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Dofs type :\tMIN\tMAX\tMEAN\n",(PetscInt)gmin[1],(PetscInt)gmax[1],(PetscInt)(gsum[1]/gsum[0]));CHKERRQ(ierr); 186 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); 187 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); 188 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); 189 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); 190 ierr = PetscViewerASCIIPrintf(viewer," BDDC: ********************************** COARSE PROBLEM DETAILS *********************************\n",pcbddc->current_level);CHKERRQ(ierr); 191 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 192 193 /* the coarse problem can be handled with a different communicator */ 194 if (pcbddc->coarse_ksp) color = 1; 195 else color = 0; 196 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);CHKERRQ(ierr); 197 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr); 198 ierr = PetscSubcommCreate(PetscObjectComm((PetscObject)pc),&subcomm);CHKERRQ(ierr); 199 ierr = PetscSubcommSetNumber(subcomm,PetscMin(size,2));CHKERRQ(ierr); 200 ierr = PetscSubcommSetTypeGeneral(subcomm,color,rank);CHKERRQ(ierr); 201 ierr = PetscViewerGetSubViewer(viewer,PetscSubcommChild(subcomm),&subviewer);CHKERRQ(ierr); 202 if (color == 1) { 203 ierr = KSPView(pcbddc->coarse_ksp,subviewer);CHKERRQ(ierr); 204 ierr = PetscViewerFlush(subviewer);CHKERRQ(ierr); 205 } 206 ierr = PetscViewerRestoreSubViewer(viewer,PetscSubcommChild(subcomm),&subviewer);CHKERRQ(ierr); 207 ierr = PetscSubcommDestroy(&subcomm);CHKERRQ(ierr); 208 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 209 } 210 PetscFunctionReturn(0); 211 } 212 213 #undef __FUNCT__ 214 #define __FUNCT__ "PCBDDCSetDiscreteGradient_BDDC" 215 static PetscErrorCode PCBDDCSetDiscreteGradient_BDDC(PC pc, Mat G, PetscInt order, PetscInt field, PetscBool conforming) 216 { 217 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 218 PetscErrorCode ierr; 219 220 PetscFunctionBegin; 221 ierr = PetscObjectReference((PetscObject)G);CHKERRQ(ierr); 222 ierr = MatDestroy(&pcbddc->discretegradient);CHKERRQ(ierr); 223 pcbddc->discretegradient = G; 224 pcbddc->nedorder = order > 0 ? order : -order; 225 pcbddc->conforming = conforming; 226 pcbddc->nedfield = field; 227 PetscFunctionReturn(0); 228 } 229 230 #undef __FUNCT__ 231 #define __FUNCT__ "PCBDDCSetDiscreteGradient" 232 /*@ 233 PCBDDCSetDiscreteGradient - Sets the discrete gradient 234 235 Collective on PC 236 237 Input Parameters: 238 + pc - the preconditioning context 239 . G - the discrete gradient matrix (should be in AIJ format) 240 . order - the order of the Nedelec space (1 for the lowest order) 241 . field - the field id of the Nedelec dofs (not used if the fields have not been specified) 242 - conforming - whether the mesh is conforming or not 243 244 Level: advanced 245 246 Notes: The discrete gradient matrix G is used to analyze the subdomain edges, and it should not contain any zero entry. 247 For variable order spaces, the order should be set to zero. 248 The discrete gradient matrix is modified by PCBDDC. 249 250 .seealso: PCBDDC,PCBDDCSetDofsSplitting(),PCBDDCSetDofsSplittingLocal() 251 @*/ 252 PetscErrorCode PCBDDCSetDiscreteGradient(PC pc, Mat G, PetscInt order, PetscInt field, PetscBool conforming) 253 { 254 PetscErrorCode ierr; 255 256 PetscFunctionBegin; 257 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 258 PetscValidHeaderSpecific(G,MAT_CLASSID,2); 259 PetscCheckSameComm(pc,1,G,2); 260 PetscValidLogicalCollectiveInt(pc,order,3); 261 PetscValidLogicalCollectiveBool(pc,conforming,4); 262 ierr = PetscTryMethod(pc,"PCBDDCSetDiscreteGradient_C",(PC,Mat,PetscInt,PetscInt,PetscBool),(pc,G,order,field,conforming));CHKERRQ(ierr); 263 PetscFunctionReturn(0); 264 } 265 266 #undef __FUNCT__ 267 #define __FUNCT__ "PCBDDCSetDivergenceMat_BDDC" 268 static PetscErrorCode PCBDDCSetDivergenceMat_BDDC(PC pc, Mat divudotp, PetscBool trans, IS vl2l) 269 { 270 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 271 PetscErrorCode ierr; 272 273 PetscFunctionBegin; 274 ierr = PetscObjectReference((PetscObject)divudotp);CHKERRQ(ierr); 275 ierr = MatDestroy(&pcbddc->divudotp);CHKERRQ(ierr); 276 pcbddc->divudotp = divudotp; 277 pcbddc->divudotp_trans = trans; 278 pcbddc->compute_nonetflux = PETSC_TRUE; 279 if (vl2l) { 280 ierr = PetscObjectReference((PetscObject)vl2l);CHKERRQ(ierr); 281 ierr = ISDestroy(&pcbddc->divudotp_vl2l);CHKERRQ(ierr); 282 pcbddc->divudotp_vl2l = vl2l; 283 } 284 PetscFunctionReturn(0); 285 } 286 287 #undef __FUNCT__ 288 #define __FUNCT__ "PCBDDCSetDivergenceMat" 289 /*@ 290 PCBDDCSetDivergenceMat - Sets the linear operator representing \int_\Omega \div {\bf u} \cdot p dx 291 292 Collective on PC 293 294 Input Parameters: 295 + pc - the preconditioning context 296 . divudotp - the matrix (must be of type MATIS) 297 . trans - if trans if false (resp. true), then pressures are in the test (trial) space and velocities are in the trial (test) space. 298 - 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 299 300 Level: advanced 301 302 Notes: This auxiliary matrix is used to compute quadrature weights representing the net-flux across subdomain boundaries 303 304 .seealso: PCBDDC 305 @*/ 306 PetscErrorCode PCBDDCSetDivergenceMat(PC pc, Mat divudotp, PetscBool trans, IS vl2l) 307 { 308 PetscBool ismatis; 309 PetscErrorCode ierr; 310 311 PetscFunctionBegin; 312 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 313 PetscValidHeaderSpecific(divudotp,MAT_CLASSID,2); 314 PetscCheckSameComm(pc,1,divudotp,2); 315 PetscValidLogicalCollectiveBool(pc,trans,3); 316 if (vl2l) PetscValidHeaderSpecific(divudotp,IS_CLASSID,4); 317 ierr = PetscObjectTypeCompare((PetscObject)divudotp,MATIS,&ismatis);CHKERRQ(ierr); 318 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Divergence matrix needs to be of type MATIS"); 319 ierr = PetscTryMethod(pc,"PCBDDCSetDivergenceMat_C",(PC,Mat,PetscBool,IS),(pc,divudotp,trans,vl2l));CHKERRQ(ierr); 320 PetscFunctionReturn(0); 321 } 322 323 #undef __FUNCT__ 324 #define __FUNCT__ "PCBDDCSetChangeOfBasisMat_BDDC" 325 static PetscErrorCode PCBDDCSetChangeOfBasisMat_BDDC(PC pc, Mat change, PetscBool interior) 326 { 327 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 328 PetscErrorCode ierr; 329 330 PetscFunctionBegin; 331 ierr = PetscObjectReference((PetscObject)change);CHKERRQ(ierr); 332 ierr = MatDestroy(&pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 333 pcbddc->user_ChangeOfBasisMatrix = change; 334 pcbddc->change_interior = interior; 335 PetscFunctionReturn(0); 336 } 337 #undef __FUNCT__ 338 #define __FUNCT__ "PCBDDCSetChangeOfBasisMat" 339 /*@ 340 PCBDDCSetChangeOfBasisMat - Set user defined change of basis for dofs 341 342 Collective on PC 343 344 Input Parameters: 345 + pc - the preconditioning context 346 . change - the change of basis matrix 347 - interior - whether or not the change of basis modifies interior dofs 348 349 Level: intermediate 350 351 Notes: 352 353 .seealso: PCBDDC 354 @*/ 355 PetscErrorCode PCBDDCSetChangeOfBasisMat(PC pc, Mat change, PetscBool interior) 356 { 357 PetscErrorCode ierr; 358 359 PetscFunctionBegin; 360 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 361 PetscValidHeaderSpecific(change,MAT_CLASSID,2); 362 PetscCheckSameComm(pc,1,change,2); 363 if (pc->mat) { 364 PetscInt rows_c,cols_c,rows,cols; 365 ierr = MatGetSize(pc->mat,&rows,&cols);CHKERRQ(ierr); 366 ierr = MatGetSize(change,&rows_c,&cols_c);CHKERRQ(ierr); 367 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); 368 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); 369 ierr = MatGetLocalSize(pc->mat,&rows,&cols);CHKERRQ(ierr); 370 ierr = MatGetLocalSize(change,&rows_c,&cols_c);CHKERRQ(ierr); 371 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); 372 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); 373 } 374 ierr = PetscTryMethod(pc,"PCBDDCSetChangeOfBasisMat_C",(PC,Mat,PetscBool),(pc,change,interior));CHKERRQ(ierr); 375 PetscFunctionReturn(0); 376 } 377 378 #undef __FUNCT__ 379 #define __FUNCT__ "PCBDDCSetPrimalVerticesIS_BDDC" 380 static PetscErrorCode PCBDDCSetPrimalVerticesIS_BDDC(PC pc, IS PrimalVertices) 381 { 382 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 383 PetscBool isequal = PETSC_FALSE; 384 PetscErrorCode ierr; 385 386 PetscFunctionBegin; 387 ierr = PetscObjectReference((PetscObject)PrimalVertices);CHKERRQ(ierr); 388 if (pcbddc->user_primal_vertices) { 389 ierr = ISEqual(PrimalVertices,pcbddc->user_primal_vertices,&isequal);CHKERRQ(ierr); 390 } 391 ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr); 392 ierr = ISDestroy(&pcbddc->user_primal_vertices_local);CHKERRQ(ierr); 393 pcbddc->user_primal_vertices = PrimalVertices; 394 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 395 PetscFunctionReturn(0); 396 } 397 #undef __FUNCT__ 398 #define __FUNCT__ "PCBDDCSetPrimalVerticesIS" 399 /*@ 400 PCBDDCSetPrimalVerticesIS - Set additional user defined primal vertices in PCBDDC 401 402 Collective 403 404 Input Parameters: 405 + pc - the preconditioning context 406 - PrimalVertices - index set of primal vertices in global numbering (can be empty) 407 408 Level: intermediate 409 410 Notes: 411 Any process can list any global node 412 413 .seealso: PCBDDC 414 @*/ 415 PetscErrorCode PCBDDCSetPrimalVerticesIS(PC pc, IS PrimalVertices) 416 { 417 PetscErrorCode ierr; 418 419 PetscFunctionBegin; 420 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 421 PetscValidHeaderSpecific(PrimalVertices,IS_CLASSID,2); 422 PetscCheckSameComm(pc,1,PrimalVertices,2); 423 ierr = PetscTryMethod(pc,"PCBDDCSetPrimalVerticesIS_C",(PC,IS),(pc,PrimalVertices));CHKERRQ(ierr); 424 PetscFunctionReturn(0); 425 } 426 427 #undef __FUNCT__ 428 #define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS_BDDC" 429 static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices) 430 { 431 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 432 PetscBool isequal = PETSC_FALSE; 433 PetscErrorCode ierr; 434 435 PetscFunctionBegin; 436 ierr = PetscObjectReference((PetscObject)PrimalVertices);CHKERRQ(ierr); 437 if (pcbddc->user_primal_vertices_local) { 438 ierr = ISEqual(PrimalVertices,pcbddc->user_primal_vertices_local,&isequal);CHKERRQ(ierr); 439 } 440 ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr); 441 ierr = ISDestroy(&pcbddc->user_primal_vertices_local);CHKERRQ(ierr); 442 pcbddc->user_primal_vertices_local = PrimalVertices; 443 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 444 PetscFunctionReturn(0); 445 } 446 #undef __FUNCT__ 447 #define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS" 448 /*@ 449 PCBDDCSetPrimalVerticesLocalIS - Set additional user defined primal vertices in PCBDDC 450 451 Collective 452 453 Input Parameters: 454 + pc - the preconditioning context 455 - PrimalVertices - index set of primal vertices in local numbering (can be empty) 456 457 Level: intermediate 458 459 Notes: 460 461 .seealso: PCBDDC 462 @*/ 463 PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices) 464 { 465 PetscErrorCode ierr; 466 467 PetscFunctionBegin; 468 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 469 PetscValidHeaderSpecific(PrimalVertices,IS_CLASSID,2); 470 PetscCheckSameComm(pc,1,PrimalVertices,2); 471 ierr = PetscTryMethod(pc,"PCBDDCSetPrimalVerticesLocalIS_C",(PC,IS),(pc,PrimalVertices));CHKERRQ(ierr); 472 PetscFunctionReturn(0); 473 } 474 475 #undef __FUNCT__ 476 #define __FUNCT__ "PCBDDCSetCoarseningRatio_BDDC" 477 static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc,PetscInt k) 478 { 479 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 480 481 PetscFunctionBegin; 482 pcbddc->coarsening_ratio = k; 483 PetscFunctionReturn(0); 484 } 485 486 #undef __FUNCT__ 487 #define __FUNCT__ "PCBDDCSetCoarseningRatio" 488 /*@ 489 PCBDDCSetCoarseningRatio - Set coarsening ratio used in multilevel 490 491 Logically collective on PC 492 493 Input Parameters: 494 + pc - the preconditioning context 495 - k - coarsening ratio (H/h at the coarser level) 496 497 Options Database Keys: 498 . -pc_bddc_coarsening_ratio 499 500 Level: intermediate 501 502 Notes: 503 Approximatively k subdomains at the finer level will be aggregated into a single subdomain at the coarser level 504 505 .seealso: PCBDDC, PCBDDCSetLevels() 506 @*/ 507 PetscErrorCode PCBDDCSetCoarseningRatio(PC pc,PetscInt k) 508 { 509 PetscErrorCode ierr; 510 511 PetscFunctionBegin; 512 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 513 PetscValidLogicalCollectiveInt(pc,k,2); 514 ierr = PetscTryMethod(pc,"PCBDDCSetCoarseningRatio_C",(PC,PetscInt),(pc,k));CHKERRQ(ierr); 515 PetscFunctionReturn(0); 516 } 517 518 /* The following functions (PCBDDCSetUseExactDirichlet PCBDDCSetLevel) are not public */ 519 #undef __FUNCT__ 520 #define __FUNCT__ "PCBDDCSetUseExactDirichlet_BDDC" 521 static PetscErrorCode PCBDDCSetUseExactDirichlet_BDDC(PC pc,PetscBool flg) 522 { 523 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 524 525 PetscFunctionBegin; 526 pcbddc->use_exact_dirichlet_trick = flg; 527 PetscFunctionReturn(0); 528 } 529 530 #undef __FUNCT__ 531 #define __FUNCT__ "PCBDDCSetUseExactDirichlet" 532 PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool flg) 533 { 534 PetscErrorCode ierr; 535 536 PetscFunctionBegin; 537 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 538 PetscValidLogicalCollectiveBool(pc,flg,2); 539 ierr = PetscTryMethod(pc,"PCBDDCSetUseExactDirichlet_C",(PC,PetscBool),(pc,flg));CHKERRQ(ierr); 540 PetscFunctionReturn(0); 541 } 542 543 #undef __FUNCT__ 544 #define __FUNCT__ "PCBDDCSetLevel_BDDC" 545 static PetscErrorCode PCBDDCSetLevel_BDDC(PC pc,PetscInt level) 546 { 547 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 548 549 PetscFunctionBegin; 550 pcbddc->current_level = level; 551 PetscFunctionReturn(0); 552 } 553 554 #undef __FUNCT__ 555 #define __FUNCT__ "PCBDDCSetLevel" 556 PetscErrorCode PCBDDCSetLevel(PC pc,PetscInt level) 557 { 558 PetscErrorCode ierr; 559 560 PetscFunctionBegin; 561 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 562 PetscValidLogicalCollectiveInt(pc,level,2); 563 ierr = PetscTryMethod(pc,"PCBDDCSetLevel_C",(PC,PetscInt),(pc,level));CHKERRQ(ierr); 564 PetscFunctionReturn(0); 565 } 566 567 #undef __FUNCT__ 568 #define __FUNCT__ "PCBDDCSetLevels_BDDC" 569 static PetscErrorCode PCBDDCSetLevels_BDDC(PC pc,PetscInt levels) 570 { 571 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 572 573 PetscFunctionBegin; 574 pcbddc->max_levels = levels; 575 PetscFunctionReturn(0); 576 } 577 578 #undef __FUNCT__ 579 #define __FUNCT__ "PCBDDCSetLevels" 580 /*@ 581 PCBDDCSetLevels - Sets the maximum number of levels for multilevel 582 583 Logically collective on PC 584 585 Input Parameters: 586 + pc - the preconditioning context 587 - levels - the maximum number of levels (max 9) 588 589 Options Database Keys: 590 . -pc_bddc_levels 591 592 Level: intermediate 593 594 Notes: 595 Default value is 0, i.e. traditional one-level BDDC 596 597 .seealso: PCBDDC, PCBDDCSetCoarseningRatio() 598 @*/ 599 PetscErrorCode PCBDDCSetLevels(PC pc,PetscInt levels) 600 { 601 PetscErrorCode ierr; 602 603 PetscFunctionBegin; 604 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 605 PetscValidLogicalCollectiveInt(pc,levels,2); 606 if (levels > 99) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Maximum number of levels for bddc is 99\n"); 607 ierr = PetscTryMethod(pc,"PCBDDCSetLevels_C",(PC,PetscInt),(pc,levels));CHKERRQ(ierr); 608 PetscFunctionReturn(0); 609 } 610 611 #undef __FUNCT__ 612 #define __FUNCT__ "PCBDDCSetDirichletBoundaries_BDDC" 613 static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc,IS DirichletBoundaries) 614 { 615 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 616 PetscBool isequal = PETSC_FALSE; 617 PetscErrorCode ierr; 618 619 PetscFunctionBegin; 620 ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr); 621 if (pcbddc->DirichletBoundaries) { 622 ierr = ISEqual(DirichletBoundaries,pcbddc->DirichletBoundaries,&isequal);CHKERRQ(ierr); 623 } 624 /* last user setting takes precendence -> destroy any other customization */ 625 ierr = ISDestroy(&pcbddc->DirichletBoundariesLocal);CHKERRQ(ierr); 626 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 627 pcbddc->DirichletBoundaries = DirichletBoundaries; 628 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 629 PetscFunctionReturn(0); 630 } 631 632 #undef __FUNCT__ 633 #define __FUNCT__ "PCBDDCSetDirichletBoundaries" 634 /*@ 635 PCBDDCSetDirichletBoundaries - Set IS defining Dirichlet boundaries for the global problem. 636 637 Collective 638 639 Input Parameters: 640 + pc - the preconditioning context 641 - DirichletBoundaries - parallel IS defining the Dirichlet boundaries 642 643 Level: intermediate 644 645 Notes: 646 Provide the information if you used MatZeroRows/Columns routines. Any process can list any global node 647 648 .seealso: PCBDDC, PCBDDCSetDirichletBoundariesLocal(), MatZeroRows(), MatZeroRowsColumns() 649 @*/ 650 PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc,IS DirichletBoundaries) 651 { 652 PetscErrorCode ierr; 653 654 PetscFunctionBegin; 655 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 656 PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2); 657 PetscCheckSameComm(pc,1,DirichletBoundaries,2); 658 ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundaries_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr); 659 PetscFunctionReturn(0); 660 } 661 662 #undef __FUNCT__ 663 #define __FUNCT__ "PCBDDCSetDirichletBoundariesLocal_BDDC" 664 static PetscErrorCode PCBDDCSetDirichletBoundariesLocal_BDDC(PC pc,IS DirichletBoundaries) 665 { 666 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 667 PetscBool isequal = PETSC_FALSE; 668 PetscErrorCode ierr; 669 670 PetscFunctionBegin; 671 ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr); 672 if (pcbddc->DirichletBoundariesLocal) { 673 ierr = ISEqual(DirichletBoundaries,pcbddc->DirichletBoundariesLocal,&isequal);CHKERRQ(ierr); 674 } 675 /* last user setting takes precendence -> destroy any other customization */ 676 ierr = ISDestroy(&pcbddc->DirichletBoundariesLocal);CHKERRQ(ierr); 677 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 678 pcbddc->DirichletBoundariesLocal = DirichletBoundaries; 679 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 680 PetscFunctionReturn(0); 681 } 682 683 #undef __FUNCT__ 684 #define __FUNCT__ "PCBDDCSetDirichletBoundariesLocal" 685 /*@ 686 PCBDDCSetDirichletBoundariesLocal - Set IS defining Dirichlet boundaries for the global problem in local ordering. 687 688 Collective 689 690 Input Parameters: 691 + pc - the preconditioning context 692 - DirichletBoundaries - parallel IS defining the Dirichlet boundaries (in local ordering) 693 694 Level: intermediate 695 696 Notes: 697 698 .seealso: PCBDDC, PCBDDCSetDirichletBoundaries(), MatZeroRows(), MatZeroRowsColumns() 699 @*/ 700 PetscErrorCode PCBDDCSetDirichletBoundariesLocal(PC pc,IS DirichletBoundaries) 701 { 702 PetscErrorCode ierr; 703 704 PetscFunctionBegin; 705 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 706 PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2); 707 PetscCheckSameComm(pc,1,DirichletBoundaries,2); 708 ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundariesLocal_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr); 709 PetscFunctionReturn(0); 710 } 711 712 #undef __FUNCT__ 713 #define __FUNCT__ "PCBDDCSetNeumannBoundaries_BDDC" 714 static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc,IS NeumannBoundaries) 715 { 716 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 717 PetscBool isequal = PETSC_FALSE; 718 PetscErrorCode ierr; 719 720 PetscFunctionBegin; 721 ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr); 722 if (pcbddc->NeumannBoundaries) { 723 ierr = ISEqual(NeumannBoundaries,pcbddc->NeumannBoundaries,&isequal);CHKERRQ(ierr); 724 } 725 /* last user setting takes precendence -> destroy any other customization */ 726 ierr = ISDestroy(&pcbddc->NeumannBoundariesLocal);CHKERRQ(ierr); 727 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 728 pcbddc->NeumannBoundaries = NeumannBoundaries; 729 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 730 PetscFunctionReturn(0); 731 } 732 733 #undef __FUNCT__ 734 #define __FUNCT__ "PCBDDCSetNeumannBoundaries" 735 /*@ 736 PCBDDCSetNeumannBoundaries - Set IS defining Neumann boundaries for the global problem. 737 738 Collective 739 740 Input Parameters: 741 + pc - the preconditioning context 742 - NeumannBoundaries - parallel IS defining the Neumann boundaries 743 744 Level: intermediate 745 746 Notes: 747 Any process can list any global node 748 749 .seealso: PCBDDC, PCBDDCSetNeumannBoundariesLocal() 750 @*/ 751 PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc,IS NeumannBoundaries) 752 { 753 PetscErrorCode ierr; 754 755 PetscFunctionBegin; 756 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 757 PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2); 758 PetscCheckSameComm(pc,1,NeumannBoundaries,2); 759 ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundaries_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr); 760 PetscFunctionReturn(0); 761 } 762 763 #undef __FUNCT__ 764 #define __FUNCT__ "PCBDDCSetNeumannBoundariesLocal_BDDC" 765 static PetscErrorCode PCBDDCSetNeumannBoundariesLocal_BDDC(PC pc,IS NeumannBoundaries) 766 { 767 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 768 PetscBool isequal = PETSC_FALSE; 769 PetscErrorCode ierr; 770 771 PetscFunctionBegin; 772 ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr); 773 if (pcbddc->NeumannBoundariesLocal) { 774 ierr = ISEqual(NeumannBoundaries,pcbddc->NeumannBoundariesLocal,&isequal);CHKERRQ(ierr); 775 } 776 /* last user setting takes precendence -> destroy any other customization */ 777 ierr = ISDestroy(&pcbddc->NeumannBoundariesLocal);CHKERRQ(ierr); 778 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 779 pcbddc->NeumannBoundariesLocal = NeumannBoundaries; 780 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 781 PetscFunctionReturn(0); 782 } 783 784 #undef __FUNCT__ 785 #define __FUNCT__ "PCBDDCSetNeumannBoundariesLocal" 786 /*@ 787 PCBDDCSetNeumannBoundariesLocal - Set IS defining Neumann boundaries for the global problem in local ordering. 788 789 Collective 790 791 Input Parameters: 792 + pc - the preconditioning context 793 - NeumannBoundaries - parallel IS defining the subdomain part of Neumann boundaries (in local ordering) 794 795 Level: intermediate 796 797 Notes: 798 799 .seealso: PCBDDC, PCBDDCSetNeumannBoundaries() 800 @*/ 801 PetscErrorCode PCBDDCSetNeumannBoundariesLocal(PC pc,IS NeumannBoundaries) 802 { 803 PetscErrorCode ierr; 804 805 PetscFunctionBegin; 806 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 807 PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2); 808 PetscCheckSameComm(pc,1,NeumannBoundaries,2); 809 ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundariesLocal_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr); 810 PetscFunctionReturn(0); 811 } 812 813 #undef __FUNCT__ 814 #define __FUNCT__ "PCBDDCGetDirichletBoundaries_BDDC" 815 static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc,IS *DirichletBoundaries) 816 { 817 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 818 819 PetscFunctionBegin; 820 *DirichletBoundaries = pcbddc->DirichletBoundaries; 821 PetscFunctionReturn(0); 822 } 823 824 #undef __FUNCT__ 825 #define __FUNCT__ "PCBDDCGetDirichletBoundaries" 826 /*@ 827 PCBDDCGetDirichletBoundaries - Get parallel IS for Dirichlet boundaries 828 829 Collective 830 831 Input Parameters: 832 . pc - the preconditioning context 833 834 Output Parameters: 835 . DirichletBoundaries - index set defining the Dirichlet boundaries 836 837 Level: intermediate 838 839 Notes: 840 The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetDirichletBoundaries 841 842 .seealso: PCBDDC 843 @*/ 844 PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc,IS *DirichletBoundaries) 845 { 846 PetscErrorCode ierr; 847 848 PetscFunctionBegin; 849 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 850 ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundaries_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr); 851 PetscFunctionReturn(0); 852 } 853 854 #undef __FUNCT__ 855 #define __FUNCT__ "PCBDDCGetDirichletBoundariesLocal_BDDC" 856 static PetscErrorCode PCBDDCGetDirichletBoundariesLocal_BDDC(PC pc,IS *DirichletBoundaries) 857 { 858 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 859 860 PetscFunctionBegin; 861 *DirichletBoundaries = pcbddc->DirichletBoundariesLocal; 862 PetscFunctionReturn(0); 863 } 864 865 #undef __FUNCT__ 866 #define __FUNCT__ "PCBDDCGetDirichletBoundariesLocal" 867 /*@ 868 PCBDDCGetDirichletBoundariesLocal - Get parallel IS for Dirichlet boundaries (in local ordering) 869 870 Collective 871 872 Input Parameters: 873 . pc - the preconditioning context 874 875 Output Parameters: 876 . DirichletBoundaries - index set defining the subdomain part of Dirichlet boundaries 877 878 Level: intermediate 879 880 Notes: 881 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). 882 In the latter case, the IS will be available after PCSetUp. 883 884 .seealso: PCBDDC 885 @*/ 886 PetscErrorCode PCBDDCGetDirichletBoundariesLocal(PC pc,IS *DirichletBoundaries) 887 { 888 PetscErrorCode ierr; 889 890 PetscFunctionBegin; 891 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 892 ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundariesLocal_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr); 893 PetscFunctionReturn(0); 894 } 895 896 #undef __FUNCT__ 897 #define __FUNCT__ "PCBDDCGetNeumannBoundaries_BDDC" 898 static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc,IS *NeumannBoundaries) 899 { 900 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 901 902 PetscFunctionBegin; 903 *NeumannBoundaries = pcbddc->NeumannBoundaries; 904 PetscFunctionReturn(0); 905 } 906 907 #undef __FUNCT__ 908 #define __FUNCT__ "PCBDDCGetNeumannBoundaries" 909 /*@ 910 PCBDDCGetNeumannBoundaries - Get parallel IS for Neumann boundaries 911 912 Collective 913 914 Input Parameters: 915 . pc - the preconditioning context 916 917 Output Parameters: 918 . NeumannBoundaries - index set defining the Neumann boundaries 919 920 Level: intermediate 921 922 Notes: 923 The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetNeumannBoundaries 924 925 .seealso: PCBDDC 926 @*/ 927 PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc,IS *NeumannBoundaries) 928 { 929 PetscErrorCode ierr; 930 931 PetscFunctionBegin; 932 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 933 ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundaries_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr); 934 PetscFunctionReturn(0); 935 } 936 937 #undef __FUNCT__ 938 #define __FUNCT__ "PCBDDCGetNeumannBoundariesLocal_BDDC" 939 static PetscErrorCode PCBDDCGetNeumannBoundariesLocal_BDDC(PC pc,IS *NeumannBoundaries) 940 { 941 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 942 943 PetscFunctionBegin; 944 *NeumannBoundaries = pcbddc->NeumannBoundariesLocal; 945 PetscFunctionReturn(0); 946 } 947 948 #undef __FUNCT__ 949 #define __FUNCT__ "PCBDDCGetNeumannBoundariesLocal" 950 /*@ 951 PCBDDCGetNeumannBoundariesLocal - Get parallel IS for Neumann boundaries (in local ordering) 952 953 Collective 954 955 Input Parameters: 956 . pc - the preconditioning context 957 958 Output Parameters: 959 . NeumannBoundaries - index set defining the subdomain part of Neumann boundaries 960 961 Level: intermediate 962 963 Notes: 964 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). 965 In the latter case, the IS will be available after PCSetUp. 966 967 .seealso: PCBDDC 968 @*/ 969 PetscErrorCode PCBDDCGetNeumannBoundariesLocal(PC pc,IS *NeumannBoundaries) 970 { 971 PetscErrorCode ierr; 972 973 PetscFunctionBegin; 974 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 975 ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundariesLocal_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr); 976 PetscFunctionReturn(0); 977 } 978 979 #undef __FUNCT__ 980 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph_BDDC" 981 static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode) 982 { 983 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 984 PCBDDCGraph mat_graph = pcbddc->mat_graph; 985 PetscBool same_data = PETSC_FALSE; 986 PetscErrorCode ierr; 987 988 PetscFunctionBegin; 989 if (!nvtxs) { 990 ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr); 991 PetscFunctionReturn(0); 992 } 993 if (mat_graph->nvtxs_csr == nvtxs && mat_graph->freecsr) { /* we own the data */ 994 if (mat_graph->xadj == xadj && mat_graph->adjncy == adjncy) same_data = PETSC_TRUE; 995 if (!same_data && mat_graph->xadj[nvtxs] == xadj[nvtxs]) { 996 ierr = PetscMemcmp(xadj,mat_graph->xadj,(nvtxs+1)*sizeof(PetscInt),&same_data);CHKERRQ(ierr); 997 if (same_data) { 998 ierr = PetscMemcmp(adjncy,mat_graph->adjncy,xadj[nvtxs]*sizeof(PetscInt),&same_data);CHKERRQ(ierr); 999 } 1000 } 1001 } 1002 if (!same_data) { 1003 /* free old CSR */ 1004 ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr); 1005 /* get CSR into graph structure */ 1006 if (copymode == PETSC_COPY_VALUES) { 1007 ierr = PetscMalloc1(nvtxs+1,&mat_graph->xadj);CHKERRQ(ierr); 1008 ierr = PetscMalloc1(xadj[nvtxs],&mat_graph->adjncy);CHKERRQ(ierr); 1009 ierr = PetscMemcpy(mat_graph->xadj,xadj,(nvtxs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1010 ierr = PetscMemcpy(mat_graph->adjncy,adjncy,xadj[nvtxs]*sizeof(PetscInt));CHKERRQ(ierr); 1011 mat_graph->freecsr = PETSC_TRUE; 1012 } else if (copymode == PETSC_OWN_POINTER) { 1013 mat_graph->xadj = (PetscInt*)xadj; 1014 mat_graph->adjncy = (PetscInt*)adjncy; 1015 mat_graph->freecsr = PETSC_TRUE; 1016 } else if (copymode == PETSC_USE_POINTER) { 1017 mat_graph->xadj = (PetscInt*)xadj; 1018 mat_graph->adjncy = (PetscInt*)adjncy; 1019 mat_graph->freecsr = PETSC_FALSE; 1020 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported copy mode %d",copymode); 1021 mat_graph->nvtxs_csr = nvtxs; 1022 pcbddc->recompute_topography = PETSC_TRUE; 1023 } 1024 PetscFunctionReturn(0); 1025 } 1026 1027 #undef __FUNCT__ 1028 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph" 1029 /*@ 1030 PCBDDCSetLocalAdjacencyGraph - Set adjacency structure (CSR graph) of the local degrees of freedom. 1031 1032 Not collective 1033 1034 Input Parameters: 1035 + pc - the preconditioning context. 1036 . nvtxs - number of local vertices of the graph (i.e., the number of local dofs). 1037 . xadj, adjncy - the connectivity of the dofs in CSR format. 1038 - copymode - supported modes are PETSC_COPY_VALUES, PETSC_USE_POINTER or PETSC_OWN_POINTER. 1039 1040 Level: intermediate 1041 1042 Notes: A dof is considered connected with all local dofs if xadj[dof+1]-xadj[dof] == 1 and adjncy[xadj[dof]] is negative. 1043 1044 .seealso: PCBDDC,PetscCopyMode 1045 @*/ 1046 PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc,PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode) 1047 { 1048 void (*f)(void) = 0; 1049 PetscErrorCode ierr; 1050 1051 PetscFunctionBegin; 1052 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1053 if (nvtxs) { 1054 PetscValidIntPointer(xadj,3); 1055 PetscValidIntPointer(adjncy,4); 1056 } 1057 ierr = PetscTryMethod(pc,"PCBDDCSetLocalAdjacencyGraph_C",(PC,PetscInt,const PetscInt[],const PetscInt[],PetscCopyMode),(pc,nvtxs,xadj,adjncy,copymode));CHKERRQ(ierr); 1058 /* free arrays if PCBDDC is not the PC type */ 1059 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",&f);CHKERRQ(ierr); 1060 if (!f && copymode == PETSC_OWN_POINTER) { 1061 ierr = PetscFree(xadj);CHKERRQ(ierr); 1062 ierr = PetscFree(adjncy);CHKERRQ(ierr); 1063 } 1064 PetscFunctionReturn(0); 1065 } 1066 1067 #undef __FUNCT__ 1068 #define __FUNCT__ "PCBDDCSetDofsSplittingLocal_BDDC" 1069 static PetscErrorCode PCBDDCSetDofsSplittingLocal_BDDC(PC pc,PetscInt n_is, IS ISForDofs[]) 1070 { 1071 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1072 PetscInt i; 1073 PetscBool isequal = PETSC_FALSE; 1074 PetscErrorCode ierr; 1075 1076 PetscFunctionBegin; 1077 if (pcbddc->n_ISForDofsLocal == n_is) { 1078 for (i=0;i<n_is;i++) { 1079 PetscBool isequalt; 1080 ierr = ISEqual(ISForDofs[i],pcbddc->ISForDofsLocal[i],&isequalt);CHKERRQ(ierr); 1081 if (!isequalt) break; 1082 } 1083 if (i == n_is) isequal = PETSC_TRUE; 1084 } 1085 for (i=0;i<n_is;i++) { 1086 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 1087 } 1088 /* Destroy ISes if they were already set */ 1089 for (i=0;i<pcbddc->n_ISForDofsLocal;i++) { 1090 ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr); 1091 } 1092 ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr); 1093 /* last user setting takes precendence -> destroy any other customization */ 1094 for (i=0;i<pcbddc->n_ISForDofs;i++) { 1095 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 1096 } 1097 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 1098 pcbddc->n_ISForDofs = 0; 1099 /* allocate space then set */ 1100 if (n_is) { 1101 ierr = PetscMalloc1(n_is,&pcbddc->ISForDofsLocal);CHKERRQ(ierr); 1102 } 1103 for (i=0;i<n_is;i++) { 1104 pcbddc->ISForDofsLocal[i] = ISForDofs[i]; 1105 } 1106 pcbddc->n_ISForDofsLocal = n_is; 1107 if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE; 1108 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 1109 PetscFunctionReturn(0); 1110 } 1111 1112 #undef __FUNCT__ 1113 #define __FUNCT__ "PCBDDCSetDofsSplittingLocal" 1114 /*@ 1115 PCBDDCSetDofsSplittingLocal - Set index sets defining fields of the local subdomain matrix 1116 1117 Collective 1118 1119 Input Parameters: 1120 + pc - the preconditioning context 1121 . n_is - number of index sets defining the fields 1122 - ISForDofs - array of IS describing the fields in local ordering 1123 1124 Level: intermediate 1125 1126 Notes: 1127 n_is should be the same among processes. Not all nodes need to be listed: unlisted nodes will belong to the complement field. 1128 1129 .seealso: PCBDDC 1130 @*/ 1131 PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc,PetscInt n_is, IS ISForDofs[]) 1132 { 1133 PetscInt i; 1134 PetscErrorCode ierr; 1135 1136 PetscFunctionBegin; 1137 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1138 PetscValidLogicalCollectiveInt(pc,n_is,2); 1139 for (i=0;i<n_is;i++) { 1140 PetscCheckSameComm(pc,1,ISForDofs[i],3); 1141 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 1142 } 1143 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplittingLocal_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 1144 PetscFunctionReturn(0); 1145 } 1146 1147 #undef __FUNCT__ 1148 #define __FUNCT__ "PCBDDCSetDofsSplitting_BDDC" 1149 static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc,PetscInt n_is, IS ISForDofs[]) 1150 { 1151 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1152 PetscInt i; 1153 PetscBool isequal = PETSC_FALSE; 1154 PetscErrorCode ierr; 1155 1156 PetscFunctionBegin; 1157 if (pcbddc->n_ISForDofs == n_is) { 1158 for (i=0;i<n_is;i++) { 1159 PetscBool isequalt; 1160 ierr = ISEqual(ISForDofs[i],pcbddc->ISForDofs[i],&isequalt);CHKERRQ(ierr); 1161 if (!isequalt) break; 1162 } 1163 if (i == n_is) isequal = PETSC_TRUE; 1164 } 1165 for (i=0;i<n_is;i++) { 1166 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 1167 } 1168 /* Destroy ISes if they were already set */ 1169 for (i=0;i<pcbddc->n_ISForDofs;i++) { 1170 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 1171 } 1172 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 1173 /* last user setting takes precendence -> destroy any other customization */ 1174 for (i=0;i<pcbddc->n_ISForDofsLocal;i++) { 1175 ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr); 1176 } 1177 ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr); 1178 pcbddc->n_ISForDofsLocal = 0; 1179 /* allocate space then set */ 1180 if (n_is) { 1181 ierr = PetscMalloc1(n_is,&pcbddc->ISForDofs);CHKERRQ(ierr); 1182 } 1183 for (i=0;i<n_is;i++) { 1184 pcbddc->ISForDofs[i] = ISForDofs[i]; 1185 } 1186 pcbddc->n_ISForDofs = n_is; 1187 if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE; 1188 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 1189 PetscFunctionReturn(0); 1190 } 1191 1192 #undef __FUNCT__ 1193 #define __FUNCT__ "PCBDDCSetDofsSplitting" 1194 /*@ 1195 PCBDDCSetDofsSplitting - Set index sets defining fields of the global matrix 1196 1197 Collective 1198 1199 Input Parameters: 1200 + pc - the preconditioning context 1201 . n_is - number of index sets defining the fields 1202 - ISForDofs - array of IS describing the fields in global ordering 1203 1204 Level: intermediate 1205 1206 Notes: 1207 Any process can list any global node. Not all nodes need to be listed: unlisted nodes will belong to the complement field. 1208 1209 .seealso: PCBDDC 1210 @*/ 1211 PetscErrorCode PCBDDCSetDofsSplitting(PC pc,PetscInt n_is, IS ISForDofs[]) 1212 { 1213 PetscInt i; 1214 PetscErrorCode ierr; 1215 1216 PetscFunctionBegin; 1217 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1218 PetscValidLogicalCollectiveInt(pc,n_is,2); 1219 for (i=0;i<n_is;i++) { 1220 PetscCheckSameComm(pc,1,ISForDofs[i],3); 1221 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 1222 } 1223 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 1224 PetscFunctionReturn(0); 1225 } 1226 1227 #undef __FUNCT__ 1228 #define __FUNCT__ "PCPreSolve_BDDC" 1229 /* 1230 PCPreSolve_BDDC - Changes the right hand side and (if necessary) the initial 1231 guess if a transformation of basis approach has been selected. 1232 1233 Input Parameter: 1234 + pc - the preconditioner contex 1235 1236 Application Interface Routine: PCPreSolve() 1237 1238 Notes: 1239 The interface routine PCPreSolve() is not usually called directly by 1240 the user, but instead is called by KSPSolve(). 1241 */ 1242 static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 1243 { 1244 PetscErrorCode ierr; 1245 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1246 PC_IS *pcis = (PC_IS*)(pc->data); 1247 Vec used_vec; 1248 PetscBool save_rhs = PETSC_TRUE, benign_correction_computed; 1249 1250 PetscFunctionBegin; 1251 /* if we are working with CG, one dirichlet solve can be avoided during Krylov iterations */ 1252 if (ksp) { 1253 PetscBool iscg, isgroppcg, ispipecg, ispipecgrr; 1254 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPCG,&iscg);CHKERRQ(ierr); 1255 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPGROPPCG,&isgroppcg);CHKERRQ(ierr); 1256 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPIPECG,&ispipecg);CHKERRQ(ierr); 1257 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPIPECGRR,&ispipecgrr);CHKERRQ(ierr); 1258 if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || (!iscg && !isgroppcg && !ispipecg && !ispipecgrr)) { 1259 ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr); 1260 } 1261 } 1262 if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static) { 1263 ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr); 1264 } 1265 1266 /* Creates parallel work vectors used in presolve */ 1267 if (!pcbddc->original_rhs) { 1268 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 1269 } 1270 if (!pcbddc->temp_solution) { 1271 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);CHKERRQ(ierr); 1272 } 1273 1274 pcbddc->temp_solution_used = PETSC_FALSE; 1275 if (x) { 1276 ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr); 1277 used_vec = x; 1278 } else { /* it can only happen when calling PCBDDCMatFETIDPGetRHS */ 1279 ierr = PetscObjectReference((PetscObject)pcbddc->temp_solution);CHKERRQ(ierr); 1280 used_vec = pcbddc->temp_solution; 1281 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 1282 pcbddc->temp_solution_used = PETSC_TRUE; 1283 ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1284 save_rhs = PETSC_FALSE; 1285 pcbddc->eliminate_dirdofs = PETSC_TRUE; 1286 } 1287 1288 /* hack into ksp data structure since PCPreSolve comes earlier than setting to zero the guess in src/ksp/ksp/interface/itfunc.c */ 1289 if (ksp) { 1290 /* store the flag for the initial guess since it will be restored back during PCPostSolve_BDDC */ 1291 ierr = KSPGetInitialGuessNonzero(ksp,&pcbddc->ksp_guess_nonzero);CHKERRQ(ierr); 1292 if (!pcbddc->ksp_guess_nonzero) { 1293 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 1294 } 1295 } 1296 1297 pcbddc->rhs_change = PETSC_FALSE; 1298 /* Take into account zeroed rows -> change rhs and store solution removed */ 1299 if (rhs && pcbddc->eliminate_dirdofs) { 1300 IS dirIS = NULL; 1301 1302 /* DirichletBoundariesLocal may not be consistent among neighbours; gets a dirichlet dofs IS from graph (may be cached) */ 1303 ierr = PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph,&dirIS);CHKERRQ(ierr); 1304 if (dirIS) { 1305 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 1306 PetscInt dirsize,i,*is_indices; 1307 PetscScalar *array_x; 1308 const PetscScalar *array_diagonal; 1309 1310 ierr = MatGetDiagonal(pc->pmat,pcis->vec1_global);CHKERRQ(ierr); 1311 ierr = VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);CHKERRQ(ierr); 1312 ierr = VecScatterBegin(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1313 ierr = VecScatterEnd(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1314 ierr = VecScatterBegin(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1315 ierr = VecScatterEnd(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1316 ierr = ISGetLocalSize(dirIS,&dirsize);CHKERRQ(ierr); 1317 ierr = VecGetArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 1318 ierr = VecGetArrayRead(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 1319 ierr = ISGetIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1320 for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]]; 1321 ierr = ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1322 ierr = VecRestoreArrayRead(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 1323 ierr = VecRestoreArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 1324 ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1325 ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1326 pcbddc->rhs_change = PETSC_TRUE; 1327 ierr = ISDestroy(&dirIS);CHKERRQ(ierr); 1328 } 1329 } 1330 1331 /* remove the computed solution or the initial guess from the rhs */ 1332 if (pcbddc->rhs_change || (ksp && pcbddc->ksp_guess_nonzero) ) { 1333 /* save the original rhs */ 1334 if (save_rhs) { 1335 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1336 save_rhs = PETSC_FALSE; 1337 } 1338 pcbddc->rhs_change = PETSC_TRUE; 1339 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 1340 ierr = MatMultAdd(pc->mat,used_vec,pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1341 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 1342 ierr = VecCopy(used_vec,pcbddc->temp_solution);CHKERRQ(ierr); 1343 pcbddc->temp_solution_used = PETSC_TRUE; 1344 if (ksp) { 1345 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_FALSE);CHKERRQ(ierr); 1346 } 1347 } 1348 ierr = VecDestroy(&used_vec);CHKERRQ(ierr); 1349 1350 /* compute initial vector in benign space if needed 1351 and remove non-benign solution from the rhs */ 1352 benign_correction_computed = PETSC_FALSE; 1353 if (rhs && pcbddc->benign_compute_correction && pcbddc->benign_have_null) { 1354 /* compute u^*_h using ideas similar to those in Xuemin Tu's PhD thesis (see Section 4.8.1) 1355 Recursively apply BDDC in the multilevel case */ 1356 if (!pcbddc->benign_vec) { 1357 ierr = VecDuplicate(rhs,&pcbddc->benign_vec);CHKERRQ(ierr); 1358 } 1359 pcbddc->benign_apply_coarse_only = PETSC_TRUE; 1360 if (!pcbddc->benign_skip_correction) { 1361 ierr = PCApply_BDDC(pc,rhs,pcbddc->benign_vec);CHKERRQ(ierr); 1362 benign_correction_computed = PETSC_TRUE; 1363 if (pcbddc->temp_solution_used) { 1364 ierr = VecAXPY(pcbddc->temp_solution,1.0,pcbddc->benign_vec);CHKERRQ(ierr); 1365 } 1366 ierr = VecScale(pcbddc->benign_vec,-1.0);CHKERRQ(ierr); 1367 /* store the original rhs if not done earlier */ 1368 if (save_rhs) { 1369 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 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 pcbddc->work_change = r; 1918 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 1919 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 1920 } 1921 PetscFunctionReturn(0); 1922 } 1923 1924 /* 1925 PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector. 1926 1927 Input Parameters: 1928 + pc - the preconditioner context 1929 - r - input vector (global) 1930 1931 Output Parameter: 1932 . z - output vector (global) 1933 1934 Application Interface Routine: PCApplyTranspose() 1935 */ 1936 #undef __FUNCT__ 1937 #define __FUNCT__ "PCApplyTranspose_BDDC" 1938 PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z) 1939 { 1940 PC_IS *pcis = (PC_IS*)(pc->data); 1941 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1942 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1943 PetscErrorCode ierr; 1944 const PetscScalar one = 1.0; 1945 const PetscScalar m_one = -1.0; 1946 const PetscScalar zero = 0.0; 1947 1948 PetscFunctionBegin; 1949 if (pcbddc->ChangeOfBasisMatrix) { 1950 Vec swap; 1951 1952 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);CHKERRQ(ierr); 1953 swap = pcbddc->work_change; 1954 pcbddc->work_change = r; 1955 r = swap; 1956 /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */ 1957 if (pcbddc->benign_apply_coarse_only && pcbddc->exact_dirichlet_trick_app && pcbddc->change_interior) { 1958 ierr = VecCopy(r,pcis->vec1_global);CHKERRQ(ierr); 1959 ierr = VecLockPush(pcis->vec1_global);CHKERRQ(ierr); 1960 } 1961 } 1962 if (pcbddc->benign_have_null) { /* get p0 from r */ 1963 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 1964 } 1965 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 1966 ierr = VecCopy(r,z);CHKERRQ(ierr); 1967 /* First Dirichlet solve */ 1968 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1969 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1970 /* 1971 Assembling right hand side for BDDC operator 1972 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1973 - pcis->vec1_B the interface part of the global vector z 1974 */ 1975 if (n_D) { 1976 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1977 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1978 if (pcbddc->switch_static) { 1979 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1980 1981 ierr = VecSet(pcis->vec1_N,0.);CHKERRQ(ierr); 1982 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1983 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1984 if (!pcbddc->switch_static_change) { 1985 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1986 } else { 1987 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1988 ierr = MatMultTranspose(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1989 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1990 } 1991 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1992 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1993 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1994 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1995 } else { 1996 ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1997 } 1998 } else { 1999 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 2000 } 2001 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2002 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2003 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 2004 } else { 2005 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 2006 } 2007 2008 /* Apply interface preconditioner 2009 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 2010 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr); 2011 2012 /* Apply transpose of partition of unity operator */ 2013 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 2014 2015 /* Second Dirichlet solve and assembling of output */ 2016 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2017 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2018 if (n_B) { 2019 if (pcbddc->switch_static) { 2020 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 2021 2022 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2023 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2024 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2025 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2026 if (!pcbddc->switch_static_change) { 2027 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2028 } else { 2029 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2030 ierr = MatMultTranspose(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 2031 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2032 } 2033 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2034 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2035 } else { 2036 ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 2037 } 2038 } else if (pcbddc->switch_static) { /* n_B is zero */ 2039 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 2040 2041 if (!pcbddc->switch_static_change) { 2042 ierr = MatMultTranspose(matis->A,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 2043 } else { 2044 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);CHKERRQ(ierr); 2045 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2046 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);CHKERRQ(ierr); 2047 } 2048 } 2049 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 2050 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 2051 if (pcbddc->switch_static) { 2052 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 2053 } else { 2054 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 2055 } 2056 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2057 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2058 } else { 2059 if (pcbddc->switch_static) { 2060 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 2061 } else { 2062 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 2063 } 2064 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2065 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2066 } 2067 if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */ 2068 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 2069 } 2070 if (pcbddc->ChangeOfBasisMatrix) { 2071 pcbddc->work_change = r; 2072 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 2073 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 2074 } 2075 PetscFunctionReturn(0); 2076 } 2077 2078 #undef __FUNCT__ 2079 #define __FUNCT__ "PCDestroy_BDDC" 2080 PetscErrorCode PCDestroy_BDDC(PC pc) 2081 { 2082 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 2083 PetscErrorCode ierr; 2084 2085 PetscFunctionBegin; 2086 /* free BDDC custom data */ 2087 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 2088 /* destroy objects related to topography */ 2089 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 2090 /* free allocated graph structure */ 2091 ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr); 2092 /* free allocated sub schurs structure */ 2093 ierr = PetscFree(pcbddc->sub_schurs);CHKERRQ(ierr); 2094 /* destroy objects for scaling operator */ 2095 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 2096 ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr); 2097 /* free solvers stuff */ 2098 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 2099 /* free global vectors needed in presolve */ 2100 ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr); 2101 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 2102 /* free data created by PCIS */ 2103 ierr = PCISDestroy(pc);CHKERRQ(ierr); 2104 /* remove functions */ 2105 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDiscreteGradient_C",NULL);CHKERRQ(ierr); 2106 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",NULL);CHKERRQ(ierr); 2107 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);CHKERRQ(ierr); 2108 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 2109 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",NULL);CHKERRQ(ierr); 2110 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 2111 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr); 2112 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr); 2113 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr); 2114 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2115 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2116 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2117 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2118 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2119 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2120 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2121 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2122 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 2123 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr); 2124 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 2125 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 2126 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 2127 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 2128 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",NULL);CHKERRQ(ierr); 2129 /* Free the private data structure */ 2130 ierr = PetscFree(pc->data);CHKERRQ(ierr); 2131 PetscFunctionReturn(0); 2132 } 2133 2134 #undef __FUNCT__ 2135 #define __FUNCT__ "PCPreSolveChangeRHS_BDDC" 2136 static PetscErrorCode PCPreSolveChangeRHS_BDDC(PC pc, PetscBool* change) 2137 { 2138 PetscFunctionBegin; 2139 *change = PETSC_TRUE; 2140 PetscFunctionReturn(0); 2141 } 2142 2143 #undef __FUNCT__ 2144 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC" 2145 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 2146 { 2147 FETIDPMat_ctx mat_ctx; 2148 Vec work; 2149 PC_IS* pcis; 2150 PC_BDDC* pcbddc; 2151 PetscErrorCode ierr; 2152 2153 PetscFunctionBegin; 2154 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2155 pcis = (PC_IS*)mat_ctx->pc->data; 2156 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2157 2158 /* 2159 change of basis for physical rhs if needed 2160 It also changes the rhs in case of dirichlet boundaries 2161 */ 2162 if (!pcbddc->original_rhs) { 2163 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 2164 } 2165 ierr = VecCopy(standard_rhs,pcbddc->original_rhs);CHKERRQ(ierr); 2166 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,pcbddc->original_rhs,NULL);CHKERRQ(ierr); 2167 if (pcbddc->ChangeOfBasisMatrix) { 2168 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,pcbddc->original_rhs,pcbddc->work_change);CHKERRQ(ierr); 2169 work = pcbddc->work_change; 2170 } else { 2171 work = pcbddc->original_rhs; 2172 } 2173 /* store vectors for computation of fetidp final solution */ 2174 ierr = VecScatterBegin(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2175 ierr = VecScatterEnd(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2176 /* scale rhs since it should be unassembled */ 2177 /* TODO use counter scaling? (also below) */ 2178 ierr = VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2179 ierr = VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2180 /* Apply partition of unity */ 2181 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2182 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2183 if (!pcbddc->switch_static) { 2184 /* compute partially subassembled Schur complement right-hand side */ 2185 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2186 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 2187 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 2188 ierr = VecSet(work,0.0);CHKERRQ(ierr); 2189 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2190 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2191 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2192 ierr = VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2193 ierr = VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2194 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2195 } 2196 /* BDDC rhs */ 2197 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 2198 if (pcbddc->switch_static) { 2199 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2200 } 2201 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2202 /* apply BDDC */ 2203 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2204 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2205 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 2206 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 2207 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 2208 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2209 ierr = VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2210 PetscFunctionReturn(0); 2211 } 2212 2213 #undef __FUNCT__ 2214 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS" 2215 /*@ 2216 PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETI-DP linear system using the physical right-hand side 2217 2218 Collective 2219 2220 Input Parameters: 2221 + fetidp_mat - the FETI-DP matrix object obtained by a call to PCBDDCCreateFETIDPOperators 2222 - standard_rhs - the right-hand side of the original linear system 2223 2224 Output Parameters: 2225 . fetidp_flux_rhs - the right-hand side for the FETI-DP linear system 2226 2227 Level: developer 2228 2229 Notes: 2230 2231 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetSolution 2232 @*/ 2233 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 2234 { 2235 FETIDPMat_ctx mat_ctx; 2236 PetscErrorCode ierr; 2237 2238 PetscFunctionBegin; 2239 PetscValidHeaderSpecific(fetidp_mat,MAT_CLASSID,1); 2240 PetscValidHeaderSpecific(standard_rhs,VEC_CLASSID,2); 2241 PetscValidHeaderSpecific(fetidp_flux_rhs,VEC_CLASSID,3); 2242 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2243 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 2244 PetscFunctionReturn(0); 2245 } 2246 2247 #undef __FUNCT__ 2248 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC" 2249 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2250 { 2251 FETIDPMat_ctx mat_ctx; 2252 PC_IS* pcis; 2253 PC_BDDC* pcbddc; 2254 PetscErrorCode ierr; 2255 2256 PetscFunctionBegin; 2257 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2258 pcis = (PC_IS*)mat_ctx->pc->data; 2259 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2260 2261 /* apply B_delta^T */ 2262 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2263 ierr = VecScatterEnd(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2264 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 2265 /* compute rhs for BDDC application */ 2266 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2267 if (pcbddc->switch_static) { 2268 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2269 } 2270 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2271 /* apply BDDC */ 2272 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2273 /* put values into standard global vector */ 2274 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2275 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2276 if (!pcbddc->switch_static) { 2277 /* compute values into the interior if solved for the partially subassembled Schur complement */ 2278 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 2279 ierr = VecAYPX(pcis->vec1_D,-1.0,mat_ctx->temp_solution_D);CHKERRQ(ierr); 2280 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); 2281 } 2282 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2283 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2284 /* add p0 solution to final solution */ 2285 ierr = PCBDDCBenignGetOrSetP0(mat_ctx->pc,standard_sol,PETSC_FALSE);CHKERRQ(ierr); 2286 if (pcbddc->ChangeOfBasisMatrix) { 2287 Vec v2; 2288 ierr = VecDuplicate(standard_sol,&v2);CHKERRQ(ierr); 2289 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,standard_sol,v2);CHKERRQ(ierr); 2290 ierr = VecCopy(v2,standard_sol);CHKERRQ(ierr); 2291 ierr = VecDestroy(&v2);CHKERRQ(ierr); 2292 } 2293 ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr); 2294 PetscFunctionReturn(0); 2295 } 2296 2297 #undef __FUNCT__ 2298 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution" 2299 /*@ 2300 PCBDDCMatFETIDPGetSolution - Compute the physical solution using the solution of the FETI-DP linear system 2301 2302 Collective 2303 2304 Input Parameters: 2305 + fetidp_mat - the FETI-DP matrix obtained by a call to PCBDDCCreateFETIDPOperators 2306 - fetidp_flux_sol - the solution of the FETI-DP linear system 2307 2308 Output Parameters: 2309 . standard_sol - the solution defined on the physical domain 2310 2311 Level: developer 2312 2313 Notes: 2314 2315 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetRHS 2316 @*/ 2317 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2318 { 2319 FETIDPMat_ctx mat_ctx; 2320 PetscErrorCode ierr; 2321 2322 PetscFunctionBegin; 2323 PetscValidHeaderSpecific(fetidp_mat,MAT_CLASSID,1); 2324 PetscValidHeaderSpecific(fetidp_flux_sol,VEC_CLASSID,2); 2325 PetscValidHeaderSpecific(standard_sol,VEC_CLASSID,3); 2326 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2327 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 2328 PetscFunctionReturn(0); 2329 } 2330 2331 extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec); 2332 extern PetscErrorCode FETIDPMatMultTranspose(Mat,Vec,Vec); 2333 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat); 2334 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec); 2335 extern PetscErrorCode FETIDPPCApplyTranspose(PC,Vec,Vec); 2336 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC); 2337 2338 #undef __FUNCT__ 2339 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC" 2340 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, PetscBool fully_redundant, Mat *fetidp_mat, PC *fetidp_pc) 2341 { 2342 2343 FETIDPMat_ctx fetidpmat_ctx; 2344 Mat newmat; 2345 FETIDPPC_ctx fetidppc_ctx; 2346 PC newpc; 2347 MPI_Comm comm; 2348 PetscErrorCode ierr; 2349 2350 PetscFunctionBegin; 2351 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 2352 /* FETIDP linear matrix */ 2353 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 2354 fetidpmat_ctx->fully_redundant = fully_redundant; 2355 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 2356 ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 2357 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 2358 ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr); 2359 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 2360 ierr = MatSetUp(newmat);CHKERRQ(ierr); 2361 /* FETIDP preconditioner */ 2362 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 2363 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 2364 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 2365 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 2366 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 2367 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 2368 ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 2369 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 2370 ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr); 2371 ierr = PCSetUp(newpc);CHKERRQ(ierr); 2372 /* return pointers for objects created */ 2373 *fetidp_mat=newmat; 2374 *fetidp_pc=newpc; 2375 PetscFunctionReturn(0); 2376 } 2377 2378 #undef __FUNCT__ 2379 #define __FUNCT__ "PCBDDCCreateFETIDPOperators" 2380 /*@ 2381 PCBDDCCreateFETIDPOperators - Create FETI-DP operators 2382 2383 Collective 2384 2385 Input Parameters: 2386 + pc - the BDDC preconditioning context (setup should have been called before) 2387 - fully_redundant - true for a fully redundant set of Lagrange multipliers 2388 2389 Output Parameters: 2390 + fetidp_mat - shell FETI-DP matrix object 2391 - fetidp_pc - shell Dirichlet preconditioner for FETI-DP matrix 2392 2393 Level: developer 2394 2395 Notes: 2396 Currently the only operations provided for FETI-DP matrix are MatMult and MatMultTranspose 2397 2398 .seealso: PCBDDC, PCBDDCMatFETIDPGetRHS, PCBDDCMatFETIDPGetSolution 2399 @*/ 2400 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, PetscBool fully_redundant, Mat *fetidp_mat, PC *fetidp_pc) 2401 { 2402 PetscErrorCode ierr; 2403 2404 PetscFunctionBegin; 2405 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2406 if (pc->setupcalled) { 2407 ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,PetscBool,Mat*,PC*),(pc,fully_redundant,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 2408 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n"); 2409 PetscFunctionReturn(0); 2410 } 2411 /* -------------------------------------------------------------------------- */ 2412 /*MC 2413 PCBDDC - Balancing Domain Decomposition by Constraints. 2414 2415 An implementation of the BDDC preconditioner based on 2416 2417 .vb 2418 [1] C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 2419 [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 2420 [3] J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", http://arxiv.org/abs/0712.3977 2421 [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 2422 .ve 2423 2424 The matrix to be preconditioned (Pmat) must be of type MATIS. 2425 2426 Currently works with MATIS matrices with local matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers. 2427 2428 It also works with unsymmetric and indefinite problems. 2429 2430 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. 2431 2432 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). 2433 2434 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() 2435 Additional information on dofs can be provided by using PCBDDCSetDofsSplitting(), PCBDDCSetDirichletBoundaries(), PCBDDCSetNeumannBoundaries(), and PCBDDCSetPrimalVerticesIS() and their local counterparts. 2436 2437 Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace(). Non-singular modes are retained via SVD. 2438 2439 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. 2440 User defined change of basis can be passed to PCBDDC by using PCBDDCSetChangeOfBasisMat() 2441 2442 The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using a MatPartitioning object. 2443 2444 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. 2445 2446 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. 2447 Deluxe scaling is not supported yet for FETI-DP. 2448 2449 Options Database Keys (some of them, run with -h for a complete list): 2450 2451 . -pc_bddc_use_vertices <true> - use or not vertices in primal space 2452 . -pc_bddc_use_edges <true> - use or not edges in primal space 2453 . -pc_bddc_use_faces <false> - use or not faces in primal space 2454 . -pc_bddc_symmetric <true> - symmetric computation of primal basis functions. Specify false for unsymmetric problems 2455 . -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only) 2456 . -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested 2457 . -pc_bddc_switch_static <false> - switches from M_2 (default) to M_3 operator (see reference article [1]) 2458 . -pc_bddc_levels <0> - maximum number of levels for multilevel 2459 . -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) 2460 . -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) 2461 . -pc_bddc_use_deluxe_scaling <false> - use deluxe scaling 2462 . -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) 2463 . -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) 2464 - -pc_bddc_check_level <0> - set verbosity level of debugging output 2465 2466 Options for Dirichlet, Neumann or coarse solver can be set with 2467 .vb 2468 -pc_bddc_dirichlet_ 2469 -pc_bddc_neumann_ 2470 -pc_bddc_coarse_ 2471 .ve 2472 e.g -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. PCBDDC uses by default KPSPREONLY and PCLU. 2473 2474 When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified as 2475 .vb 2476 -pc_bddc_dirichlet_lN_ 2477 -pc_bddc_neumann_lN_ 2478 -pc_bddc_coarse_lN_ 2479 .ve 2480 Note that level number ranges from the finest (0) to the coarsest (N). 2481 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. 2482 .vb 2483 -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3 2484 .ve 2485 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 2486 2487 Level: intermediate 2488 2489 Developer notes: 2490 2491 Contributed by Stefano Zampini 2492 2493 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 2494 M*/ 2495 2496 #undef __FUNCT__ 2497 #define __FUNCT__ "PCCreate_BDDC" 2498 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 2499 { 2500 PetscErrorCode ierr; 2501 PC_BDDC *pcbddc; 2502 2503 PetscFunctionBegin; 2504 /* Creates the private data structure for this preconditioner and attach it to the PC object. */ 2505 ierr = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr); 2506 pc->data = (void*)pcbddc; 2507 2508 /* create PCIS data structure */ 2509 ierr = PCISCreate(pc);CHKERRQ(ierr); 2510 2511 /* BDDC customization */ 2512 pcbddc->use_local_adj = PETSC_TRUE; 2513 pcbddc->use_vertices = PETSC_TRUE; 2514 pcbddc->use_edges = PETSC_TRUE; 2515 pcbddc->use_faces = PETSC_FALSE; 2516 pcbddc->use_change_of_basis = PETSC_FALSE; 2517 pcbddc->use_change_on_faces = PETSC_FALSE; 2518 pcbddc->switch_static = PETSC_FALSE; 2519 pcbddc->use_nnsp_true = PETSC_FALSE; 2520 pcbddc->use_qr_single = PETSC_FALSE; 2521 pcbddc->symmetric_primal = PETSC_TRUE; 2522 pcbddc->benign_saddle_point = PETSC_FALSE; 2523 pcbddc->benign_have_null = PETSC_FALSE; 2524 pcbddc->vertex_size = 1; 2525 pcbddc->dbg_flag = 0; 2526 /* private */ 2527 pcbddc->local_primal_size = 0; 2528 pcbddc->local_primal_size_cc = 0; 2529 pcbddc->local_primal_ref_node = 0; 2530 pcbddc->local_primal_ref_mult = 0; 2531 pcbddc->n_vertices = 0; 2532 pcbddc->primal_indices_local_idxs = 0; 2533 pcbddc->recompute_topography = PETSC_TRUE; 2534 pcbddc->coarse_size = -1; 2535 pcbddc->new_primal_space = PETSC_FALSE; 2536 pcbddc->new_primal_space_local = PETSC_FALSE; 2537 pcbddc->global_primal_indices = 0; 2538 pcbddc->onearnullspace = 0; 2539 pcbddc->onearnullvecs_state = 0; 2540 pcbddc->user_primal_vertices = 0; 2541 pcbddc->user_primal_vertices_local = 0; 2542 pcbddc->temp_solution = 0; 2543 pcbddc->original_rhs = 0; 2544 pcbddc->local_mat = 0; 2545 pcbddc->ChangeOfBasisMatrix = 0; 2546 pcbddc->user_ChangeOfBasisMatrix = 0; 2547 pcbddc->coarse_vec = 0; 2548 pcbddc->coarse_ksp = 0; 2549 pcbddc->coarse_phi_B = 0; 2550 pcbddc->coarse_phi_D = 0; 2551 pcbddc->coarse_psi_B = 0; 2552 pcbddc->coarse_psi_D = 0; 2553 pcbddc->vec1_P = 0; 2554 pcbddc->vec1_R = 0; 2555 pcbddc->vec2_R = 0; 2556 pcbddc->local_auxmat1 = 0; 2557 pcbddc->local_auxmat2 = 0; 2558 pcbddc->R_to_B = 0; 2559 pcbddc->R_to_D = 0; 2560 pcbddc->ksp_D = 0; 2561 pcbddc->ksp_R = 0; 2562 pcbddc->NeumannBoundaries = 0; 2563 pcbddc->NeumannBoundariesLocal = 0; 2564 pcbddc->DirichletBoundaries = 0; 2565 pcbddc->DirichletBoundariesLocal = 0; 2566 pcbddc->user_provided_isfordofs = PETSC_FALSE; 2567 pcbddc->n_ISForDofs = 0; 2568 pcbddc->n_ISForDofsLocal = 0; 2569 pcbddc->ISForDofs = 0; 2570 pcbddc->ISForDofsLocal = 0; 2571 pcbddc->ConstraintMatrix = 0; 2572 pcbddc->use_exact_dirichlet_trick = PETSC_TRUE; 2573 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 2574 pcbddc->coarse_loc_to_glob = 0; 2575 pcbddc->coarsening_ratio = 8; 2576 pcbddc->coarse_adj_red = 0; 2577 pcbddc->current_level = 0; 2578 pcbddc->max_levels = 0; 2579 pcbddc->use_coarse_estimates = PETSC_FALSE; 2580 pcbddc->coarse_eqs_per_proc = 1; 2581 pcbddc->coarse_subassembling = 0; 2582 pcbddc->detect_disconnected = PETSC_FALSE; 2583 pcbddc->n_local_subs = 0; 2584 pcbddc->local_subs = NULL; 2585 2586 /* benign subspace trick */ 2587 pcbddc->benign_change = 0; 2588 pcbddc->benign_compute_correction = PETSC_TRUE; 2589 pcbddc->benign_vec = 0; 2590 pcbddc->benign_original_mat = 0; 2591 pcbddc->benign_sf = 0; 2592 pcbddc->benign_B0 = 0; 2593 pcbddc->benign_n = 0; 2594 pcbddc->benign_p0 = NULL; 2595 pcbddc->benign_p0_lidx = NULL; 2596 pcbddc->benign_p0_gidx = NULL; 2597 pcbddc->benign_null = PETSC_FALSE; 2598 2599 /* create local graph structure */ 2600 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 2601 pcbddc->graphmaxcount = PETSC_MAX_INT; 2602 2603 /* scaling */ 2604 pcbddc->work_scaling = 0; 2605 pcbddc->use_deluxe_scaling = PETSC_FALSE; 2606 2607 /* sub schurs options */ 2608 pcbddc->sub_schurs_rebuild = PETSC_FALSE; 2609 pcbddc->sub_schurs_layers = -1; 2610 pcbddc->sub_schurs_use_useradj = PETSC_FALSE; 2611 2612 pcbddc->computed_rowadj = PETSC_FALSE; 2613 2614 /* adaptivity */ 2615 pcbddc->adaptive_threshold = 0.0; 2616 pcbddc->adaptive_nmax = 0; 2617 pcbddc->adaptive_nmin = 0; 2618 2619 /* function pointers */ 2620 pc->ops->apply = PCApply_BDDC; 2621 pc->ops->applytranspose = PCApplyTranspose_BDDC; 2622 pc->ops->setup = PCSetUp_BDDC; 2623 pc->ops->destroy = PCDestroy_BDDC; 2624 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 2625 pc->ops->view = PCView_BDDC; 2626 pc->ops->applyrichardson = 0; 2627 pc->ops->applysymmetricleft = 0; 2628 pc->ops->applysymmetricright = 0; 2629 pc->ops->presolve = PCPreSolve_BDDC; 2630 pc->ops->postsolve = PCPostSolve_BDDC; 2631 2632 /* composing function */ 2633 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDiscreteGradient_C",PCBDDCSetDiscreteGradient_BDDC);CHKERRQ(ierr); 2634 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",PCBDDCSetDivergenceMat_BDDC);CHKERRQ(ierr); 2635 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);CHKERRQ(ierr); 2636 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 2637 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",PCBDDCSetPrimalVerticesIS_BDDC);CHKERRQ(ierr); 2638 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 2639 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr); 2640 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr); 2641 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr); 2642 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2643 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2644 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2645 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2646 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2647 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2648 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2649 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2650 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 2651 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr); 2652 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 2653 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 2654 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 2655 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 2656 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",PCPreSolveChangeRHS_BDDC);CHKERRQ(ierr); 2657 PetscFunctionReturn(0); 2658 } 2659 2660