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 if (copymode == PETSC_OWN_POINTER) { 991 ierr = PetscFree(xadj);CHKERRQ(ierr); 992 ierr = PetscFree(adjncy);CHKERRQ(ierr); 993 } 994 ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr); 995 PetscFunctionReturn(0); 996 } 997 if (mat_graph->nvtxs_csr == nvtxs && mat_graph->freecsr) { /* we own the data */ 998 if (mat_graph->xadj == xadj && mat_graph->adjncy == adjncy) same_data = PETSC_TRUE; 999 if (!same_data && mat_graph->xadj[nvtxs] == xadj[nvtxs]) { 1000 ierr = PetscMemcmp(xadj,mat_graph->xadj,(nvtxs+1)*sizeof(PetscInt),&same_data);CHKERRQ(ierr); 1001 if (same_data) { 1002 ierr = PetscMemcmp(adjncy,mat_graph->adjncy,xadj[nvtxs]*sizeof(PetscInt),&same_data);CHKERRQ(ierr); 1003 } 1004 } 1005 } 1006 if (!same_data) { 1007 /* free old CSR */ 1008 ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr); 1009 /* get CSR into graph structure */ 1010 if (copymode == PETSC_COPY_VALUES) { 1011 ierr = PetscMalloc1(nvtxs+1,&mat_graph->xadj);CHKERRQ(ierr); 1012 ierr = PetscMalloc1(xadj[nvtxs],&mat_graph->adjncy);CHKERRQ(ierr); 1013 ierr = PetscMemcpy(mat_graph->xadj,xadj,(nvtxs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1014 ierr = PetscMemcpy(mat_graph->adjncy,adjncy,xadj[nvtxs]*sizeof(PetscInt));CHKERRQ(ierr); 1015 mat_graph->freecsr = PETSC_TRUE; 1016 } else if (copymode == PETSC_OWN_POINTER) { 1017 mat_graph->xadj = (PetscInt*)xadj; 1018 mat_graph->adjncy = (PetscInt*)adjncy; 1019 mat_graph->freecsr = PETSC_TRUE; 1020 } else if (copymode == PETSC_USE_POINTER) { 1021 mat_graph->xadj = (PetscInt*)xadj; 1022 mat_graph->adjncy = (PetscInt*)adjncy; 1023 mat_graph->freecsr = PETSC_FALSE; 1024 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported copy mode %d",copymode); 1025 mat_graph->nvtxs_csr = nvtxs; 1026 pcbddc->recompute_topography = PETSC_TRUE; 1027 } 1028 PetscFunctionReturn(0); 1029 } 1030 1031 #undef __FUNCT__ 1032 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph" 1033 /*@ 1034 PCBDDCSetLocalAdjacencyGraph - Set adjacency structure (CSR graph) of the local degrees of freedom. 1035 1036 Not collective 1037 1038 Input Parameters: 1039 + pc - the preconditioning context. 1040 . nvtxs - number of local vertices of the graph (i.e., the number of local dofs). 1041 . xadj, adjncy - the connectivity of the dofs in CSR format. 1042 - copymode - supported modes are PETSC_COPY_VALUES, PETSC_USE_POINTER or PETSC_OWN_POINTER. 1043 1044 Level: intermediate 1045 1046 Notes: A dof is considered connected with all local dofs if xadj[dof+1]-xadj[dof] == 1 and adjncy[xadj[dof]] is negative. 1047 1048 .seealso: PCBDDC,PetscCopyMode 1049 @*/ 1050 PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc,PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode) 1051 { 1052 void (*f)(void) = 0; 1053 PetscErrorCode ierr; 1054 1055 PetscFunctionBegin; 1056 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1057 if (nvtxs) { 1058 PetscValidIntPointer(xadj,3); 1059 PetscValidIntPointer(adjncy,4); 1060 } 1061 ierr = PetscTryMethod(pc,"PCBDDCSetLocalAdjacencyGraph_C",(PC,PetscInt,const PetscInt[],const PetscInt[],PetscCopyMode),(pc,nvtxs,xadj,adjncy,copymode));CHKERRQ(ierr); 1062 /* free arrays if PCBDDC is not the PC type */ 1063 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",&f);CHKERRQ(ierr); 1064 if (!f && copymode == PETSC_OWN_POINTER) { 1065 ierr = PetscFree(xadj);CHKERRQ(ierr); 1066 ierr = PetscFree(adjncy);CHKERRQ(ierr); 1067 } 1068 PetscFunctionReturn(0); 1069 } 1070 1071 #undef __FUNCT__ 1072 #define __FUNCT__ "PCBDDCSetDofsSplittingLocal_BDDC" 1073 static PetscErrorCode PCBDDCSetDofsSplittingLocal_BDDC(PC pc,PetscInt n_is, IS ISForDofs[]) 1074 { 1075 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1076 PetscInt i; 1077 PetscBool isequal = PETSC_FALSE; 1078 PetscErrorCode ierr; 1079 1080 PetscFunctionBegin; 1081 if (pcbddc->n_ISForDofsLocal == n_is) { 1082 for (i=0;i<n_is;i++) { 1083 PetscBool isequalt; 1084 ierr = ISEqual(ISForDofs[i],pcbddc->ISForDofsLocal[i],&isequalt);CHKERRQ(ierr); 1085 if (!isequalt) break; 1086 } 1087 if (i == n_is) isequal = PETSC_TRUE; 1088 } 1089 for (i=0;i<n_is;i++) { 1090 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 1091 } 1092 /* Destroy ISes if they were already set */ 1093 for (i=0;i<pcbddc->n_ISForDofsLocal;i++) { 1094 ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr); 1095 } 1096 ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr); 1097 /* last user setting takes precendence -> destroy any other customization */ 1098 for (i=0;i<pcbddc->n_ISForDofs;i++) { 1099 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 1100 } 1101 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 1102 pcbddc->n_ISForDofs = 0; 1103 /* allocate space then set */ 1104 if (n_is) { 1105 ierr = PetscMalloc1(n_is,&pcbddc->ISForDofsLocal);CHKERRQ(ierr); 1106 } 1107 for (i=0;i<n_is;i++) { 1108 pcbddc->ISForDofsLocal[i] = ISForDofs[i]; 1109 } 1110 pcbddc->n_ISForDofsLocal = n_is; 1111 if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE; 1112 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 1113 PetscFunctionReturn(0); 1114 } 1115 1116 #undef __FUNCT__ 1117 #define __FUNCT__ "PCBDDCSetDofsSplittingLocal" 1118 /*@ 1119 PCBDDCSetDofsSplittingLocal - Set index sets defining fields of the local subdomain matrix 1120 1121 Collective 1122 1123 Input Parameters: 1124 + pc - the preconditioning context 1125 . n_is - number of index sets defining the fields 1126 - ISForDofs - array of IS describing the fields in local ordering 1127 1128 Level: intermediate 1129 1130 Notes: 1131 n_is should be the same among processes. Not all nodes need to be listed: unlisted nodes will belong to the complement field. 1132 1133 .seealso: PCBDDC 1134 @*/ 1135 PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc,PetscInt n_is, IS ISForDofs[]) 1136 { 1137 PetscInt i; 1138 PetscErrorCode ierr; 1139 1140 PetscFunctionBegin; 1141 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1142 PetscValidLogicalCollectiveInt(pc,n_is,2); 1143 for (i=0;i<n_is;i++) { 1144 PetscCheckSameComm(pc,1,ISForDofs[i],3); 1145 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 1146 } 1147 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplittingLocal_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 1148 PetscFunctionReturn(0); 1149 } 1150 1151 #undef __FUNCT__ 1152 #define __FUNCT__ "PCBDDCSetDofsSplitting_BDDC" 1153 static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc,PetscInt n_is, IS ISForDofs[]) 1154 { 1155 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1156 PetscInt i; 1157 PetscBool isequal = PETSC_FALSE; 1158 PetscErrorCode ierr; 1159 1160 PetscFunctionBegin; 1161 if (pcbddc->n_ISForDofs == n_is) { 1162 for (i=0;i<n_is;i++) { 1163 PetscBool isequalt; 1164 ierr = ISEqual(ISForDofs[i],pcbddc->ISForDofs[i],&isequalt);CHKERRQ(ierr); 1165 if (!isequalt) break; 1166 } 1167 if (i == n_is) isequal = PETSC_TRUE; 1168 } 1169 for (i=0;i<n_is;i++) { 1170 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 1171 } 1172 /* Destroy ISes if they were already set */ 1173 for (i=0;i<pcbddc->n_ISForDofs;i++) { 1174 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 1175 } 1176 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 1177 /* last user setting takes precendence -> destroy any other customization */ 1178 for (i=0;i<pcbddc->n_ISForDofsLocal;i++) { 1179 ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr); 1180 } 1181 ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr); 1182 pcbddc->n_ISForDofsLocal = 0; 1183 /* allocate space then set */ 1184 if (n_is) { 1185 ierr = PetscMalloc1(n_is,&pcbddc->ISForDofs);CHKERRQ(ierr); 1186 } 1187 for (i=0;i<n_is;i++) { 1188 pcbddc->ISForDofs[i] = ISForDofs[i]; 1189 } 1190 pcbddc->n_ISForDofs = n_is; 1191 if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE; 1192 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 1193 PetscFunctionReturn(0); 1194 } 1195 1196 #undef __FUNCT__ 1197 #define __FUNCT__ "PCBDDCSetDofsSplitting" 1198 /*@ 1199 PCBDDCSetDofsSplitting - Set index sets defining fields of the global matrix 1200 1201 Collective 1202 1203 Input Parameters: 1204 + pc - the preconditioning context 1205 . n_is - number of index sets defining the fields 1206 - ISForDofs - array of IS describing the fields in global ordering 1207 1208 Level: intermediate 1209 1210 Notes: 1211 Any process can list any global node. Not all nodes need to be listed: unlisted nodes will belong to the complement field. 1212 1213 .seealso: PCBDDC 1214 @*/ 1215 PetscErrorCode PCBDDCSetDofsSplitting(PC pc,PetscInt n_is, IS ISForDofs[]) 1216 { 1217 PetscInt i; 1218 PetscErrorCode ierr; 1219 1220 PetscFunctionBegin; 1221 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1222 PetscValidLogicalCollectiveInt(pc,n_is,2); 1223 for (i=0;i<n_is;i++) { 1224 PetscCheckSameComm(pc,1,ISForDofs[i],3); 1225 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 1226 } 1227 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 1228 PetscFunctionReturn(0); 1229 } 1230 1231 #undef __FUNCT__ 1232 #define __FUNCT__ "PCPreSolve_BDDC" 1233 /* 1234 PCPreSolve_BDDC - Changes the right hand side and (if necessary) the initial 1235 guess if a transformation of basis approach has been selected. 1236 1237 Input Parameter: 1238 + pc - the preconditioner contex 1239 1240 Application Interface Routine: PCPreSolve() 1241 1242 Notes: 1243 The interface routine PCPreSolve() is not usually called directly by 1244 the user, but instead is called by KSPSolve(). 1245 */ 1246 static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 1247 { 1248 PetscErrorCode ierr; 1249 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1250 PC_IS *pcis = (PC_IS*)(pc->data); 1251 Vec used_vec; 1252 PetscBool save_rhs = PETSC_TRUE, benign_correction_computed; 1253 1254 PetscFunctionBegin; 1255 /* if we are working with CG, one dirichlet solve can be avoided during Krylov iterations */ 1256 if (ksp) { 1257 PetscBool iscg, isgroppcg, ispipecg, ispipecgrr; 1258 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPCG,&iscg);CHKERRQ(ierr); 1259 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPGROPPCG,&isgroppcg);CHKERRQ(ierr); 1260 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPIPECG,&ispipecg);CHKERRQ(ierr); 1261 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPIPECGRR,&ispipecgrr);CHKERRQ(ierr); 1262 if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || (!iscg && !isgroppcg && !ispipecg && !ispipecgrr)) { 1263 ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr); 1264 } 1265 } 1266 if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static) { 1267 ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr); 1268 } 1269 1270 /* Creates parallel work vectors used in presolve */ 1271 if (!pcbddc->original_rhs) { 1272 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 1273 } 1274 if (!pcbddc->temp_solution) { 1275 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);CHKERRQ(ierr); 1276 } 1277 1278 pcbddc->temp_solution_used = PETSC_FALSE; 1279 if (x) { 1280 ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr); 1281 used_vec = x; 1282 } else { /* it can only happen when calling PCBDDCMatFETIDPGetRHS */ 1283 ierr = PetscObjectReference((PetscObject)pcbddc->temp_solution);CHKERRQ(ierr); 1284 used_vec = pcbddc->temp_solution; 1285 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 1286 pcbddc->temp_solution_used = PETSC_TRUE; 1287 ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1288 save_rhs = PETSC_FALSE; 1289 pcbddc->eliminate_dirdofs = PETSC_TRUE; 1290 } 1291 1292 /* hack into ksp data structure since PCPreSolve comes earlier than setting to zero the guess in src/ksp/ksp/interface/itfunc.c */ 1293 if (ksp) { 1294 /* store the flag for the initial guess since it will be restored back during PCPostSolve_BDDC */ 1295 ierr = KSPGetInitialGuessNonzero(ksp,&pcbddc->ksp_guess_nonzero);CHKERRQ(ierr); 1296 if (!pcbddc->ksp_guess_nonzero) { 1297 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 1298 } 1299 } 1300 1301 pcbddc->rhs_change = PETSC_FALSE; 1302 /* Take into account zeroed rows -> change rhs and store solution removed */ 1303 if (rhs && pcbddc->eliminate_dirdofs) { 1304 IS dirIS = NULL; 1305 1306 /* DirichletBoundariesLocal may not be consistent among neighbours; gets a dirichlet dofs IS from graph (may be cached) */ 1307 ierr = PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph,&dirIS);CHKERRQ(ierr); 1308 if (dirIS) { 1309 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 1310 PetscInt dirsize,i,*is_indices; 1311 PetscScalar *array_x; 1312 const PetscScalar *array_diagonal; 1313 1314 ierr = MatGetDiagonal(pc->pmat,pcis->vec1_global);CHKERRQ(ierr); 1315 ierr = VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);CHKERRQ(ierr); 1316 ierr = VecScatterBegin(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1317 ierr = VecScatterEnd(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1318 ierr = VecScatterBegin(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1319 ierr = VecScatterEnd(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1320 ierr = ISGetLocalSize(dirIS,&dirsize);CHKERRQ(ierr); 1321 ierr = VecGetArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 1322 ierr = VecGetArrayRead(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 1323 ierr = ISGetIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1324 for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]]; 1325 ierr = ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1326 ierr = VecRestoreArrayRead(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 1327 ierr = VecRestoreArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 1328 ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1329 ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1330 pcbddc->rhs_change = PETSC_TRUE; 1331 ierr = ISDestroy(&dirIS);CHKERRQ(ierr); 1332 } 1333 } 1334 1335 /* remove the computed solution or the initial guess from the rhs */ 1336 if (pcbddc->rhs_change || (ksp && pcbddc->ksp_guess_nonzero) ) { 1337 /* save the original rhs */ 1338 if (save_rhs) { 1339 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1340 save_rhs = PETSC_FALSE; 1341 } 1342 pcbddc->rhs_change = PETSC_TRUE; 1343 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 1344 ierr = MatMultAdd(pc->mat,used_vec,pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1345 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 1346 ierr = VecCopy(used_vec,pcbddc->temp_solution);CHKERRQ(ierr); 1347 pcbddc->temp_solution_used = PETSC_TRUE; 1348 if (ksp) { 1349 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_FALSE);CHKERRQ(ierr); 1350 } 1351 } 1352 ierr = VecDestroy(&used_vec);CHKERRQ(ierr); 1353 1354 /* compute initial vector in benign space if needed 1355 and remove non-benign solution from the rhs */ 1356 benign_correction_computed = PETSC_FALSE; 1357 if (rhs && pcbddc->benign_compute_correction && pcbddc->benign_have_null) { 1358 /* compute u^*_h using ideas similar to those in Xuemin Tu's PhD thesis (see Section 4.8.1) 1359 Recursively apply BDDC in the multilevel case */ 1360 if (!pcbddc->benign_vec) { 1361 ierr = VecDuplicate(rhs,&pcbddc->benign_vec);CHKERRQ(ierr); 1362 } 1363 pcbddc->benign_apply_coarse_only = PETSC_TRUE; 1364 if (!pcbddc->benign_skip_correction) { 1365 ierr = PCApply_BDDC(pc,rhs,pcbddc->benign_vec);CHKERRQ(ierr); 1366 benign_correction_computed = PETSC_TRUE; 1367 if (pcbddc->temp_solution_used) { 1368 ierr = VecAXPY(pcbddc->temp_solution,1.0,pcbddc->benign_vec);CHKERRQ(ierr); 1369 } 1370 ierr = VecScale(pcbddc->benign_vec,-1.0);CHKERRQ(ierr); 1371 /* store the original rhs if not done earlier */ 1372 if (save_rhs) { 1373 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1374 } 1375 if (pcbddc->rhs_change) { 1376 ierr = MatMultAdd(pc->mat,pcbddc->benign_vec,rhs,rhs);CHKERRQ(ierr); 1377 } else { 1378 ierr = MatMultAdd(pc->mat,pcbddc->benign_vec,pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1379 } 1380 pcbddc->rhs_change = PETSC_TRUE; 1381 } 1382 pcbddc->benign_apply_coarse_only = PETSC_FALSE; 1383 } 1384 1385 /* dbg output */ 1386 if (pcbddc->dbg_flag && benign_correction_computed) { 1387 Vec v; 1388 ierr = VecDuplicate(pcis->vec1_global,&v);CHKERRQ(ierr); 1389 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,rhs,v);CHKERRQ(ierr); 1390 ierr = PCBDDCBenignGetOrSetP0(pc,v,PETSC_TRUE);CHKERRQ(ierr); 1391 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"LEVEL %d: is the correction benign?\n",pcbddc->current_level);CHKERRQ(ierr); 1392 ierr = PetscScalarView(pcbddc->benign_n,pcbddc->benign_p0,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)));CHKERRQ(ierr); 1393 ierr = VecDestroy(&v);CHKERRQ(ierr); 1394 } 1395 1396 /* set initial guess if using PCG */ 1397 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 1398 if (x && pcbddc->use_exact_dirichlet_trick) { 1399 ierr = VecSet(x,0.0);CHKERRQ(ierr); 1400 if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) { 1401 if (benign_correction_computed) { /* we have already saved the changed rhs */ 1402 ierr = VecLockPop(pcis->vec1_global);CHKERRQ(ierr); 1403 } else { 1404 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,rhs,pcis->vec1_global);CHKERRQ(ierr); 1405 } 1406 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_global,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1407 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec1_global,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1408 } else { 1409 ierr = VecScatterBegin(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1410 ierr = VecScatterEnd(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1411 } 1412 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1413 if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) { 1414 ierr = VecSet(pcis->vec1_global,0.);CHKERRQ(ierr); 1415 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1416 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1417 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcis->vec1_global,x);CHKERRQ(ierr); 1418 } else { 1419 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1420 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1421 } 1422 if (ksp) { 1423 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); 1424 } 1425 pcbddc->exact_dirichlet_trick_app = PETSC_TRUE; 1426 } else if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior && benign_correction_computed && pcbddc->use_exact_dirichlet_trick) { 1427 ierr = VecLockPop(pcis->vec1_global);CHKERRQ(ierr); 1428 } 1429 PetscFunctionReturn(0); 1430 } 1431 1432 #undef __FUNCT__ 1433 #define __FUNCT__ "PCPostSolve_BDDC" 1434 /* 1435 PCPostSolve_BDDC - Changes the computed solution if a transformation of basis 1436 approach has been selected. Also, restores rhs to its original state. 1437 1438 Input Parameter: 1439 + pc - the preconditioner contex 1440 1441 Application Interface Routine: PCPostSolve() 1442 1443 Notes: 1444 The interface routine PCPostSolve() is not usually called directly by 1445 the user, but instead is called by KSPSolve(). 1446 */ 1447 static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 1448 { 1449 PetscErrorCode ierr; 1450 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1451 1452 PetscFunctionBegin; 1453 /* add solution removed in presolve */ 1454 if (x && pcbddc->rhs_change) { 1455 if (pcbddc->temp_solution_used) { 1456 ierr = VecAXPY(x,1.0,pcbddc->temp_solution);CHKERRQ(ierr); 1457 } else if (pcbddc->benign_compute_correction) { 1458 ierr = VecAXPY(x,-1.0,pcbddc->benign_vec);CHKERRQ(ierr); 1459 } 1460 } 1461 pcbddc->temp_solution_used = PETSC_FALSE; 1462 1463 /* restore rhs to its original state (not needed for FETI-DP) */ 1464 if (rhs && pcbddc->rhs_change) { 1465 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1466 } 1467 pcbddc->rhs_change = PETSC_FALSE; 1468 /* restore ksp guess state */ 1469 if (ksp) { 1470 ierr = KSPSetInitialGuessNonzero(ksp,pcbddc->ksp_guess_nonzero);CHKERRQ(ierr); 1471 } 1472 /* reset flag for exact dirichlet trick */ 1473 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 1474 PetscFunctionReturn(0); 1475 } 1476 #undef __FUNCT__ 1477 #define __FUNCT__ "PCSetUp_BDDC" 1478 /* 1479 PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner 1480 by setting data structures and options. 1481 1482 Input Parameter: 1483 + pc - the preconditioner context 1484 1485 Application Interface Routine: PCSetUp() 1486 1487 Notes: 1488 The interface routine PCSetUp() is not usually called directly by 1489 the user, but instead is called by PCApply() if necessary. 1490 */ 1491 PetscErrorCode PCSetUp_BDDC(PC pc) 1492 { 1493 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 1494 PCBDDCSubSchurs sub_schurs; 1495 Mat_IS* matis; 1496 MatNullSpace nearnullspace; 1497 IS zerodiag = NULL; 1498 PetscInt nrows,ncols; 1499 PetscBool computesubschurs; 1500 PetscBool computeconstraintsmatrix; 1501 PetscBool new_nearnullspace_provided,ismatis; 1502 PetscErrorCode ierr; 1503 1504 PetscFunctionBegin; 1505 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&ismatis);CHKERRQ(ierr); 1506 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PCBDDC preconditioner requires matrix of type MATIS"); 1507 ierr = MatGetSize(pc->pmat,&nrows,&ncols);CHKERRQ(ierr); 1508 if (nrows != ncols) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCBDDC preconditioner requires a square preconditioning matrix"); 1509 matis = (Mat_IS*)pc->pmat->data; 1510 /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */ 1511 /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup 1512 Also, BDDC builds its own KSP for the Dirichlet problem */ 1513 if (pc->setupcalled && pc->flag != SAME_NONZERO_PATTERN) pcbddc->recompute_topography = PETSC_TRUE; 1514 if (pcbddc->recompute_topography) pcbddc->graphanalyzed = PETSC_FALSE; 1515 computeconstraintsmatrix = PETSC_FALSE; 1516 1517 /* check parameters' compatibility */ 1518 if (!pcbddc->use_deluxe_scaling) pcbddc->deluxe_zerorows = PETSC_FALSE; 1519 pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_threshold > 0.0); 1520 pcbddc->adaptive_userdefined = (PetscBool)(pcbddc->adaptive_selection && pcbddc->adaptive_userdefined); 1521 if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE; 1522 1523 computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling); 1524 if (pcbddc->switch_static) { 1525 PetscBool ismatis; 1526 ierr = PetscObjectTypeCompare((PetscObject)pc->mat,MATIS,&ismatis);CHKERRQ(ierr); 1527 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"When the static switch is one, the iteration matrix should be of type MATIS"); 1528 } 1529 1530 /* activate all connected components if the netflux has been requested */ 1531 if (pcbddc->compute_nonetflux) { 1532 pcbddc->use_vertices = PETSC_TRUE; 1533 pcbddc->use_edges = PETSC_TRUE; 1534 pcbddc->use_faces = PETSC_TRUE; 1535 } 1536 1537 /* Get stdout for dbg */ 1538 if (pcbddc->dbg_flag) { 1539 if (!pcbddc->dbg_viewer) { 1540 pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)); 1541 ierr = PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer);CHKERRQ(ierr); 1542 } 1543 ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1544 } 1545 1546 /* process topology information */ 1547 if (pcbddc->recompute_topography) { 1548 ierr = PCBDDCComputeLocalTopologyInfo(pc);CHKERRQ(ierr); 1549 if (pcbddc->discretegradient) { 1550 ierr = PCBDDCNedelecSupport(pc);CHKERRQ(ierr); 1551 } 1552 } 1553 1554 /* change basis if requested by the user */ 1555 if (pcbddc->user_ChangeOfBasisMatrix) { 1556 /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */ 1557 pcbddc->use_change_of_basis = PETSC_FALSE; 1558 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 1559 } else { 1560 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1561 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1562 pcbddc->local_mat = matis->A; 1563 } 1564 1565 /* detect local disconnected subdomains if requested */ 1566 if (pcbddc->detect_disconnected && pcbddc->recompute_topography) { 1567 PetscInt i; 1568 for (i=0;i<pcbddc->n_local_subs;i++) { 1569 ierr = ISDestroy(&pcbddc->local_subs[i]);CHKERRQ(ierr); 1570 } 1571 ierr = PetscFree(pcbddc->local_subs);CHKERRQ(ierr); 1572 ierr = MatDetectDisconnectedComponents(pcbddc->local_mat,PETSC_FALSE,&pcbddc->n_local_subs,&pcbddc->local_subs);CHKERRQ(ierr); 1573 } 1574 1575 /* 1576 Compute change of basis on local pressures (aka zerodiag dofs) with the benign trick 1577 This should come earlier then PCISSetUp for extracting the correct subdomain matrices 1578 */ 1579 ierr = PCBDDCBenignShellMat(pc,PETSC_TRUE);CHKERRQ(ierr); 1580 if (pcbddc->benign_saddle_point) { 1581 PC_IS* pcis = (PC_IS*)pc->data; 1582 1583 if (pcbddc->user_ChangeOfBasisMatrix || pcbddc->use_change_of_basis || !computesubschurs) pcbddc->benign_change_explicit = PETSC_TRUE; 1584 /* detect local saddle point and change the basis in pcbddc->local_mat (TODO: reuse case) */ 1585 ierr = PCBDDCBenignDetectSaddlePoint(pc,&zerodiag);CHKERRQ(ierr); 1586 /* pop B0 mat from local mat */ 1587 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 1588 /* give pcis a hint to not reuse submatrices during PCISCreate */ 1589 if (pc->flag == SAME_NONZERO_PATTERN && pcis->reusesubmatrices == PETSC_TRUE) { 1590 if (pcbddc->benign_n && (pcbddc->benign_change_explicit || pcbddc->dbg_flag)) { 1591 pcis->reusesubmatrices = PETSC_FALSE; 1592 } else { 1593 pcis->reusesubmatrices = PETSC_TRUE; 1594 } 1595 } else { 1596 pcis->reusesubmatrices = PETSC_FALSE; 1597 } 1598 } 1599 1600 /* propagate relevant information -> TODO remove*/ 1601 #if !defined(PETSC_USE_COMPLEX) /* workaround for reals */ 1602 if (matis->A->symmetric_set) { 1603 ierr = MatSetOption(pcbddc->local_mat,MAT_HERMITIAN,matis->A->symmetric);CHKERRQ(ierr); 1604 } 1605 #endif 1606 if (matis->A->symmetric_set) { 1607 ierr = MatSetOption(pcbddc->local_mat,MAT_SYMMETRIC,matis->A->symmetric);CHKERRQ(ierr); 1608 } 1609 if (matis->A->spd_set) { 1610 ierr = MatSetOption(pcbddc->local_mat,MAT_SPD,matis->A->spd);CHKERRQ(ierr); 1611 } 1612 1613 /* Set up all the "iterative substructuring" common block without computing solvers */ 1614 { 1615 Mat temp_mat; 1616 1617 temp_mat = matis->A; 1618 matis->A = pcbddc->local_mat; 1619 ierr = PCISSetUp(pc,PETSC_FALSE);CHKERRQ(ierr); 1620 pcbddc->local_mat = matis->A; 1621 matis->A = temp_mat; 1622 } 1623 1624 /* Analyze interface */ 1625 if (!pcbddc->graphanalyzed) { 1626 ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr); 1627 computeconstraintsmatrix = PETSC_TRUE; 1628 if (pcbddc->adaptive_selection && !pcbddc->use_deluxe_scaling && !pcbddc->mat_graph->twodim) { 1629 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot compute the adaptive primal space for a problem with 3D edges without deluxe scaling"); 1630 } 1631 if (pcbddc->compute_nonetflux) { 1632 MatNullSpace nnfnnsp; 1633 1634 ierr = PCBDDCComputeNoNetFlux(pc->pmat,pcbddc->divudotp,pcbddc->divudotp_trans,pcbddc->divudotp_vl2l,pcbddc->mat_graph,&nnfnnsp);CHKERRQ(ierr); 1635 /* TODO what if a nearnullspace is already attached? */ 1636 ierr = MatSetNearNullSpace(pc->pmat,nnfnnsp);CHKERRQ(ierr); 1637 ierr = MatNullSpaceDestroy(&nnfnnsp);CHKERRQ(ierr); 1638 } 1639 } 1640 1641 /* check existence of a divergence free extension, i.e. 1642 b(v_I,p_0) = 0 for all v_I (raise error if not). 1643 Also, check that PCBDDCBenignGetOrSetP0 works */ 1644 #if defined(PETSC_USE_DEBUG) 1645 if (pcbddc->benign_saddle_point) { 1646 ierr = PCBDDCBenignCheck(pc,zerodiag);CHKERRQ(ierr); 1647 } 1648 #endif 1649 ierr = ISDestroy(&zerodiag);CHKERRQ(ierr); 1650 1651 /* Setup local dirichlet solver ksp_D and sub_schurs solvers */ 1652 if (computesubschurs && pcbddc->recompute_topography) { 1653 ierr = PCBDDCInitSubSchurs(pc);CHKERRQ(ierr); 1654 } 1655 /* SetUp Scaling operator (scaling matrices could be needed in SubSchursSetUp)*/ 1656 if (!pcbddc->use_deluxe_scaling) { 1657 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1658 } 1659 1660 /* finish setup solvers and do adaptive selection of constraints */ 1661 sub_schurs = pcbddc->sub_schurs; 1662 if (sub_schurs && sub_schurs->schur_explicit) { 1663 if (computesubschurs) { 1664 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1665 } 1666 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1667 } else { 1668 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1669 if (computesubschurs) { 1670 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1671 } 1672 } 1673 if (pcbddc->adaptive_selection) { 1674 ierr = PCBDDCAdaptiveSelection(pc);CHKERRQ(ierr); 1675 computeconstraintsmatrix = PETSC_TRUE; 1676 } 1677 1678 /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */ 1679 new_nearnullspace_provided = PETSC_FALSE; 1680 ierr = MatGetNearNullSpace(pc->pmat,&nearnullspace);CHKERRQ(ierr); 1681 if (pcbddc->onearnullspace) { /* already used nearnullspace */ 1682 if (!nearnullspace) { /* near null space attached to mat has been destroyed */ 1683 new_nearnullspace_provided = PETSC_TRUE; 1684 } else { 1685 /* determine if the two nullspaces are different (should be lightweight) */ 1686 if (nearnullspace != pcbddc->onearnullspace) { 1687 new_nearnullspace_provided = PETSC_TRUE; 1688 } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */ 1689 PetscInt i; 1690 const Vec *nearnullvecs; 1691 PetscObjectState state; 1692 PetscInt nnsp_size; 1693 ierr = MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 1694 for (i=0;i<nnsp_size;i++) { 1695 ierr = PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);CHKERRQ(ierr); 1696 if (pcbddc->onearnullvecs_state[i] != state) { 1697 new_nearnullspace_provided = PETSC_TRUE; 1698 break; 1699 } 1700 } 1701 } 1702 } 1703 } else { 1704 if (!nearnullspace) { /* both nearnullspaces are null */ 1705 new_nearnullspace_provided = PETSC_FALSE; 1706 } else { /* nearnullspace attached later */ 1707 new_nearnullspace_provided = PETSC_TRUE; 1708 } 1709 } 1710 1711 /* Setup constraints and related work vectors */ 1712 /* reset primal space flags */ 1713 pcbddc->new_primal_space = PETSC_FALSE; 1714 pcbddc->new_primal_space_local = PETSC_FALSE; 1715 if (computeconstraintsmatrix || new_nearnullspace_provided) { 1716 /* It also sets the primal space flags */ 1717 ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr); 1718 /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */ 1719 ierr = PCBDDCSetUpLocalWorkVectors(pc);CHKERRQ(ierr); 1720 } 1721 1722 if (pcbddc->use_change_of_basis) { 1723 PC_IS *pcis = (PC_IS*)(pc->data); 1724 1725 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 1726 if (pcbddc->benign_change) { 1727 ierr = MatDestroy(&pcbddc->benign_B0);CHKERRQ(ierr); 1728 /* pop B0 from pcbddc->local_mat */ 1729 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 1730 } 1731 /* get submatrices */ 1732 ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr); 1733 ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr); 1734 ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr); 1735 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);CHKERRQ(ierr); 1736 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);CHKERRQ(ierr); 1737 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);CHKERRQ(ierr); 1738 /* set flag in pcis to not reuse submatrices during PCISCreate */ 1739 pcis->reusesubmatrices = PETSC_FALSE; 1740 } else if (!pcbddc->user_ChangeOfBasisMatrix && !pcbddc->benign_change) { 1741 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1742 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1743 pcbddc->local_mat = matis->A; 1744 } 1745 /* SetUp coarse and local Neumann solvers */ 1746 ierr = PCBDDCSetUpSolvers(pc);CHKERRQ(ierr); 1747 /* SetUp Scaling operator */ 1748 if (pcbddc->use_deluxe_scaling) { 1749 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1750 } 1751 1752 /* mark topography as done */ 1753 pcbddc->recompute_topography = PETSC_FALSE; 1754 1755 /* wrap pcis->A_IB and pcis->A_BI if we did not change explicitly the variables on the pressures */ 1756 ierr = PCBDDCBenignShellMat(pc,PETSC_FALSE);CHKERRQ(ierr); 1757 1758 if (pcbddc->dbg_flag) { 1759 ierr = PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1760 } 1761 PetscFunctionReturn(0); 1762 } 1763 1764 /* 1765 PCApply_BDDC - Applies the BDDC operator to a vector. 1766 1767 Input Parameters: 1768 + pc - the preconditioner context 1769 - r - input vector (global) 1770 1771 Output Parameter: 1772 . z - output vector (global) 1773 1774 Application Interface Routine: PCApply() 1775 */ 1776 #undef __FUNCT__ 1777 #define __FUNCT__ "PCApply_BDDC" 1778 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z) 1779 { 1780 PC_IS *pcis = (PC_IS*)(pc->data); 1781 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1782 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1783 PetscErrorCode ierr; 1784 const PetscScalar one = 1.0; 1785 const PetscScalar m_one = -1.0; 1786 const PetscScalar zero = 0.0; 1787 1788 /* This code is similar to that provided in nn.c for PCNN 1789 NN interface preconditioner changed to BDDC 1790 Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */ 1791 1792 PetscFunctionBegin; 1793 if (pcbddc->ChangeOfBasisMatrix) { 1794 Vec swap; 1795 1796 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);CHKERRQ(ierr); 1797 swap = pcbddc->work_change; 1798 pcbddc->work_change = r; 1799 r = swap; 1800 /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */ 1801 if (pcbddc->benign_apply_coarse_only && pcbddc->exact_dirichlet_trick_app && pcbddc->change_interior) { 1802 ierr = VecCopy(r,pcis->vec1_global);CHKERRQ(ierr); 1803 ierr = VecLockPush(pcis->vec1_global);CHKERRQ(ierr); 1804 } 1805 } 1806 if (pcbddc->benign_have_null) { /* get p0 from r */ 1807 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 1808 } 1809 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 1810 ierr = VecCopy(r,z);CHKERRQ(ierr); 1811 /* First Dirichlet solve */ 1812 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1813 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1814 /* 1815 Assembling right hand side for BDDC operator 1816 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1817 - pcis->vec1_B the interface part of the global vector z 1818 */ 1819 if (n_D) { 1820 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1821 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1822 if (pcbddc->switch_static) { 1823 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1824 1825 ierr = VecSet(pcis->vec1_N,0.);CHKERRQ(ierr); 1826 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1827 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1828 if (!pcbddc->switch_static_change) { 1829 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1830 } else { 1831 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1832 ierr = MatMult(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1833 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1834 } 1835 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1836 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1837 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1838 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1839 } else { 1840 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1841 } 1842 } else { 1843 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1844 } 1845 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1846 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1847 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1848 } else { 1849 if (!pcbddc->benign_apply_coarse_only) { 1850 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1851 } 1852 } 1853 1854 /* Apply interface preconditioner 1855 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1856 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);CHKERRQ(ierr); 1857 1858 /* Apply transpose of partition of unity operator */ 1859 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1860 1861 /* Second Dirichlet solve and assembling of output */ 1862 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1863 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1864 if (n_B) { 1865 if (pcbddc->switch_static) { 1866 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1867 1868 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1869 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1870 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1871 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1872 if (!pcbddc->switch_static_change) { 1873 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1874 } else { 1875 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1876 ierr = MatMult(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1877 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1878 } 1879 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1880 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1881 } else { 1882 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1883 } 1884 } else if (pcbddc->switch_static) { /* n_B is zero */ 1885 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1886 1887 if (!pcbddc->switch_static_change) { 1888 ierr = MatMult(matis->A,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 1889 } else { 1890 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);CHKERRQ(ierr); 1891 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1892 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);CHKERRQ(ierr); 1893 } 1894 } 1895 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1896 1897 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 1898 if (pcbddc->switch_static) { 1899 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 1900 } else { 1901 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 1902 } 1903 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1904 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1905 } else { 1906 if (pcbddc->switch_static) { 1907 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 1908 } else { 1909 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1910 } 1911 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1912 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1913 } 1914 if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */ 1915 if (pcbddc->benign_apply_coarse_only) { 1916 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 1917 } 1918 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 1919 } 1920 1921 if (pcbddc->ChangeOfBasisMatrix) { 1922 pcbddc->work_change = r; 1923 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 1924 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 1925 } 1926 PetscFunctionReturn(0); 1927 } 1928 1929 /* 1930 PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector. 1931 1932 Input Parameters: 1933 + pc - the preconditioner context 1934 - r - input vector (global) 1935 1936 Output Parameter: 1937 . z - output vector (global) 1938 1939 Application Interface Routine: PCApplyTranspose() 1940 */ 1941 #undef __FUNCT__ 1942 #define __FUNCT__ "PCApplyTranspose_BDDC" 1943 PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z) 1944 { 1945 PC_IS *pcis = (PC_IS*)(pc->data); 1946 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1947 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1948 PetscErrorCode ierr; 1949 const PetscScalar one = 1.0; 1950 const PetscScalar m_one = -1.0; 1951 const PetscScalar zero = 0.0; 1952 1953 PetscFunctionBegin; 1954 if (pcbddc->ChangeOfBasisMatrix) { 1955 Vec swap; 1956 1957 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);CHKERRQ(ierr); 1958 swap = pcbddc->work_change; 1959 pcbddc->work_change = r; 1960 r = swap; 1961 /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */ 1962 if (pcbddc->benign_apply_coarse_only && pcbddc->exact_dirichlet_trick_app && pcbddc->change_interior) { 1963 ierr = VecCopy(r,pcis->vec1_global);CHKERRQ(ierr); 1964 ierr = VecLockPush(pcis->vec1_global);CHKERRQ(ierr); 1965 } 1966 } 1967 if (pcbddc->benign_have_null) { /* get p0 from r */ 1968 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 1969 } 1970 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 1971 ierr = VecCopy(r,z);CHKERRQ(ierr); 1972 /* First Dirichlet solve */ 1973 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1974 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1975 /* 1976 Assembling right hand side for BDDC operator 1977 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1978 - pcis->vec1_B the interface part of the global vector z 1979 */ 1980 if (n_D) { 1981 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1982 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1983 if (pcbddc->switch_static) { 1984 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1985 1986 ierr = VecSet(pcis->vec1_N,0.);CHKERRQ(ierr); 1987 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1988 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1989 if (!pcbddc->switch_static_change) { 1990 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1991 } else { 1992 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1993 ierr = MatMultTranspose(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1994 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1995 } 1996 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1997 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1998 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1999 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2000 } else { 2001 ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 2002 } 2003 } else { 2004 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 2005 } 2006 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2007 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2008 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 2009 } else { 2010 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 2011 } 2012 2013 /* Apply interface preconditioner 2014 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 2015 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr); 2016 2017 /* Apply transpose of partition of unity operator */ 2018 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 2019 2020 /* Second Dirichlet solve and assembling of output */ 2021 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2022 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2023 if (n_B) { 2024 if (pcbddc->switch_static) { 2025 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 2026 2027 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2028 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2029 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2030 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2031 if (!pcbddc->switch_static_change) { 2032 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2033 } else { 2034 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2035 ierr = MatMultTranspose(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 2036 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2037 } 2038 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2039 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2040 } else { 2041 ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 2042 } 2043 } else if (pcbddc->switch_static) { /* n_B is zero */ 2044 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 2045 2046 if (!pcbddc->switch_static_change) { 2047 ierr = MatMultTranspose(matis->A,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 2048 } else { 2049 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);CHKERRQ(ierr); 2050 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2051 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);CHKERRQ(ierr); 2052 } 2053 } 2054 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 2055 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 2056 if (pcbddc->switch_static) { 2057 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 2058 } else { 2059 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 2060 } 2061 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2062 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2063 } else { 2064 if (pcbddc->switch_static) { 2065 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 2066 } else { 2067 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 2068 } 2069 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2070 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2071 } 2072 if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */ 2073 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 2074 } 2075 if (pcbddc->ChangeOfBasisMatrix) { 2076 pcbddc->work_change = r; 2077 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 2078 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 2079 } 2080 PetscFunctionReturn(0); 2081 } 2082 2083 #undef __FUNCT__ 2084 #define __FUNCT__ "PCDestroy_BDDC" 2085 PetscErrorCode PCDestroy_BDDC(PC pc) 2086 { 2087 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 2088 PetscErrorCode ierr; 2089 2090 PetscFunctionBegin; 2091 /* free BDDC custom data */ 2092 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 2093 /* destroy objects related to topography */ 2094 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 2095 /* free allocated graph structure */ 2096 ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr); 2097 /* free allocated sub schurs structure */ 2098 ierr = PetscFree(pcbddc->sub_schurs);CHKERRQ(ierr); 2099 /* destroy objects for scaling operator */ 2100 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 2101 ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr); 2102 /* free solvers stuff */ 2103 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 2104 /* free global vectors needed in presolve */ 2105 ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr); 2106 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 2107 /* free data created by PCIS */ 2108 ierr = PCISDestroy(pc);CHKERRQ(ierr); 2109 /* remove functions */ 2110 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDiscreteGradient_C",NULL);CHKERRQ(ierr); 2111 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",NULL);CHKERRQ(ierr); 2112 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);CHKERRQ(ierr); 2113 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 2114 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",NULL);CHKERRQ(ierr); 2115 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 2116 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr); 2117 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr); 2118 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr); 2119 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2120 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2121 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2122 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2123 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2124 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2125 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2126 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2127 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 2128 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr); 2129 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 2130 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 2131 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 2132 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 2133 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",NULL);CHKERRQ(ierr); 2134 /* Free the private data structure */ 2135 ierr = PetscFree(pc->data);CHKERRQ(ierr); 2136 PetscFunctionReturn(0); 2137 } 2138 2139 #undef __FUNCT__ 2140 #define __FUNCT__ "PCPreSolveChangeRHS_BDDC" 2141 static PetscErrorCode PCPreSolveChangeRHS_BDDC(PC pc, PetscBool* change) 2142 { 2143 PetscFunctionBegin; 2144 *change = PETSC_TRUE; 2145 PetscFunctionReturn(0); 2146 } 2147 2148 #undef __FUNCT__ 2149 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC" 2150 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 2151 { 2152 FETIDPMat_ctx mat_ctx; 2153 Vec work; 2154 PC_IS* pcis; 2155 PC_BDDC* pcbddc; 2156 PetscErrorCode ierr; 2157 2158 PetscFunctionBegin; 2159 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2160 pcis = (PC_IS*)mat_ctx->pc->data; 2161 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2162 2163 /* 2164 change of basis for physical rhs if needed 2165 It also changes the rhs in case of dirichlet boundaries 2166 */ 2167 if (!pcbddc->original_rhs) { 2168 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 2169 } 2170 ierr = VecCopy(standard_rhs,pcbddc->original_rhs);CHKERRQ(ierr); 2171 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,pcbddc->original_rhs,NULL);CHKERRQ(ierr); 2172 if (pcbddc->ChangeOfBasisMatrix) { 2173 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,pcbddc->original_rhs,pcbddc->work_change);CHKERRQ(ierr); 2174 work = pcbddc->work_change; 2175 } else { 2176 work = pcbddc->original_rhs; 2177 } 2178 /* store vectors for computation of fetidp final solution */ 2179 ierr = VecScatterBegin(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2180 ierr = VecScatterEnd(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2181 /* scale rhs since it should be unassembled */ 2182 /* TODO use counter scaling? (also below) */ 2183 ierr = VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2184 ierr = VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2185 /* Apply partition of unity */ 2186 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2187 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2188 if (!pcbddc->switch_static) { 2189 /* compute partially subassembled Schur complement right-hand side */ 2190 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2191 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 2192 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 2193 ierr = VecSet(work,0.0);CHKERRQ(ierr); 2194 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2195 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2196 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2197 ierr = VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2198 ierr = VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2199 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2200 } 2201 /* BDDC rhs */ 2202 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 2203 if (pcbddc->switch_static) { 2204 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2205 } 2206 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2207 /* apply BDDC */ 2208 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2209 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2210 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 2211 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 2212 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 2213 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2214 ierr = VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2215 PetscFunctionReturn(0); 2216 } 2217 2218 #undef __FUNCT__ 2219 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS" 2220 /*@ 2221 PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETI-DP linear system using the physical right-hand side 2222 2223 Collective 2224 2225 Input Parameters: 2226 + fetidp_mat - the FETI-DP matrix object obtained by a call to PCBDDCCreateFETIDPOperators 2227 - standard_rhs - the right-hand side of the original linear system 2228 2229 Output Parameters: 2230 . fetidp_flux_rhs - the right-hand side for the FETI-DP linear system 2231 2232 Level: developer 2233 2234 Notes: 2235 2236 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetSolution 2237 @*/ 2238 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 2239 { 2240 FETIDPMat_ctx mat_ctx; 2241 PetscErrorCode ierr; 2242 2243 PetscFunctionBegin; 2244 PetscValidHeaderSpecific(fetidp_mat,MAT_CLASSID,1); 2245 PetscValidHeaderSpecific(standard_rhs,VEC_CLASSID,2); 2246 PetscValidHeaderSpecific(fetidp_flux_rhs,VEC_CLASSID,3); 2247 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2248 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 2249 PetscFunctionReturn(0); 2250 } 2251 2252 #undef __FUNCT__ 2253 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC" 2254 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2255 { 2256 FETIDPMat_ctx mat_ctx; 2257 PC_IS* pcis; 2258 PC_BDDC* pcbddc; 2259 PetscErrorCode ierr; 2260 2261 PetscFunctionBegin; 2262 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2263 pcis = (PC_IS*)mat_ctx->pc->data; 2264 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2265 2266 /* apply B_delta^T */ 2267 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2268 ierr = VecScatterEnd(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2269 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 2270 /* compute rhs for BDDC application */ 2271 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2272 if (pcbddc->switch_static) { 2273 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2274 } 2275 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2276 /* apply BDDC */ 2277 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2278 /* put values into standard global vector */ 2279 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2280 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2281 if (!pcbddc->switch_static) { 2282 /* compute values into the interior if solved for the partially subassembled Schur complement */ 2283 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 2284 ierr = VecAYPX(pcis->vec1_D,-1.0,mat_ctx->temp_solution_D);CHKERRQ(ierr); 2285 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); 2286 } 2287 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2288 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2289 /* add p0 solution to final solution */ 2290 ierr = PCBDDCBenignGetOrSetP0(mat_ctx->pc,standard_sol,PETSC_FALSE);CHKERRQ(ierr); 2291 if (pcbddc->ChangeOfBasisMatrix) { 2292 Vec v2; 2293 ierr = VecDuplicate(standard_sol,&v2);CHKERRQ(ierr); 2294 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,standard_sol,v2);CHKERRQ(ierr); 2295 ierr = VecCopy(v2,standard_sol);CHKERRQ(ierr); 2296 ierr = VecDestroy(&v2);CHKERRQ(ierr); 2297 } 2298 ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr); 2299 PetscFunctionReturn(0); 2300 } 2301 2302 #undef __FUNCT__ 2303 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution" 2304 /*@ 2305 PCBDDCMatFETIDPGetSolution - Compute the physical solution using the solution of the FETI-DP linear system 2306 2307 Collective 2308 2309 Input Parameters: 2310 + fetidp_mat - the FETI-DP matrix obtained by a call to PCBDDCCreateFETIDPOperators 2311 - fetidp_flux_sol - the solution of the FETI-DP linear system 2312 2313 Output Parameters: 2314 . standard_sol - the solution defined on the physical domain 2315 2316 Level: developer 2317 2318 Notes: 2319 2320 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetRHS 2321 @*/ 2322 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2323 { 2324 FETIDPMat_ctx mat_ctx; 2325 PetscErrorCode ierr; 2326 2327 PetscFunctionBegin; 2328 PetscValidHeaderSpecific(fetidp_mat,MAT_CLASSID,1); 2329 PetscValidHeaderSpecific(fetidp_flux_sol,VEC_CLASSID,2); 2330 PetscValidHeaderSpecific(standard_sol,VEC_CLASSID,3); 2331 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2332 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 2333 PetscFunctionReturn(0); 2334 } 2335 2336 extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec); 2337 extern PetscErrorCode FETIDPMatMultTranspose(Mat,Vec,Vec); 2338 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat); 2339 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec); 2340 extern PetscErrorCode FETIDPPCApplyTranspose(PC,Vec,Vec); 2341 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC); 2342 2343 #undef __FUNCT__ 2344 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC" 2345 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, PetscBool fully_redundant, Mat *fetidp_mat, PC *fetidp_pc) 2346 { 2347 2348 FETIDPMat_ctx fetidpmat_ctx; 2349 Mat newmat; 2350 FETIDPPC_ctx fetidppc_ctx; 2351 PC newpc; 2352 MPI_Comm comm; 2353 PetscErrorCode ierr; 2354 2355 PetscFunctionBegin; 2356 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 2357 /* FETIDP linear matrix */ 2358 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 2359 fetidpmat_ctx->fully_redundant = fully_redundant; 2360 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 2361 ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 2362 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 2363 ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr); 2364 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 2365 ierr = MatSetUp(newmat);CHKERRQ(ierr); 2366 /* FETIDP preconditioner */ 2367 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 2368 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 2369 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 2370 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 2371 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 2372 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 2373 ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 2374 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 2375 ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr); 2376 ierr = PCSetUp(newpc);CHKERRQ(ierr); 2377 /* return pointers for objects created */ 2378 *fetidp_mat=newmat; 2379 *fetidp_pc=newpc; 2380 PetscFunctionReturn(0); 2381 } 2382 2383 #undef __FUNCT__ 2384 #define __FUNCT__ "PCBDDCCreateFETIDPOperators" 2385 /*@ 2386 PCBDDCCreateFETIDPOperators - Create FETI-DP operators 2387 2388 Collective 2389 2390 Input Parameters: 2391 + pc - the BDDC preconditioning context (setup should have been called before) 2392 - fully_redundant - true for a fully redundant set of Lagrange multipliers 2393 2394 Output Parameters: 2395 + fetidp_mat - shell FETI-DP matrix object 2396 - fetidp_pc - shell Dirichlet preconditioner for FETI-DP matrix 2397 2398 Level: developer 2399 2400 Notes: 2401 Currently the only operations provided for FETI-DP matrix are MatMult and MatMultTranspose 2402 2403 .seealso: PCBDDC, PCBDDCMatFETIDPGetRHS, PCBDDCMatFETIDPGetSolution 2404 @*/ 2405 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, PetscBool fully_redundant, Mat *fetidp_mat, PC *fetidp_pc) 2406 { 2407 PetscErrorCode ierr; 2408 2409 PetscFunctionBegin; 2410 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2411 if (pc->setupcalled) { 2412 ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,PetscBool,Mat*,PC*),(pc,fully_redundant,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 2413 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n"); 2414 PetscFunctionReturn(0); 2415 } 2416 /* -------------------------------------------------------------------------- */ 2417 /*MC 2418 PCBDDC - Balancing Domain Decomposition by Constraints. 2419 2420 An implementation of the BDDC preconditioner based on 2421 2422 .vb 2423 [1] C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 2424 [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 2425 [3] J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", http://arxiv.org/abs/0712.3977 2426 [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 2427 .ve 2428 2429 The matrix to be preconditioned (Pmat) must be of type MATIS. 2430 2431 Currently works with MATIS matrices with local matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers. 2432 2433 It also works with unsymmetric and indefinite problems. 2434 2435 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. 2436 2437 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). 2438 2439 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() 2440 Additional information on dofs can be provided by using PCBDDCSetDofsSplitting(), PCBDDCSetDirichletBoundaries(), PCBDDCSetNeumannBoundaries(), and PCBDDCSetPrimalVerticesIS() and their local counterparts. 2441 2442 Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace(). Non-singular modes are retained via SVD. 2443 2444 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. 2445 User defined change of basis can be passed to PCBDDC by using PCBDDCSetChangeOfBasisMat() 2446 2447 The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using a MatPartitioning object. 2448 2449 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. 2450 2451 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. 2452 Deluxe scaling is not supported yet for FETI-DP. 2453 2454 Options Database Keys (some of them, run with -h for a complete list): 2455 2456 . -pc_bddc_use_vertices <true> - use or not vertices in primal space 2457 . -pc_bddc_use_edges <true> - use or not edges in primal space 2458 . -pc_bddc_use_faces <false> - use or not faces in primal space 2459 . -pc_bddc_symmetric <true> - symmetric computation of primal basis functions. Specify false for unsymmetric problems 2460 . -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only) 2461 . -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested 2462 . -pc_bddc_switch_static <false> - switches from M_2 (default) to M_3 operator (see reference article [1]) 2463 . -pc_bddc_levels <0> - maximum number of levels for multilevel 2464 . -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) 2465 . -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) 2466 . -pc_bddc_use_deluxe_scaling <false> - use deluxe scaling 2467 . -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) 2468 . -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) 2469 - -pc_bddc_check_level <0> - set verbosity level of debugging output 2470 2471 Options for Dirichlet, Neumann or coarse solver can be set with 2472 .vb 2473 -pc_bddc_dirichlet_ 2474 -pc_bddc_neumann_ 2475 -pc_bddc_coarse_ 2476 .ve 2477 e.g -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. PCBDDC uses by default KPSPREONLY and PCLU. 2478 2479 When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified as 2480 .vb 2481 -pc_bddc_dirichlet_lN_ 2482 -pc_bddc_neumann_lN_ 2483 -pc_bddc_coarse_lN_ 2484 .ve 2485 Note that level number ranges from the finest (0) to the coarsest (N). 2486 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. 2487 .vb 2488 -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3 2489 .ve 2490 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 2491 2492 Level: intermediate 2493 2494 Developer notes: 2495 2496 Contributed by Stefano Zampini 2497 2498 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 2499 M*/ 2500 2501 #undef __FUNCT__ 2502 #define __FUNCT__ "PCCreate_BDDC" 2503 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 2504 { 2505 PetscErrorCode ierr; 2506 PC_BDDC *pcbddc; 2507 2508 PetscFunctionBegin; 2509 /* Creates the private data structure for this preconditioner and attach it to the PC object. */ 2510 ierr = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr); 2511 pc->data = (void*)pcbddc; 2512 2513 /* create PCIS data structure */ 2514 ierr = PCISCreate(pc);CHKERRQ(ierr); 2515 2516 /* BDDC customization */ 2517 pcbddc->use_local_adj = PETSC_TRUE; 2518 pcbddc->use_vertices = PETSC_TRUE; 2519 pcbddc->use_edges = PETSC_TRUE; 2520 pcbddc->use_faces = PETSC_FALSE; 2521 pcbddc->use_change_of_basis = PETSC_FALSE; 2522 pcbddc->use_change_on_faces = PETSC_FALSE; 2523 pcbddc->switch_static = PETSC_FALSE; 2524 pcbddc->use_nnsp_true = PETSC_FALSE; 2525 pcbddc->use_qr_single = PETSC_FALSE; 2526 pcbddc->symmetric_primal = PETSC_TRUE; 2527 pcbddc->benign_saddle_point = PETSC_FALSE; 2528 pcbddc->benign_have_null = PETSC_FALSE; 2529 pcbddc->vertex_size = 1; 2530 pcbddc->dbg_flag = 0; 2531 /* private */ 2532 pcbddc->local_primal_size = 0; 2533 pcbddc->local_primal_size_cc = 0; 2534 pcbddc->local_primal_ref_node = 0; 2535 pcbddc->local_primal_ref_mult = 0; 2536 pcbddc->n_vertices = 0; 2537 pcbddc->primal_indices_local_idxs = 0; 2538 pcbddc->recompute_topography = PETSC_TRUE; 2539 pcbddc->coarse_size = -1; 2540 pcbddc->new_primal_space = PETSC_FALSE; 2541 pcbddc->new_primal_space_local = PETSC_FALSE; 2542 pcbddc->global_primal_indices = 0; 2543 pcbddc->onearnullspace = 0; 2544 pcbddc->onearnullvecs_state = 0; 2545 pcbddc->user_primal_vertices = 0; 2546 pcbddc->user_primal_vertices_local = 0; 2547 pcbddc->temp_solution = 0; 2548 pcbddc->original_rhs = 0; 2549 pcbddc->local_mat = 0; 2550 pcbddc->ChangeOfBasisMatrix = 0; 2551 pcbddc->user_ChangeOfBasisMatrix = 0; 2552 pcbddc->coarse_vec = 0; 2553 pcbddc->coarse_ksp = 0; 2554 pcbddc->coarse_phi_B = 0; 2555 pcbddc->coarse_phi_D = 0; 2556 pcbddc->coarse_psi_B = 0; 2557 pcbddc->coarse_psi_D = 0; 2558 pcbddc->vec1_P = 0; 2559 pcbddc->vec1_R = 0; 2560 pcbddc->vec2_R = 0; 2561 pcbddc->local_auxmat1 = 0; 2562 pcbddc->local_auxmat2 = 0; 2563 pcbddc->R_to_B = 0; 2564 pcbddc->R_to_D = 0; 2565 pcbddc->ksp_D = 0; 2566 pcbddc->ksp_R = 0; 2567 pcbddc->NeumannBoundaries = 0; 2568 pcbddc->NeumannBoundariesLocal = 0; 2569 pcbddc->DirichletBoundaries = 0; 2570 pcbddc->DirichletBoundariesLocal = 0; 2571 pcbddc->user_provided_isfordofs = PETSC_FALSE; 2572 pcbddc->n_ISForDofs = 0; 2573 pcbddc->n_ISForDofsLocal = 0; 2574 pcbddc->ISForDofs = 0; 2575 pcbddc->ISForDofsLocal = 0; 2576 pcbddc->ConstraintMatrix = 0; 2577 pcbddc->use_exact_dirichlet_trick = PETSC_TRUE; 2578 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 2579 pcbddc->coarse_loc_to_glob = 0; 2580 pcbddc->coarsening_ratio = 8; 2581 pcbddc->coarse_adj_red = 0; 2582 pcbddc->current_level = 0; 2583 pcbddc->max_levels = 0; 2584 pcbddc->use_coarse_estimates = PETSC_FALSE; 2585 pcbddc->coarse_eqs_per_proc = 1; 2586 pcbddc->coarse_subassembling = 0; 2587 pcbddc->detect_disconnected = PETSC_FALSE; 2588 pcbddc->n_local_subs = 0; 2589 pcbddc->local_subs = NULL; 2590 2591 /* benign subspace trick */ 2592 pcbddc->benign_change = 0; 2593 pcbddc->benign_compute_correction = PETSC_TRUE; 2594 pcbddc->benign_vec = 0; 2595 pcbddc->benign_original_mat = 0; 2596 pcbddc->benign_sf = 0; 2597 pcbddc->benign_B0 = 0; 2598 pcbddc->benign_n = 0; 2599 pcbddc->benign_p0 = NULL; 2600 pcbddc->benign_p0_lidx = NULL; 2601 pcbddc->benign_p0_gidx = NULL; 2602 pcbddc->benign_null = PETSC_FALSE; 2603 2604 /* create local graph structure */ 2605 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 2606 pcbddc->graphmaxcount = PETSC_MAX_INT; 2607 2608 /* scaling */ 2609 pcbddc->work_scaling = 0; 2610 pcbddc->use_deluxe_scaling = PETSC_FALSE; 2611 2612 /* sub schurs options */ 2613 pcbddc->sub_schurs_rebuild = PETSC_FALSE; 2614 pcbddc->sub_schurs_layers = -1; 2615 pcbddc->sub_schurs_use_useradj = PETSC_FALSE; 2616 2617 pcbddc->computed_rowadj = PETSC_FALSE; 2618 2619 /* adaptivity */ 2620 pcbddc->adaptive_threshold = 0.0; 2621 pcbddc->adaptive_nmax = 0; 2622 pcbddc->adaptive_nmin = 0; 2623 2624 /* function pointers */ 2625 pc->ops->apply = PCApply_BDDC; 2626 pc->ops->applytranspose = PCApplyTranspose_BDDC; 2627 pc->ops->setup = PCSetUp_BDDC; 2628 pc->ops->destroy = PCDestroy_BDDC; 2629 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 2630 pc->ops->view = PCView_BDDC; 2631 pc->ops->applyrichardson = 0; 2632 pc->ops->applysymmetricleft = 0; 2633 pc->ops->applysymmetricright = 0; 2634 pc->ops->presolve = PCPreSolve_BDDC; 2635 pc->ops->postsolve = PCPostSolve_BDDC; 2636 2637 /* composing function */ 2638 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDiscreteGradient_C",PCBDDCSetDiscreteGradient_BDDC);CHKERRQ(ierr); 2639 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",PCBDDCSetDivergenceMat_BDDC);CHKERRQ(ierr); 2640 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);CHKERRQ(ierr); 2641 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 2642 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",PCBDDCSetPrimalVerticesIS_BDDC);CHKERRQ(ierr); 2643 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 2644 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr); 2645 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr); 2646 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr); 2647 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2648 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2649 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2650 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2651 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2652 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2653 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2654 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2655 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 2656 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr); 2657 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 2658 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 2659 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 2660 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 2661 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",PCPreSolveChangeRHS_BDDC);CHKERRQ(ierr); 2662 PetscFunctionReturn(0); 2663 } 2664 2665