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