18 #include "libmesh/libmesh_common.h" 19 #include "libmesh/petsc_macro.h" 21 #ifdef LIBMESH_HAVE_PETSC 22 #if !PETSC_VERSION_LESS_THAN(3,7,3) 23 #if defined(LIBMESH_ENABLE_AMR) && defined(LIBMESH_HAVE_METAPHYSICL) 25 #include "libmesh/ignore_warnings.h" 27 #if PETSC_VERSION_LESS_THAN(3,12,0) 28 #include <petsc/private/dmimpl.h> 30 #include <petscdmshell.h> 31 #include "libmesh/restore_warnings.h" 33 #include "libmesh/petsc_dm_wrapper.h" 34 #include "libmesh/libmesh_logging.h" 35 #include "libmesh/system.h" 36 #include "libmesh/mesh.h" 37 #include "libmesh/mesh_base.h" 38 #include "libmesh/mesh_refinement.h" 39 #include "libmesh/mesh_tools.h" 40 #include "libmesh/partitioner.h" 41 #include "libmesh/dof_map.h" 42 #include "libmesh/elem.h" 43 #include "libmesh/petsc_matrix.h" 59 #if PETSC_VERSION_LESS_THAN(3,9,0) 90 #if PETSC_VERSION_LESS_THAN(3,12,0) 94 PetscErrorCode (*coarsen)(DM,MPI_Comm,DM*) =
nullptr;
101 #if PETSC_VERSION_LESS_THAN(3,12,0) 105 PetscErrorCode (*refine)(DM,MPI_Comm,DM*) =
nullptr;
113 #if PETSC_VERSION_LESS_THAN(3,12,0) 114 if (dm->ops->createinterpolation)
117 PetscErrorCode (*interp)(DM,DM,Mat*,Vec*) =
nullptr;
124 #if PETSC_VERSION_LESS_THAN(3,12,0) 125 if (dm->ops->createrestriction)
128 PetscErrorCode (*createrestriction)(DM,DM,Mat*) =
nullptr;
130 if (createrestriction)
135 #if PETSC_VERSION_LESS_THAN(3,12,0) 136 if (dm->ops->createsubdm)
139 PetscErrorCode (*createsubdm)(DM,PetscInt,
const PetscInt[],IS*,DM*) =
nullptr;
147 #if PETSC_VERSION_LESS_THAN(3,11,0) 150 #elif PETSC_RELEASE_LESS_THAN(3, 21, 0) 169 void * ctx_c =
nullptr;
191 void * ctx_f =
nullptr;
209 PetscInt nfieldsf, nfieldsg;
222 if ( nfieldsf < nfieldsg )
232 for (
int i = 0; i < nfieldsf ; i++)
234 for (
int j = 0; j < nfieldsg ;j++)
235 if ( strcmp( fieldnamesg[j], fieldnamesf[i] ) == 0 )
246 void * ctx_c =
nullptr;
269 Mat * mat ,Vec * vec)
283 void * ctx_c =
nullptr;
289 void * ctx_f =
nullptr;
300 PetscInt nfieldsf, nfieldsg;
311 if ( nfieldsf < nfieldsg)
314 std::vector<std::vector<numeric_index_type>> allrows,allcols;
315 std::vector<numeric_index_type> rows,cols;
322 const int n_subfields = p_ctx_f->
subfields.size();
323 if ( n_subfields >= 1 )
327 rows.insert(rows.end(), allrows[i].begin(), allrows[i].end());
328 cols.insert(cols.end(), allcols[i].begin(), allcols[i].end());
330 std::sort(rows.begin(),rows.end());
331 std::sort(cols.begin(),cols.end());
349 *(vec) = LIBMESH_PETSC_NULLPTR;
369 void * ctx_f =
nullptr;
409 LOG_SCOPE (
"init_and_attach_petscdm()",
"PetscDMWrapper");
415 mesh.allow_renumbering(
false);
416 mesh.allow_remote_element_removal(
false);
417 mesh.partitioner() =
nullptr;
436 unsigned int usr_requested_mg_lvls = 0;
437 usr_requested_mg_lvls =
command_line_next(
"-pc_mg_levels", usr_requested_mg_lvls);
440 if ( usr_requested_mg_lvls != 0 )
443 libmesh_assert_less_equal( usr_requested_mg_lvls,
n_levels );
464 for(
unsigned int level =
n_levels; level >= 1; level--)
471 DM & dm = this->
get_dm(level-1);
472 PetscSection & section = this->
get_section(level-1);
476 LibmeshPetscCall2(system.
comm(), DMShellCreate(system.
comm().
get(), &dm));
479 LibmeshPetscCall2(system.
comm(), DMSetCoordinateDim(dm,
mesh.mesh_dimension()));
483 #if PETSC_VERSION_LESS_THAN(3,10,0) 484 LibmeshPetscCall2(system.
comm(), DMSetDefaultSection(dm, section));
485 #elif PETSC_VERSION_LESS_THAN(3,23,0) 486 LibmeshPetscCall2(system.
comm(), DMSetSection(dm, section));
488 LibmeshPetscCall2(system.
comm(), DMSetLocalSection(dm, section));
495 this->
build_sf(system, star_forest);
496 #if PETSC_VERSION_LESS_THAN(3,12,0) 497 LibmeshPetscCall2(system.
comm(), DMSetDefaultSF(dm, star_forest));
499 LibmeshPetscCall2(system.
comm(), DMSetSectionSF(dm, star_forest));
507 bool supply_restriction =
false;
508 if (supply_restriction)
527 for(
unsigned int i=1; i <=
n_levels; i++ )
563 for (
unsigned int i = 1; i <=
n_levels ; ++i)
565 DM & dm = this->
get_dm(i-1);
567 LibmeshPetscCall2(system.
comm(), DMShellSetGlobalVector( dm,
_ctx_vec[i-1].current_vec->vec() ));
569 LibmeshPetscCall2(system.
comm(), DMShellSetContext( dm, &
_ctx_vec[i-1] ));
582 for(
unsigned int v = 0; v <
n_vars; v++ )
584 std::vector<numeric_index_type> di;
589 LOG_CALL (
"PDM_refine",
"PetscDMWrapper", mesh_refinement.
uniformly_refine(1));
596 for (
unsigned int i = 1 ; i <
n_levels ; ++i )
604 for(
unsigned int v = 0; v <
n_vars; v++ )
606 std::vector<numeric_index_type> di;
621 unsigned int ndofs_old_size = ndofs_old_end - ndofs_old_first;
624 _ctx_vec[i-1].K_interp_ptr->init(ndofs_f, ndofs_c, ndofs_local, ndofs_old_size, 30 , 20);
627 _ctx_vec[i-1].K_interp_ptr->set_destroy_mat_on_exit(
false);
628 _ctx_vec[i-1].K_sub_interp_ptr->set_destroy_mat_on_exit(
false);
637 _ctx_vec[i-1].K_interp_ptr->close();
643 LOG_CALL (
"PDM_refine",
"PetscDMWrapper", mesh_refinement.
uniformly_refine(1));
659 LibmeshPetscCall2(system.
comm(), SNESSetDM(snes, dm));
668 LibmeshPetscCall2(system.
comm(), KSPSetDM(ksp, dm));
673 LOG_SCOPE (
"build_section()",
"PetscDMWrapper");
675 LibmeshPetscCall2(system.
comm(), PetscSectionCreate(system.
comm().
get(),§ion));
678 LibmeshPetscCall2(system.
comm(), PetscSectionSetNumFields(section,system.
n_vars()));
682 LibmeshPetscCall2(system.
comm(), PetscSectionSetFieldName( section, v, system.
variable_name(v).c_str() ));
693 std::unordered_map<dof_id_type,dof_id_type> node_map;
694 std::unordered_map<dof_id_type,dof_id_type> elem_map;
695 std::map<dof_id_type,unsigned int> scalar_map;
710 "ERROR: Must use --node-major-dofs with PetscSection!");
722 LibmeshPetscCall2(system.
comm(), PetscSectionSetUp(section));
730 LOG_SCOPE (
"build_sf()",
"PetscDMWrapper");
734 const std::vector<dof_id_type> & send_list = dof_map.
get_send_list();
737 const PetscInt n_leaves = cast_int<PetscInt>(send_list.size());
740 const PetscInt n_roots = dof_map.
n_local_dofs() + n_leaves;
745 static_assert(
sizeof(PetscInt) ==
sizeof(
dof_id_type),
"PetscInt is not a dof_id_type!");
746 PetscInt * local_dofs =
reinterpret_cast<PetscInt *
>(
const_cast<dof_id_type *
>(send_list.data()));
752 std::set<dof_id_type> send_set(send_list.begin(), send_list.end());
753 libmesh_assert_equal_to(send_list.size(), send_set.size());
762 std::vector<PetscSFNode> sf_nodes(send_list.size());
772 PetscInt index = incoming_dof - dof_map.
first_dof(rank);
774 sf_nodes[i].rank = rank;
775 sf_nodes[i].index = index;
778 PetscSFNode * remote_dofs = sf_nodes.data();
780 LibmeshPetscCall2(system.
comm(), PetscSFCreate(system.
comm().
get(), &star_forest));
785 LibmeshPetscCall2(system.
comm(), PetscSFSetGraph(star_forest,
795 PetscSection & section,
796 std::unordered_map<dof_id_type,dof_id_type> & node_map,
797 std::unordered_map<dof_id_type,dof_id_type> & elem_map,
798 std::map<dof_id_type,unsigned int> & scalar_map)
836 node_map.reserve(
mesh.n_local_nodes());
837 elem_map.reserve(
mesh.n_active_local_elem());
843 for (
const auto & elem :
mesh.active_local_element_ptr_range())
845 for (
const Node & node : elem->node_ref_range())
853 std::vector<dof_id_type> node_dof_indices;
855 if( !node_dof_indices.empty() && dof_map.
local_index(node_dof_indices[0]) )
860 for(
auto dof : node_dof_indices )
866 if( node_map.count(node_id) == 0 )
868 node_map.emplace(node_id, pend);
879 if( elem->n_dofs(system.
number()) > 0 )
882 elem_map.emplace(elem_id, pend);
897 scalar_map.emplace(pend, v);
905 LibmeshPetscCall2(system.
comm(), PetscSectionSetChart(section, pstart, pend));
909 PetscSection & section,
910 const std::unordered_map<dof_id_type,dof_id_type> & node_map,
911 const std::unordered_map<dof_id_type,dof_id_type> & elem_map,
912 const std::map<dof_id_type,unsigned int> & scalar_map)
922 for (
const auto [global_node_id, local_node_id] : node_map )
926 const Node & node =
mesh.node_ref(global_node_id);
934 for (
const auto [global_elem_id, local_elem_id] : elem_map )
938 const Elem & elem =
mesh.elem_ref(global_elem_id);
947 for (
const auto [local_id, scalar_var] : scalar_map )
952 LibmeshPetscCall2(system.
comm(), PetscSectionSetFieldDof( section, local_id, scalar_var, n_dofs ));
956 LibmeshPetscCall2(system.
comm(), PetscSectionSetDof( section, local_id, n_dofs ));
965 PetscSection & section)
967 unsigned int total_n_dofs_at_dofobject = 0;
972 unsigned int n_dofs_at_dofobject = dof_object.
n_dofs(system.
number(), v);
974 if( n_dofs_at_dofobject > 0 )
976 LibmeshPetscCall2(system.
comm(), PetscSectionSetFieldDof( section,
979 n_dofs_at_dofobject ));
981 total_n_dofs_at_dofobject += n_dofs_at_dofobject;
985 libmesh_assert_equal_to(total_n_dofs_at_dofobject, dof_object.
n_dofs(system.
number()));
987 LibmeshPetscCall2(system.
comm(), PetscSectionSetDof( section, local_id, total_n_dofs_at_dofobject ));
993 PetscInt n_local_dofs = 0;
996 PetscInt pstart, pend;
997 PetscErrorCode ierr = PetscSectionGetChart(section, &pstart, &pend);
1002 for( PetscInt p = pstart; p < pend; p++ )
1005 ierr = PetscSectionGetDof(section,p,&n_dofs);
1008 n_local_dofs += n_dofs;
1011 static_assert(
sizeof(PetscInt) ==
sizeof(
dof_id_type),
"PetscInt is not a dof_id_type!");
1012 return n_local_dofs;
1027 for(
unsigned int i = 0; i <
n_levels; i++ )
1030 _pmtx_vec[i] = std::make_unique<PetscMatrix<Number>>(comm);
1031 _subpmtx_vec[i] = std::make_unique<PetscMatrix<Number>>(comm);
1032 _vec_vec[i] = std::make_unique<PetscVector<Number>>(comm);
1038 #endif // #if LIBMESH_ENABLE_AMR && LIBMESH_HAVE_METAPHYSICL 1039 #endif // PETSC_VERSION 1040 #endif // LIBMESH_HAVE_PETSC void init_dm_data(unsigned int n_levels, const Parallel::Communicator &comm)
Init all the n_mesh_level dependent data structures.
PetscErrorCode PetscInt const PetscInt IS DM * subdm
FEFamily family
The type of finite element.
dof_id_type check_section_n_dofs(PetscSection §ion)
Helper function to sanity check PetscSection construction The PetscSection contains local dof informa...
std::vector< std::unique_ptr< PetscVector< Number > > > _vec_vec
Vector of solution vectors for all grid levels.
T command_line_next(std::string name, T default_value)
Use GetPot's search()/next() functions to get following arguments from the command line...
virtual void create_submatrix(SparseMatrix< T > &submatrix, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols) const
This function creates a matrix called "submatrix" which is defined by the row and column indices give...
std::vector< PetscDMContext > _ctx_vec
Vector of internal PetscDM context structs for all grid levels Pointers to these C++ objects are pass...
A Node is like a Point, but with more information.
dof_id_type n_SCALAR_dofs() const
const Variable & variable(unsigned int var) const
Return a constant reference to Variable var.
void set_point_range_in_section(const System &system, PetscSection §ion, std::unordered_map< dof_id_type, dof_id_type > &node_map, std::unordered_map< dof_id_type, dof_id_type > &elem_map, std::map< dof_id_type, unsigned int > &scalar_map)
Helper function for build_section.
processor_id_type dof_owner(const dof_id_type dof) const
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
std::size_t distribute_dofs(MeshBase &)
Distribute dofs on the current mesh.
DM & get_dm(unsigned int level)
Get reference to DM for the given mesh level.
void local_variable_indices(T &idx, const MeshBase &mesh, unsigned int var_num) const
If T == dof_id_type, counts, if T == std::vector<dof_id_type>, fills an array of, those dof indices w...
dof_id_type end_old_dof(const processor_id_type proc) const
std::vector< WrappedPetsc< DM > > _dms
Vector of DMs for all grid levels.
void projection_matrix(SparseMatrix< Number > &proj_mat) const
This method creates a projection matrix which corresponds to the operation of project_vector between ...
unsigned int init_petscdm(System &system)
Initialize the PETSc DM and return the number of geometric multigrid levels.
This is the base class from which all geometric element types are derived.
std::vector< PetscInt > subfields
Stores subfield ids for use in subprojection matrixes on coarser DMs.
const Parallel::Communicator & comm() const
void add_dofs_to_section(const System &system, PetscSection §ion, const std::unordered_map< dof_id_type, dof_id_type > &node_map, const std::unordered_map< dof_id_type, dof_id_type > &elem_map, const std::map< dof_id_type, unsigned int > &scalar_map)
Helper function for build_section.
PetscMatrixBase< libMesh::Number > * K_interp_ptr
OrderWrapper order
The approximation order of the element.
dof_id_type n_dofs(const unsigned int vn) const
The libMesh namespace provides an interface to certain functionality in the library.
dof_id_type n_local_dofs() const
dof_id_type n_local_dofs(const unsigned int vn) const
PetscErrorCode libmesh_petsc_DMCreateInterpolation(DM dmc, DM dmf, Mat *mat, Vec *vec)
Function to give PETSc that sets the Interpolation Matrix between two DMs.
const MeshBase & get_mesh() const
PetscErrorCode libmesh_petsc_DMRefine(DM dmc, MPI_Comm, DM *dmf)
Help PETSc identify the finer DM given a dmc.
PetscErrorCode PetscInt numFields
PetscErrorCode libmesh_petsc_DMCreateRestriction(DM dmc, DM dmf, Mat *mat)
Function to give PETSc that sets the Restriction Matrix between two DMs.
PetscMatrixBase< libMesh::Number > * K_restrict_ptr
uint8_t processor_id_type
This is the MeshBase class.
unsigned int n_dofs(const unsigned int s, const unsigned int var=libMesh::invalid_uint) const
This class handles the numbering of degrees of freedom on a mesh.
processor_id_type n_processors() const
unsigned int number() const
Struct to house data regarding where in the mesh hierarchy we are located.
std::vector< std::unique_ptr< PetscMatrixBase< Number > > > _pmtx_vec
Vector of projection matrixes for all grid levels.
Implements (adaptive) mesh refinement algorithms for a MeshBase.
std::vector< WrappedPetsc< PetscSection > > _sections
Vector of PETScSections for all grid levels.
Manages consistently variables, degrees of freedom, and coefficient vectors.
void uniformly_coarsen(unsigned int n=1)
Attempts to uniformly coarsen the mesh n times.
PetscErrorCode PetscInt const PetscInt IS * is
const std::string & variable_name(const unsigned int i) const
std::vector< std::vector< numeric_index_type > > dof_vec
Stores local dofs for each var for use in subprojection matrixes.
PetscSection & get_section(unsigned int level)
Get reference to PetscSection for the given mesh level.
PetscErrorCode libmesh_petsc_DMCoarsen(DM dmf, MPI_Comm, DM *dmc)
Help PETSc identify the coarser DM dmc given the fine DM dmf.
PetscErrorCode libmesh_petsc_DMCreateSubDM(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm) PetscErrorCode libmesh_petsc_DMCreateSubDM(DM dm
Help PETSc create a subDM given a global dm when using fieldsplit.
std::vector< unsigned int > _mesh_dof_loc_sizes
Stores n_local_dofs for each grid level, to be used for projection vector sizing. ...
virtual void reinit_constraints()
Reinitializes the constraints for this system.
std::vector< PetscSF > _star_forests
Vector of star forests for all grid levels.
int get_order() const
Explicitly request the order as an int.
PetscSF & get_star_forest(unsigned int level)
Get reference to PetscSF for the given mesh level.
PetscErrorCode PetscInt const PetscInt fields[]
LibmeshPetscCallQ(DMShellGetContext(dm, &ctx))
void add_dofs_helper(const System &system, const DofObject &dof_object, dof_id_type local_id, PetscSection §ion)
PetscMatrixBase< libMesh::Number > * K_sub_interp_ptr
void init_and_attach_petscdm(System &system, SNES snes)
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
std::vector< unsigned int > _mesh_dof_sizes
Stores n_dofs for each grid level, to be used for projection matrix sizing.
The DofObject defines an abstract base class for objects that have degrees of freedom associated with...
void prepare_send_list()
Takes the _send_list vector (which may have duplicate entries) and sorts it.
bool on_command_line(std::string arg)
dof_id_type first_old_dof(const processor_id_type proc) const
dof_id_type first_dof(const processor_id_type proc) const
void build_section(const System &system, PetscSection §ion)
Takes System, empty PetscSection and populates the PetscSection.
void build_sf(const System &system, PetscSF &star_forest)
Takes System, empty PetscSF and populates the PetscSF.
unsigned int n_vars() const
PetscFunctionReturn(LIBMESH_PETSC_SUCCESS)
processor_id_type processor_id() const
std::vector< std::unique_ptr< PetscMatrixBase< Number > > > _subpmtx_vec
Vector of sub projection matrixes for all grid levels for fieldsplit.
void clear()
Destroys and clears all build DM-related data.
const DofMap & get_dof_map() const
const std::vector< dof_id_type > & get_send_list() const
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
bool local_index(dof_id_type dof_index) const
const FEType & type() const
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.