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