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