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) pcbddc->graphanalyzed = PETSC_FALSE; 1522 computeconstraintsmatrix = PETSC_FALSE; 1523 1524 /* check parameters' compatibility */ 1525 if (!pcbddc->use_deluxe_scaling) pcbddc->deluxe_zerorows = PETSC_FALSE; 1526 pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_threshold > 0.0); 1527 pcbddc->adaptive_userdefined = (PetscBool)(pcbddc->adaptive_selection && pcbddc->adaptive_userdefined); 1528 if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE; 1529 1530 computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling); 1531 if (pcbddc->switch_static) { 1532 PetscBool ismatis; 1533 ierr = PetscObjectTypeCompare((PetscObject)pc->mat,MATIS,&ismatis);CHKERRQ(ierr); 1534 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"When the static switch is one, the iteration matrix should be of type MATIS"); 1535 } 1536 1537 /* activate all connected components if the netflux has been requested */ 1538 if (pcbddc->compute_nonetflux) { 1539 pcbddc->use_vertices = PETSC_TRUE; 1540 pcbddc->use_edges = PETSC_TRUE; 1541 pcbddc->use_faces = PETSC_TRUE; 1542 } 1543 1544 /* Get stdout for dbg */ 1545 if (pcbddc->dbg_flag) { 1546 if (!pcbddc->dbg_viewer) { 1547 pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)); 1548 ierr = PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer);CHKERRQ(ierr); 1549 } 1550 ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1551 } 1552 1553 /* process topology information */ 1554 if (pcbddc->recompute_topography) { 1555 ierr = PCBDDCComputeLocalTopologyInfo(pc);CHKERRQ(ierr); 1556 if (pcbddc->discretegradient) { 1557 ierr = PCBDDCNedelecSupport(pc);CHKERRQ(ierr); 1558 } 1559 } 1560 1561 /* change basis if requested by the user */ 1562 if (pcbddc->user_ChangeOfBasisMatrix) { 1563 /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */ 1564 pcbddc->use_change_of_basis = PETSC_FALSE; 1565 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 1566 } else { 1567 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1568 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1569 pcbddc->local_mat = matis->A; 1570 } 1571 1572 /* detect local disconnected subdomains if requested */ 1573 if (pcbddc->detect_disconnected && pcbddc->recompute_topography) { 1574 PetscInt i; 1575 for (i=0;i<pcbddc->n_local_subs;i++) { 1576 ierr = ISDestroy(&pcbddc->local_subs[i]);CHKERRQ(ierr); 1577 } 1578 ierr = PetscFree(pcbddc->local_subs);CHKERRQ(ierr); 1579 ierr = MatDetectDisconnectedComponents(pcbddc->local_mat,PETSC_FALSE,&pcbddc->n_local_subs,&pcbddc->local_subs);CHKERRQ(ierr); 1580 } 1581 1582 /* 1583 Compute change of basis on local pressures (aka zerodiag dofs) with the benign trick 1584 This should come earlier then PCISSetUp for extracting the correct subdomain matrices 1585 */ 1586 ierr = PCBDDCBenignShellMat(pc,PETSC_TRUE);CHKERRQ(ierr); 1587 if (pcbddc->benign_saddle_point) { 1588 PC_IS* pcis = (PC_IS*)pc->data; 1589 1590 if (pcbddc->user_ChangeOfBasisMatrix || pcbddc->use_change_of_basis || !computesubschurs) pcbddc->benign_change_explicit = PETSC_TRUE; 1591 /* detect local saddle point and change the basis in pcbddc->local_mat (TODO: reuse case) */ 1592 ierr = PCBDDCBenignDetectSaddlePoint(pc,&zerodiag);CHKERRQ(ierr); 1593 /* pop B0 mat from local mat */ 1594 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 1595 /* give pcis a hint to not reuse submatrices during PCISCreate */ 1596 if (pc->flag == SAME_NONZERO_PATTERN && pcis->reusesubmatrices == PETSC_TRUE) { 1597 if (pcbddc->benign_n && (pcbddc->benign_change_explicit || pcbddc->dbg_flag)) { 1598 pcis->reusesubmatrices = PETSC_FALSE; 1599 } else { 1600 pcis->reusesubmatrices = PETSC_TRUE; 1601 } 1602 } else { 1603 pcis->reusesubmatrices = PETSC_FALSE; 1604 } 1605 } 1606 1607 /* propagate relevant information -> TODO remove*/ 1608 #if !defined(PETSC_USE_COMPLEX) /* workaround for reals */ 1609 if (matis->A->symmetric_set) { 1610 ierr = MatSetOption(pcbddc->local_mat,MAT_HERMITIAN,matis->A->symmetric);CHKERRQ(ierr); 1611 } 1612 #endif 1613 if (matis->A->symmetric_set) { 1614 ierr = MatSetOption(pcbddc->local_mat,MAT_SYMMETRIC,matis->A->symmetric);CHKERRQ(ierr); 1615 } 1616 if (matis->A->spd_set) { 1617 ierr = MatSetOption(pcbddc->local_mat,MAT_SPD,matis->A->spd);CHKERRQ(ierr); 1618 } 1619 1620 /* Set up all the "iterative substructuring" common block without computing solvers */ 1621 { 1622 Mat temp_mat; 1623 1624 temp_mat = matis->A; 1625 matis->A = pcbddc->local_mat; 1626 ierr = PCISSetUp(pc,PETSC_FALSE);CHKERRQ(ierr); 1627 pcbddc->local_mat = matis->A; 1628 matis->A = temp_mat; 1629 } 1630 1631 /* Analyze interface */ 1632 if (!pcbddc->graphanalyzed) { 1633 ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr); 1634 computeconstraintsmatrix = PETSC_TRUE; 1635 if (pcbddc->adaptive_selection && !pcbddc->use_deluxe_scaling && !pcbddc->mat_graph->twodim) { 1636 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot compute the adaptive primal space for a problem with 3D edges without deluxe scaling"); 1637 } 1638 if (pcbddc->compute_nonetflux) { 1639 MatNullSpace nnfnnsp; 1640 1641 ierr = PCBDDCComputeNoNetFlux(pc->pmat,pcbddc->divudotp,pcbddc->divudotp_trans,pcbddc->divudotp_vl2l,pcbddc->mat_graph,&nnfnnsp);CHKERRQ(ierr); 1642 /* TODO what if a nearnullspace is already attached? */ 1643 ierr = MatSetNearNullSpace(pc->pmat,nnfnnsp);CHKERRQ(ierr); 1644 ierr = MatNullSpaceDestroy(&nnfnnsp);CHKERRQ(ierr); 1645 } 1646 } 1647 1648 /* check existence of a divergence free extension, i.e. 1649 b(v_I,p_0) = 0 for all v_I (raise error if not). 1650 Also, check that PCBDDCBenignGetOrSetP0 works */ 1651 #if defined(PETSC_USE_DEBUG) 1652 if (pcbddc->benign_saddle_point) { 1653 ierr = PCBDDCBenignCheck(pc,zerodiag);CHKERRQ(ierr); 1654 } 1655 #endif 1656 ierr = ISDestroy(&zerodiag);CHKERRQ(ierr); 1657 1658 /* Setup local dirichlet solver ksp_D and sub_schurs solvers */ 1659 if (computesubschurs && pcbddc->recompute_topography) { 1660 ierr = PCBDDCInitSubSchurs(pc);CHKERRQ(ierr); 1661 } 1662 /* SetUp Scaling operator (scaling matrices could be needed in SubSchursSetUp)*/ 1663 if (!pcbddc->use_deluxe_scaling) { 1664 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1665 } 1666 1667 /* finish setup solvers and do adaptive selection of constraints */ 1668 sub_schurs = pcbddc->sub_schurs; 1669 if (sub_schurs && sub_schurs->schur_explicit) { 1670 if (computesubschurs) { 1671 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1672 } 1673 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1674 } else { 1675 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1676 if (computesubschurs) { 1677 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1678 } 1679 } 1680 if (pcbddc->adaptive_selection) { 1681 ierr = PCBDDCAdaptiveSelection(pc);CHKERRQ(ierr); 1682 computeconstraintsmatrix = PETSC_TRUE; 1683 } 1684 1685 /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */ 1686 new_nearnullspace_provided = PETSC_FALSE; 1687 ierr = MatGetNearNullSpace(pc->pmat,&nearnullspace);CHKERRQ(ierr); 1688 if (pcbddc->onearnullspace) { /* already used nearnullspace */ 1689 if (!nearnullspace) { /* near null space attached to mat has been destroyed */ 1690 new_nearnullspace_provided = PETSC_TRUE; 1691 } else { 1692 /* determine if the two nullspaces are different (should be lightweight) */ 1693 if (nearnullspace != pcbddc->onearnullspace) { 1694 new_nearnullspace_provided = PETSC_TRUE; 1695 } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */ 1696 PetscInt i; 1697 const Vec *nearnullvecs; 1698 PetscObjectState state; 1699 PetscInt nnsp_size; 1700 ierr = MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 1701 for (i=0;i<nnsp_size;i++) { 1702 ierr = PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);CHKERRQ(ierr); 1703 if (pcbddc->onearnullvecs_state[i] != state) { 1704 new_nearnullspace_provided = PETSC_TRUE; 1705 break; 1706 } 1707 } 1708 } 1709 } 1710 } else { 1711 if (!nearnullspace) { /* both nearnullspaces are null */ 1712 new_nearnullspace_provided = PETSC_FALSE; 1713 } else { /* nearnullspace attached later */ 1714 new_nearnullspace_provided = PETSC_TRUE; 1715 } 1716 } 1717 1718 /* Setup constraints and related work vectors */ 1719 /* reset primal space flags */ 1720 pcbddc->new_primal_space = PETSC_FALSE; 1721 pcbddc->new_primal_space_local = PETSC_FALSE; 1722 if (computeconstraintsmatrix || new_nearnullspace_provided) { 1723 /* It also sets the primal space flags */ 1724 ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr); 1725 /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */ 1726 ierr = PCBDDCSetUpLocalWorkVectors(pc);CHKERRQ(ierr); 1727 } 1728 1729 if (pcbddc->use_change_of_basis) { 1730 PC_IS *pcis = (PC_IS*)(pc->data); 1731 1732 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 1733 if (pcbddc->benign_change) { 1734 ierr = MatDestroy(&pcbddc->benign_B0);CHKERRQ(ierr); 1735 /* pop B0 from pcbddc->local_mat */ 1736 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 1737 } 1738 /* get submatrices */ 1739 ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr); 1740 ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr); 1741 ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr); 1742 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);CHKERRQ(ierr); 1743 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);CHKERRQ(ierr); 1744 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);CHKERRQ(ierr); 1745 /* set flag in pcis to not reuse submatrices during PCISCreate */ 1746 pcis->reusesubmatrices = PETSC_FALSE; 1747 } else if (!pcbddc->user_ChangeOfBasisMatrix && !pcbddc->benign_change) { 1748 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1749 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1750 pcbddc->local_mat = matis->A; 1751 } 1752 /* SetUp coarse and local Neumann solvers */ 1753 ierr = PCBDDCSetUpSolvers(pc);CHKERRQ(ierr); 1754 /* SetUp Scaling operator */ 1755 if (pcbddc->use_deluxe_scaling) { 1756 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1757 } 1758 1759 /* mark topography as done */ 1760 pcbddc->recompute_topography = PETSC_FALSE; 1761 1762 /* wrap pcis->A_IB and pcis->A_BI if we did not change explicitly the variables on the pressures */ 1763 ierr = PCBDDCBenignShellMat(pc,PETSC_FALSE);CHKERRQ(ierr); 1764 1765 if (pcbddc->dbg_flag) { 1766 ierr = PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1767 } 1768 PetscFunctionReturn(0); 1769 } 1770 1771 /* 1772 PCApply_BDDC - Applies the BDDC operator to a vector. 1773 1774 Input Parameters: 1775 + pc - the preconditioner context 1776 - r - input vector (global) 1777 1778 Output Parameter: 1779 . z - output vector (global) 1780 1781 Application Interface Routine: PCApply() 1782 */ 1783 #undef __FUNCT__ 1784 #define __FUNCT__ "PCApply_BDDC" 1785 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z) 1786 { 1787 PC_IS *pcis = (PC_IS*)(pc->data); 1788 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1789 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1790 PetscErrorCode ierr; 1791 const PetscScalar one = 1.0; 1792 const PetscScalar m_one = -1.0; 1793 const PetscScalar zero = 0.0; 1794 1795 /* This code is similar to that provided in nn.c for PCNN 1796 NN interface preconditioner changed to BDDC 1797 Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */ 1798 1799 PetscFunctionBegin; 1800 if (pcbddc->ChangeOfBasisMatrix) { 1801 Vec swap; 1802 1803 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);CHKERRQ(ierr); 1804 swap = pcbddc->work_change; 1805 pcbddc->work_change = r; 1806 r = swap; 1807 /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */ 1808 if (pcbddc->benign_apply_coarse_only && pcbddc->exact_dirichlet_trick_app && pcbddc->change_interior) { 1809 ierr = VecCopy(r,pcis->vec1_global);CHKERRQ(ierr); 1810 ierr = VecLockPush(pcis->vec1_global);CHKERRQ(ierr); 1811 } 1812 } 1813 if (pcbddc->benign_have_null) { /* get p0 from r */ 1814 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 1815 } 1816 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 1817 ierr = VecCopy(r,z);CHKERRQ(ierr); 1818 /* First Dirichlet solve */ 1819 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1820 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1821 /* 1822 Assembling right hand side for BDDC operator 1823 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1824 - pcis->vec1_B the interface part of the global vector z 1825 */ 1826 if (n_D) { 1827 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1828 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1829 if (pcbddc->switch_static) { 1830 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1831 1832 ierr = VecSet(pcis->vec1_N,0.);CHKERRQ(ierr); 1833 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1834 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1835 if (!pcbddc->switch_static_change) { 1836 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1837 } else { 1838 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1839 ierr = MatMult(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1840 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1841 } 1842 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1843 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1844 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1845 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1846 } else { 1847 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1848 } 1849 } else { 1850 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1851 } 1852 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1853 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1854 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1855 } else { 1856 if (!pcbddc->benign_apply_coarse_only) { 1857 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1858 } 1859 } 1860 1861 /* Apply interface preconditioner 1862 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1863 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);CHKERRQ(ierr); 1864 1865 /* Apply transpose of partition of unity operator */ 1866 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1867 1868 /* Second Dirichlet solve and assembling of output */ 1869 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1870 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1871 if (n_B) { 1872 if (pcbddc->switch_static) { 1873 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1874 1875 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1876 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1877 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1878 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1879 if (!pcbddc->switch_static_change) { 1880 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1881 } else { 1882 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1883 ierr = MatMult(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1884 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1885 } 1886 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1887 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1888 } else { 1889 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1890 } 1891 } else if (pcbddc->switch_static) { /* n_B is zero */ 1892 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1893 1894 if (!pcbddc->switch_static_change) { 1895 ierr = MatMult(matis->A,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 1896 } else { 1897 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);CHKERRQ(ierr); 1898 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1899 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);CHKERRQ(ierr); 1900 } 1901 } 1902 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1903 1904 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 1905 if (pcbddc->switch_static) { 1906 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 1907 } else { 1908 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 1909 } 1910 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1911 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1912 } else { 1913 if (pcbddc->switch_static) { 1914 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 1915 } else { 1916 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1917 } 1918 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1919 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1920 } 1921 if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */ 1922 if (pcbddc->benign_apply_coarse_only) { 1923 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 1924 } 1925 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 1926 } 1927 1928 if (pcbddc->ChangeOfBasisMatrix) { 1929 pcbddc->work_change = r; 1930 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 1931 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 1932 } 1933 PetscFunctionReturn(0); 1934 } 1935 1936 /* 1937 PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector. 1938 1939 Input Parameters: 1940 + pc - the preconditioner context 1941 - r - input vector (global) 1942 1943 Output Parameter: 1944 . z - output vector (global) 1945 1946 Application Interface Routine: PCApplyTranspose() 1947 */ 1948 #undef __FUNCT__ 1949 #define __FUNCT__ "PCApplyTranspose_BDDC" 1950 PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z) 1951 { 1952 PC_IS *pcis = (PC_IS*)(pc->data); 1953 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1954 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1955 PetscErrorCode ierr; 1956 const PetscScalar one = 1.0; 1957 const PetscScalar m_one = -1.0; 1958 const PetscScalar zero = 0.0; 1959 1960 PetscFunctionBegin; 1961 if (pcbddc->ChangeOfBasisMatrix) { 1962 Vec swap; 1963 1964 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);CHKERRQ(ierr); 1965 swap = pcbddc->work_change; 1966 pcbddc->work_change = r; 1967 r = swap; 1968 /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */ 1969 if (pcbddc->benign_apply_coarse_only && pcbddc->exact_dirichlet_trick_app && pcbddc->change_interior) { 1970 ierr = VecCopy(r,pcis->vec1_global);CHKERRQ(ierr); 1971 ierr = VecLockPush(pcis->vec1_global);CHKERRQ(ierr); 1972 } 1973 } 1974 if (pcbddc->benign_have_null) { /* get p0 from r */ 1975 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 1976 } 1977 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 1978 ierr = VecCopy(r,z);CHKERRQ(ierr); 1979 /* First Dirichlet solve */ 1980 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1981 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1982 /* 1983 Assembling right hand side for BDDC operator 1984 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1985 - pcis->vec1_B the interface part of the global vector z 1986 */ 1987 if (n_D) { 1988 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1989 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1990 if (pcbddc->switch_static) { 1991 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1992 1993 ierr = VecSet(pcis->vec1_N,0.);CHKERRQ(ierr); 1994 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1995 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1996 if (!pcbddc->switch_static_change) { 1997 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1998 } else { 1999 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2000 ierr = MatMultTranspose(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 2001 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2002 } 2003 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2004 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2005 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2006 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2007 } else { 2008 ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 2009 } 2010 } else { 2011 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 2012 } 2013 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2014 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2015 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 2016 } else { 2017 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 2018 } 2019 2020 /* Apply interface preconditioner 2021 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 2022 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr); 2023 2024 /* Apply transpose of partition of unity operator */ 2025 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 2026 2027 /* Second Dirichlet solve and assembling of output */ 2028 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2029 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2030 if (n_B) { 2031 if (pcbddc->switch_static) { 2032 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 2033 2034 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2035 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2036 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2037 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2038 if (!pcbddc->switch_static_change) { 2039 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2040 } else { 2041 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2042 ierr = MatMultTranspose(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 2043 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2044 } 2045 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2046 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2047 } else { 2048 ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 2049 } 2050 } else if (pcbddc->switch_static) { /* n_B is zero */ 2051 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 2052 2053 if (!pcbddc->switch_static_change) { 2054 ierr = MatMultTranspose(matis->A,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 2055 } else { 2056 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);CHKERRQ(ierr); 2057 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2058 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);CHKERRQ(ierr); 2059 } 2060 } 2061 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 2062 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 2063 if (pcbddc->switch_static) { 2064 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 2065 } else { 2066 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 2067 } 2068 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2069 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2070 } else { 2071 if (pcbddc->switch_static) { 2072 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 2073 } else { 2074 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 2075 } 2076 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2077 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2078 } 2079 if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */ 2080 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 2081 } 2082 if (pcbddc->ChangeOfBasisMatrix) { 2083 pcbddc->work_change = r; 2084 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 2085 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 2086 } 2087 PetscFunctionReturn(0); 2088 } 2089 2090 #undef __FUNCT__ 2091 #define __FUNCT__ "PCDestroy_BDDC" 2092 PetscErrorCode PCDestroy_BDDC(PC pc) 2093 { 2094 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 2095 PetscErrorCode ierr; 2096 2097 PetscFunctionBegin; 2098 /* free BDDC custom data */ 2099 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 2100 /* destroy objects related to topography */ 2101 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 2102 /* free allocated graph structure */ 2103 ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr); 2104 /* free allocated sub schurs structure */ 2105 ierr = PetscFree(pcbddc->sub_schurs);CHKERRQ(ierr); 2106 /* destroy objects for scaling operator */ 2107 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 2108 ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr); 2109 /* free solvers stuff */ 2110 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 2111 /* free global vectors needed in presolve */ 2112 ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr); 2113 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 2114 /* free data created by PCIS */ 2115 ierr = PCISDestroy(pc);CHKERRQ(ierr); 2116 /* remove functions */ 2117 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDiscreteGradient_C",NULL);CHKERRQ(ierr); 2118 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",NULL);CHKERRQ(ierr); 2119 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);CHKERRQ(ierr); 2120 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 2121 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",NULL);CHKERRQ(ierr); 2122 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 2123 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr); 2124 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr); 2125 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr); 2126 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2127 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2128 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2129 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2130 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2131 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2132 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2133 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2134 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 2135 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr); 2136 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 2137 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 2138 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 2139 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 2140 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",NULL);CHKERRQ(ierr); 2141 /* Free the private data structure */ 2142 ierr = PetscFree(pc->data);CHKERRQ(ierr); 2143 PetscFunctionReturn(0); 2144 } 2145 2146 #undef __FUNCT__ 2147 #define __FUNCT__ "PCPreSolveChangeRHS_BDDC" 2148 static PetscErrorCode PCPreSolveChangeRHS_BDDC(PC pc, PetscBool* change) 2149 { 2150 PetscFunctionBegin; 2151 *change = PETSC_TRUE; 2152 PetscFunctionReturn(0); 2153 } 2154 2155 #undef __FUNCT__ 2156 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC" 2157 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 2158 { 2159 FETIDPMat_ctx mat_ctx; 2160 Vec work; 2161 PC_IS* pcis; 2162 PC_BDDC* pcbddc; 2163 PetscErrorCode ierr; 2164 2165 PetscFunctionBegin; 2166 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2167 pcis = (PC_IS*)mat_ctx->pc->data; 2168 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2169 2170 /* 2171 change of basis for physical rhs if needed 2172 It also changes the rhs in case of dirichlet boundaries 2173 */ 2174 if (!pcbddc->original_rhs) { 2175 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 2176 } 2177 ierr = VecCopy(standard_rhs,pcbddc->original_rhs);CHKERRQ(ierr); 2178 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,pcbddc->original_rhs,NULL);CHKERRQ(ierr); 2179 if (pcbddc->ChangeOfBasisMatrix) { 2180 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,pcbddc->original_rhs,pcbddc->work_change);CHKERRQ(ierr); 2181 work = pcbddc->work_change; 2182 } else { 2183 work = pcbddc->original_rhs; 2184 } 2185 /* store vectors for computation of fetidp final solution */ 2186 ierr = VecScatterBegin(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2187 ierr = VecScatterEnd(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2188 /* scale rhs since it should be unassembled */ 2189 /* TODO use counter scaling? (also below) */ 2190 ierr = VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2191 ierr = VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2192 /* Apply partition of unity */ 2193 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2194 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2195 if (!pcbddc->switch_static) { 2196 /* compute partially subassembled Schur complement right-hand side */ 2197 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2198 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 2199 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 2200 ierr = VecSet(work,0.0);CHKERRQ(ierr); 2201 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2202 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2203 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2204 ierr = VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2205 ierr = VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2206 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2207 } 2208 /* BDDC rhs */ 2209 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 2210 if (pcbddc->switch_static) { 2211 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2212 } 2213 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2214 /* apply BDDC */ 2215 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2216 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2217 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 2218 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 2219 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 2220 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2221 ierr = VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2222 PetscFunctionReturn(0); 2223 } 2224 2225 #undef __FUNCT__ 2226 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS" 2227 /*@ 2228 PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETI-DP linear system using the physical right-hand side 2229 2230 Collective 2231 2232 Input Parameters: 2233 + fetidp_mat - the FETI-DP matrix object obtained by a call to PCBDDCCreateFETIDPOperators 2234 - standard_rhs - the right-hand side of the original linear system 2235 2236 Output Parameters: 2237 . fetidp_flux_rhs - the right-hand side for the FETI-DP linear system 2238 2239 Level: developer 2240 2241 Notes: 2242 2243 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetSolution 2244 @*/ 2245 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 2246 { 2247 FETIDPMat_ctx mat_ctx; 2248 PetscErrorCode ierr; 2249 2250 PetscFunctionBegin; 2251 PetscValidHeaderSpecific(fetidp_mat,MAT_CLASSID,1); 2252 PetscValidHeaderSpecific(standard_rhs,VEC_CLASSID,2); 2253 PetscValidHeaderSpecific(fetidp_flux_rhs,VEC_CLASSID,3); 2254 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2255 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 2256 PetscFunctionReturn(0); 2257 } 2258 2259 #undef __FUNCT__ 2260 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC" 2261 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2262 { 2263 FETIDPMat_ctx mat_ctx; 2264 PC_IS* pcis; 2265 PC_BDDC* pcbddc; 2266 PetscErrorCode ierr; 2267 2268 PetscFunctionBegin; 2269 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2270 pcis = (PC_IS*)mat_ctx->pc->data; 2271 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2272 2273 /* apply B_delta^T */ 2274 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2275 ierr = VecScatterEnd(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2276 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 2277 /* compute rhs for BDDC application */ 2278 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2279 if (pcbddc->switch_static) { 2280 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2281 } 2282 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2283 /* apply BDDC */ 2284 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2285 /* put values into standard global vector */ 2286 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2287 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2288 if (!pcbddc->switch_static) { 2289 /* compute values into the interior if solved for the partially subassembled Schur complement */ 2290 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 2291 ierr = VecAYPX(pcis->vec1_D,-1.0,mat_ctx->temp_solution_D);CHKERRQ(ierr); 2292 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); 2293 } 2294 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2295 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2296 /* add p0 solution to final solution */ 2297 ierr = PCBDDCBenignGetOrSetP0(mat_ctx->pc,standard_sol,PETSC_FALSE);CHKERRQ(ierr); 2298 if (pcbddc->ChangeOfBasisMatrix) { 2299 Vec v2; 2300 ierr = VecDuplicate(standard_sol,&v2);CHKERRQ(ierr); 2301 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,standard_sol,v2);CHKERRQ(ierr); 2302 ierr = VecCopy(v2,standard_sol);CHKERRQ(ierr); 2303 ierr = VecDestroy(&v2);CHKERRQ(ierr); 2304 } 2305 ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr); 2306 PetscFunctionReturn(0); 2307 } 2308 2309 #undef __FUNCT__ 2310 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution" 2311 /*@ 2312 PCBDDCMatFETIDPGetSolution - Compute the physical solution using the solution of the FETI-DP linear system 2313 2314 Collective 2315 2316 Input Parameters: 2317 + fetidp_mat - the FETI-DP matrix obtained by a call to PCBDDCCreateFETIDPOperators 2318 - fetidp_flux_sol - the solution of the FETI-DP linear system 2319 2320 Output Parameters: 2321 . standard_sol - the solution defined on the physical domain 2322 2323 Level: developer 2324 2325 Notes: 2326 2327 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetRHS 2328 @*/ 2329 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2330 { 2331 FETIDPMat_ctx mat_ctx; 2332 PetscErrorCode ierr; 2333 2334 PetscFunctionBegin; 2335 PetscValidHeaderSpecific(fetidp_mat,MAT_CLASSID,1); 2336 PetscValidHeaderSpecific(fetidp_flux_sol,VEC_CLASSID,2); 2337 PetscValidHeaderSpecific(standard_sol,VEC_CLASSID,3); 2338 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2339 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 2340 PetscFunctionReturn(0); 2341 } 2342 2343 extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec); 2344 extern PetscErrorCode FETIDPMatMultTranspose(Mat,Vec,Vec); 2345 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat); 2346 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec); 2347 extern PetscErrorCode FETIDPPCApplyTranspose(PC,Vec,Vec); 2348 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC); 2349 2350 #undef __FUNCT__ 2351 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC" 2352 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, PetscBool fully_redundant, Mat *fetidp_mat, PC *fetidp_pc) 2353 { 2354 2355 FETIDPMat_ctx fetidpmat_ctx; 2356 Mat newmat; 2357 FETIDPPC_ctx fetidppc_ctx; 2358 PC newpc; 2359 MPI_Comm comm; 2360 PetscErrorCode ierr; 2361 2362 PetscFunctionBegin; 2363 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 2364 /* FETIDP linear matrix */ 2365 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 2366 fetidpmat_ctx->fully_redundant = fully_redundant; 2367 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 2368 ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 2369 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 2370 ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr); 2371 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 2372 ierr = MatSetUp(newmat);CHKERRQ(ierr); 2373 /* FETIDP preconditioner */ 2374 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 2375 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 2376 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 2377 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 2378 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 2379 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 2380 ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 2381 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 2382 ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr); 2383 ierr = PCSetUp(newpc);CHKERRQ(ierr); 2384 /* return pointers for objects created */ 2385 *fetidp_mat=newmat; 2386 *fetidp_pc=newpc; 2387 PetscFunctionReturn(0); 2388 } 2389 2390 #undef __FUNCT__ 2391 #define __FUNCT__ "PCBDDCCreateFETIDPOperators" 2392 /*@ 2393 PCBDDCCreateFETIDPOperators - Create FETI-DP operators 2394 2395 Collective 2396 2397 Input Parameters: 2398 + pc - the BDDC preconditioning context (setup should have been called before) 2399 - fully_redundant - true for a fully redundant set of Lagrange multipliers 2400 2401 Output Parameters: 2402 + fetidp_mat - shell FETI-DP matrix object 2403 - fetidp_pc - shell Dirichlet preconditioner for FETI-DP matrix 2404 2405 Level: developer 2406 2407 Notes: 2408 Currently the only operations provided for FETI-DP matrix are MatMult and MatMultTranspose 2409 2410 .seealso: PCBDDC, PCBDDCMatFETIDPGetRHS, PCBDDCMatFETIDPGetSolution 2411 @*/ 2412 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, PetscBool fully_redundant, Mat *fetidp_mat, PC *fetidp_pc) 2413 { 2414 PetscErrorCode ierr; 2415 2416 PetscFunctionBegin; 2417 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2418 if (pc->setupcalled) { 2419 ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,PetscBool,Mat*,PC*),(pc,fully_redundant,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 2420 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n"); 2421 PetscFunctionReturn(0); 2422 } 2423 /* -------------------------------------------------------------------------- */ 2424 /*MC 2425 PCBDDC - Balancing Domain Decomposition by Constraints. 2426 2427 An implementation of the BDDC preconditioner based on 2428 2429 .vb 2430 [1] C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 2431 [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 2432 [3] J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", http://arxiv.org/abs/0712.3977 2433 [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 2434 .ve 2435 2436 The matrix to be preconditioned (Pmat) must be of type MATIS. 2437 2438 Currently works with MATIS matrices with local matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers. 2439 2440 It also works with unsymmetric and indefinite problems. 2441 2442 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. 2443 2444 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). 2445 2446 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() 2447 Additional information on dofs can be provided by using PCBDDCSetDofsSplitting(), PCBDDCSetDirichletBoundaries(), PCBDDCSetNeumannBoundaries(), and PCBDDCSetPrimalVerticesIS() and their local counterparts. 2448 2449 Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace(). Non-singular modes are retained via SVD. 2450 2451 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. 2452 User defined change of basis can be passed to PCBDDC by using PCBDDCSetChangeOfBasisMat() 2453 2454 The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using a MatPartitioning object. 2455 2456 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. 2457 2458 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. 2459 Deluxe scaling is not supported yet for FETI-DP. 2460 2461 Options Database Keys (some of them, run with -h for a complete list): 2462 2463 . -pc_bddc_use_vertices <true> - use or not vertices in primal space 2464 . -pc_bddc_use_edges <true> - use or not edges in primal space 2465 . -pc_bddc_use_faces <false> - use or not faces in primal space 2466 . -pc_bddc_symmetric <true> - symmetric computation of primal basis functions. Specify false for unsymmetric problems 2467 . -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only) 2468 . -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested 2469 . -pc_bddc_switch_static <false> - switches from M_2 (default) to M_3 operator (see reference article [1]) 2470 . -pc_bddc_levels <0> - maximum number of levels for multilevel 2471 . -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) 2472 . -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) 2473 . -pc_bddc_use_deluxe_scaling <false> - use deluxe scaling 2474 . -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) 2475 . -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) 2476 - -pc_bddc_check_level <0> - set verbosity level of debugging output 2477 2478 Options for Dirichlet, Neumann or coarse solver can be set with 2479 .vb 2480 -pc_bddc_dirichlet_ 2481 -pc_bddc_neumann_ 2482 -pc_bddc_coarse_ 2483 .ve 2484 e.g -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. PCBDDC uses by default KPSPREONLY and PCLU. 2485 2486 When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified as 2487 .vb 2488 -pc_bddc_dirichlet_lN_ 2489 -pc_bddc_neumann_lN_ 2490 -pc_bddc_coarse_lN_ 2491 .ve 2492 Note that level number ranges from the finest (0) to the coarsest (N). 2493 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. 2494 .vb 2495 -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3 2496 .ve 2497 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 2498 2499 Level: intermediate 2500 2501 Developer notes: 2502 2503 Contributed by Stefano Zampini 2504 2505 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 2506 M*/ 2507 2508 #undef __FUNCT__ 2509 #define __FUNCT__ "PCCreate_BDDC" 2510 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 2511 { 2512 PetscErrorCode ierr; 2513 PC_BDDC *pcbddc; 2514 2515 PetscFunctionBegin; 2516 /* Creates the private data structure for this preconditioner and attach it to the PC object. */ 2517 ierr = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr); 2518 pc->data = (void*)pcbddc; 2519 2520 /* create PCIS data structure */ 2521 ierr = PCISCreate(pc);CHKERRQ(ierr); 2522 2523 /* BDDC customization */ 2524 pcbddc->use_local_adj = PETSC_TRUE; 2525 pcbddc->use_vertices = PETSC_TRUE; 2526 pcbddc->use_edges = PETSC_TRUE; 2527 pcbddc->use_faces = PETSC_FALSE; 2528 pcbddc->use_change_of_basis = PETSC_FALSE; 2529 pcbddc->use_change_on_faces = PETSC_FALSE; 2530 pcbddc->switch_static = PETSC_FALSE; 2531 pcbddc->use_nnsp_true = PETSC_FALSE; 2532 pcbddc->use_qr_single = PETSC_FALSE; 2533 pcbddc->symmetric_primal = PETSC_TRUE; 2534 pcbddc->benign_saddle_point = PETSC_FALSE; 2535 pcbddc->benign_have_null = PETSC_FALSE; 2536 pcbddc->vertex_size = 1; 2537 pcbddc->dbg_flag = 0; 2538 /* private */ 2539 pcbddc->local_primal_size = 0; 2540 pcbddc->local_primal_size_cc = 0; 2541 pcbddc->local_primal_ref_node = 0; 2542 pcbddc->local_primal_ref_mult = 0; 2543 pcbddc->n_vertices = 0; 2544 pcbddc->primal_indices_local_idxs = 0; 2545 pcbddc->recompute_topography = PETSC_TRUE; 2546 pcbddc->coarse_size = -1; 2547 pcbddc->new_primal_space = PETSC_FALSE; 2548 pcbddc->new_primal_space_local = PETSC_FALSE; 2549 pcbddc->global_primal_indices = 0; 2550 pcbddc->onearnullspace = 0; 2551 pcbddc->onearnullvecs_state = 0; 2552 pcbddc->user_primal_vertices = 0; 2553 pcbddc->user_primal_vertices_local = 0; 2554 pcbddc->temp_solution = 0; 2555 pcbddc->original_rhs = 0; 2556 pcbddc->local_mat = 0; 2557 pcbddc->ChangeOfBasisMatrix = 0; 2558 pcbddc->user_ChangeOfBasisMatrix = 0; 2559 pcbddc->coarse_vec = 0; 2560 pcbddc->coarse_ksp = 0; 2561 pcbddc->coarse_phi_B = 0; 2562 pcbddc->coarse_phi_D = 0; 2563 pcbddc->coarse_psi_B = 0; 2564 pcbddc->coarse_psi_D = 0; 2565 pcbddc->vec1_P = 0; 2566 pcbddc->vec1_R = 0; 2567 pcbddc->vec2_R = 0; 2568 pcbddc->local_auxmat1 = 0; 2569 pcbddc->local_auxmat2 = 0; 2570 pcbddc->R_to_B = 0; 2571 pcbddc->R_to_D = 0; 2572 pcbddc->ksp_D = 0; 2573 pcbddc->ksp_R = 0; 2574 pcbddc->NeumannBoundaries = 0; 2575 pcbddc->NeumannBoundariesLocal = 0; 2576 pcbddc->DirichletBoundaries = 0; 2577 pcbddc->DirichletBoundariesLocal = 0; 2578 pcbddc->user_provided_isfordofs = PETSC_FALSE; 2579 pcbddc->n_ISForDofs = 0; 2580 pcbddc->n_ISForDofsLocal = 0; 2581 pcbddc->ISForDofs = 0; 2582 pcbddc->ISForDofsLocal = 0; 2583 pcbddc->ConstraintMatrix = 0; 2584 pcbddc->use_exact_dirichlet_trick = PETSC_TRUE; 2585 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 2586 pcbddc->coarse_loc_to_glob = 0; 2587 pcbddc->coarsening_ratio = 8; 2588 pcbddc->coarse_adj_red = 0; 2589 pcbddc->current_level = 0; 2590 pcbddc->max_levels = 0; 2591 pcbddc->use_coarse_estimates = PETSC_FALSE; 2592 pcbddc->coarse_eqs_per_proc = 1; 2593 pcbddc->coarse_subassembling = 0; 2594 pcbddc->detect_disconnected = PETSC_FALSE; 2595 pcbddc->n_local_subs = 0; 2596 pcbddc->local_subs = NULL; 2597 2598 /* benign subspace trick */ 2599 pcbddc->benign_change = 0; 2600 pcbddc->benign_compute_correction = PETSC_TRUE; 2601 pcbddc->benign_vec = 0; 2602 pcbddc->benign_original_mat = 0; 2603 pcbddc->benign_sf = 0; 2604 pcbddc->benign_B0 = 0; 2605 pcbddc->benign_n = 0; 2606 pcbddc->benign_p0 = NULL; 2607 pcbddc->benign_p0_lidx = NULL; 2608 pcbddc->benign_p0_gidx = NULL; 2609 pcbddc->benign_null = PETSC_FALSE; 2610 2611 /* Nedelec */ 2612 pcbddc->nedfield = -1; 2613 pcbddc->nedglobal = PETSC_TRUE; 2614 2615 /* create local graph structure */ 2616 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 2617 pcbddc->graphmaxcount = PETSC_MAX_INT; 2618 2619 /* scaling */ 2620 pcbddc->work_scaling = 0; 2621 pcbddc->use_deluxe_scaling = PETSC_FALSE; 2622 2623 /* sub schurs options */ 2624 pcbddc->sub_schurs_rebuild = PETSC_FALSE; 2625 pcbddc->sub_schurs_layers = -1; 2626 pcbddc->sub_schurs_use_useradj = PETSC_FALSE; 2627 2628 pcbddc->computed_rowadj = PETSC_FALSE; 2629 2630 /* adaptivity */ 2631 pcbddc->adaptive_threshold = 0.0; 2632 pcbddc->adaptive_nmax = 0; 2633 pcbddc->adaptive_nmin = 0; 2634 2635 /* function pointers */ 2636 pc->ops->apply = PCApply_BDDC; 2637 pc->ops->applytranspose = PCApplyTranspose_BDDC; 2638 pc->ops->setup = PCSetUp_BDDC; 2639 pc->ops->destroy = PCDestroy_BDDC; 2640 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 2641 pc->ops->view = PCView_BDDC; 2642 pc->ops->applyrichardson = 0; 2643 pc->ops->applysymmetricleft = 0; 2644 pc->ops->applysymmetricright = 0; 2645 pc->ops->presolve = PCPreSolve_BDDC; 2646 pc->ops->postsolve = PCPostSolve_BDDC; 2647 2648 /* composing function */ 2649 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDiscreteGradient_C",PCBDDCSetDiscreteGradient_BDDC);CHKERRQ(ierr); 2650 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",PCBDDCSetDivergenceMat_BDDC);CHKERRQ(ierr); 2651 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);CHKERRQ(ierr); 2652 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 2653 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",PCBDDCSetPrimalVerticesIS_BDDC);CHKERRQ(ierr); 2654 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 2655 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr); 2656 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr); 2657 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr); 2658 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2659 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2660 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2661 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2662 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2663 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2664 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2665 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2666 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 2667 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr); 2668 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 2669 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 2670 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 2671 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 2672 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",PCPreSolveChangeRHS_BDDC);CHKERRQ(ierr); 2673 PetscFunctionReturn(0); 2674 } 2675 2676