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