libMesh::PetscDMWrapper Class Reference

This class defines a wrapper around the PETSc DM infrastructure. More...

#include <petsc_dm_wrapper.h>

Public Member Functions
	PetscDMWrapper ()=default
	~PetscDMWrapper ()
void	clear ()
	Destroys and clears all build DM-related data. More...
void	init_and_attach_petscdm (System &system, SNES snes)
void	init_and_attach_petscdm (System &system, KSP ksp)

Private Member Functions
unsigned int	init_petscdm (System &system)
	Initialize the PETSc DM and return the number of geometric multigrid levels. More...
void	init_dm_data (unsigned int n_levels, const Parallel::Communicator &comm)
	Init all the n_mesh_level dependent data structures. More...
DM &	get_dm (unsigned int level)
	Get reference to DM for the given mesh level. More...
PetscSection &	get_section (unsigned int level)
	Get reference to PetscSection for the given mesh level. More...
PetscSF &	get_star_forest (unsigned int level)
	Get reference to PetscSF for the given mesh level. More...
void	build_section (const System &system, PetscSection &section)
	Takes System, empty PetscSection and populates the PetscSection. More...
void	build_sf (const System &system, PetscSF &star_forest)
	Takes System, empty PetscSF and populates the PetscSF. More...
void	set_point_range_in_section (const System &system, PetscSection &section, 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. More...
void	add_dofs_to_section (const System &system, PetscSection &section, 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. More...
dof_id_type	check_section_n_dofs (PetscSection &section)
	Helper function to sanity check PetscSection construction The PetscSection contains local dof information. More...
void	add_dofs_helper (const System &system, const DofObject &dof_object, dof_id_type local_id, PetscSection &section)

Private Attributes
std::vector< WrappedPetsc< DM > >	_dms
	Vector of DMs for all grid levels. More...
std::vector< WrappedPetsc< PetscSection > >	_sections
	Vector of PETScSections for all grid levels. More...
std::vector< PetscSF >	_star_forests
	Vector of star forests for all grid levels. More...
std::vector< std::unique_ptr< PetscMatrixBase< Number > > >	_pmtx_vec
	Vector of projection matrixes for all grid levels. More...
std::vector< std::unique_ptr< PetscMatrixBase< Number > > >	_subpmtx_vec
	Vector of sub projection matrixes for all grid levels for fieldsplit. More...
std::vector< PetscDMContext >	_ctx_vec
	Vector of internal PetscDM context structs for all grid levels Pointers to these C++ objects are passed to DMShellSetContext(), they are cleaned up automatically by their destructors. More...
std::vector< std::unique_ptr< PetscVector< Number > > >	_vec_vec
	Vector of solution vectors for all grid levels. More...
std::vector< unsigned int >	_mesh_dof_sizes
	Stores n_dofs for each grid level, to be used for projection matrix sizing. More...
std::vector< unsigned int >	_mesh_dof_loc_sizes
	Stores n_local_dofs for each grid level, to be used for projection vector sizing. More...

Detailed Description

This class defines a wrapper around the PETSc DM infrastructure.

By coordinating DM data structures with libMesh, we can use libMesh mesh hierarchies for geometric multigrid. Additionally, by setting the DM data, we can additionally (with or without multigrid) define recursive fieldsplits of our variables.

Author

Paul T. Bauman, Boris Boutkov

Date

2018

Definition at line 97 of file petsc_dm_wrapper.h.

Constructor & Destructor Documentation

PetscDMWrapper()

libMesh::PetscDMWrapper::PetscDMWrapper ( )

default

~PetscDMWrapper()

libMesh::PetscDMWrapper::~PetscDMWrapper ( )

default

Member Function Documentation

add_dofs_helper()

void libMesh::PetscDMWrapper::add_dofs_helper ( const System &  system, const DofObject &  dof_object, dof_id_type  local_id, PetscSection &  section  )

private

Definition at line 962 of file petsc_dm_wrapper.C.

References libMesh::ParallelObject::comm(), libMesh::make_range(), libMesh::DofObject::n_dofs(), libMesh::System::n_vars(), and libMesh::System::number().

Referenced by add_dofs_to_section().

  {
    unsigned int total_n_dofs_at_dofobject = 0;

    // We are assuming variables are also numbered 0 to n_vars()-1
    for (auto v : make_range(system.n_vars()))
      {
        unsigned int n_dofs_at_dofobject = dof_object.n_dofs(system.number(), v);

        if( n_dofs_at_dofobject > 0 )
          {
            LibmeshPetscCall2(system.comm(), PetscSectionSetFieldDof( section,
                                                                      local_id,
                                                                      v,
                                                                      n_dofs_at_dofobject ));

            total_n_dofs_at_dofobject += n_dofs_at_dofobject;
          }
      }

    libmesh_assert_equal_to(total_n_dofs_at_dofobject, dof_object.n_dofs(system.number()));

    LibmeshPetscCall2(system.comm(), PetscSectionSetDof( section, local_id, total_n_dofs_at_dofobject ));
  }

add_dofs_to_section()

void libMesh::PetscDMWrapper::add_dofs_to_section ( const System &  system, PetscSection &  section, 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  )

private

Helper function for build_section.

This function will set the DoF info for each "point" in the PetscSection.

Definition at line 908 of file petsc_dm_wrapper.C.

References add_dofs_helper(), libMesh::ParallelObject::comm(), libMesh::System::get_mesh(), libMesh::OrderWrapper::get_order(), libMesh::libmesh_assert(), mesh, libMesh::ParallelObject::n_processors(), libMesh::FEType::order, libMesh::ParallelObject::processor_id(), libMesh::Variable::type(), and libMesh::System::variable().

Referenced by build_section().

  {
    const MeshBase & mesh = system.get_mesh();

    // Now we go through and add dof information for each "point".
    //
    // In libMesh, for most finite elements, we just associate those DoFs with the
    // geometric nodes. So can we loop over the nodes we cached in the node_map and
    // the DoFs for each field for that node. We need to give PETSc the local id
    // we built up in the node map.
    for (const auto [global_node_id, local_node_id] : node_map )
      {
        libmesh_assert( mesh.query_node_ptr(global_node_id) );

        const Node & node = mesh.node_ref(global_node_id);

        this->add_dofs_helper(system,node,local_node_id,section);
      }

    // Some finite element types associate dofs with the element. So now we go through
    // any of those with the Elem as the point we add to the PetscSection with accompanying
    // dofs
    for (const auto [global_elem_id, local_elem_id] : elem_map )
      {
        libmesh_assert( mesh.query_elem_ptr(global_elem_id) );

        const Elem & elem = mesh.elem_ref(global_elem_id);

        this->add_dofs_helper(system,elem,local_elem_id,section);
      }

    // Now add any SCALAR dofs to the PetscSection
    // SCALAR dofs live on the "last" processor, so only work there if there are any
    if (system.processor_id() == (system.n_processors()-1))
      {
        for (const auto [local_id, scalar_var] : scalar_map )
          {
            // The number of SCALAR dofs comes from the variable order
            const int n_dofs = system.variable(scalar_var).type().order.get_order();

            LibmeshPetscCall2(system.comm(), PetscSectionSetFieldDof( section, local_id, scalar_var, n_dofs ));

            // In the SCALAR case, there are no other variables associate with the "point"
            // the total number of dofs on the point is the same as that for the field
            LibmeshPetscCall2(system.comm(), PetscSectionSetDof( section, local_id, n_dofs ));
          }
      }

  }

build_section()

void libMesh::PetscDMWrapper::build_section ( const System &  system, PetscSection &  section  )

private

Takes System, empty PetscSection and populates the PetscSection.

Take the System in its current state and an empty PetscSection and then populate the PetscSection. The PetscSection is comprised of global "point" numbers, where a "point" in PetscDM parlance is a geometric entity: node, edge, face, or element. Then, we also add the DoF numbering for each variable for each of the "points". The PetscSection, together the with PetscSF will allow for recursive fieldsplits from the command line using PETSc.

Definition at line 671 of file petsc_dm_wrapper.C.

References add_dofs_to_section(), check_section_n_dofs(), libMesh::ParallelObject::comm(), TIMPI::Communicator::get(), libMesh::make_range(), libMesh::System::n_local_dofs(), libMesh::System::n_vars(), libMesh::on_command_line(), set_point_range_in_section(), and libMesh::System::variable_name().

Referenced by init_petscdm().

  {
    LOG_SCOPE ("build_section()", "PetscDMWrapper");

    LibmeshPetscCall2(system.comm(), PetscSectionCreate(system.comm().get(),&section));

    // Tell the PetscSection about all of our System variables
    LibmeshPetscCall2(system.comm(), PetscSectionSetNumFields(section,system.n_vars()));

    // Set the actual names of all the field variables
    for (auto v : make_range(system.n_vars()))
      LibmeshPetscCall2(system.comm(), PetscSectionSetFieldName( section, v, system.variable_name(v).c_str() ));

    // For building the section, we need to create local-to-global map
    // of local "point" ids to the libMesh global id of that point.
    // A "point" in PETSc nomenclature is a geometric object that can have
    // dofs associated with it, e.g. Node, Edge, Face, Elem.
    // The numbering PETSc expects is continuous for the local numbering.
    // Since we're only using this interface for solvers, then we can just
    // assign whatever local id to any of the global ids. But it is local
    // so we don't need to worry about processor numbering for the local
    // point ids.
    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;

    // First we tell the PetscSection about all of our points that have
    // dofs associated with them.
    this->set_point_range_in_section(system, section, node_map, elem_map, scalar_map);

    // Now we can build up the dofs per "point" in the PetscSection
    this->add_dofs_to_section(system, section, node_map, elem_map, scalar_map);

    // Final setup of PetscSection
    // Until Matt Knepley finishes implementing the commented out function
    // below, the PetscSection will be assuming node-major ordering
    // so let's throw an error if the user tries to use this without
    // node-major order
    libmesh_error_msg_if(!libMesh::on_command_line("--node-major-dofs"),
                         "ERROR: Must use --node-major-dofs with PetscSection!");

    //else if (!system.identify_variable_groups())
    //  ierr = PetscSectionSetUseFieldOffsets(section,PETSC_TRUE);LIBMESH_CHKERR(ierr);
    //else
    //  {
    //    std::string msg = "ERROR: Only node-major or var-major ordering supported for PetscSection!\n";
    //    msg += "       var-group-major ordering not supported!\n";
    //    msg += "       Must use --node-major-dofs or set System::identify_variable_groups() = false!\n";
    //    libmesh_error_msg(msg);
    //  }

    LibmeshPetscCall2(system.comm(), PetscSectionSetUp(section));

    // Sanity checking at least that local_n_dofs match
    libmesh_assert_equal_to(system.n_local_dofs(),this->check_section_n_dofs(section));
  }

build_sf()

void libMesh::PetscDMWrapper::build_sf ( const System &  system, PetscSF &  star_forest  )

private

Takes System, empty PetscSF and populates the PetscSF.

The PetscSF (star forest) is a cousin of PetscSection. PetscSection has the DoF info, and PetscSF gives the parallel distribution of the DoF info. So PetscSF should only be necessary when we have more than one MPI rank. Essentially, we are copying the DofMap.send_list(): we are specifying the local dofs, what rank communicates that dof info (for off-processor dofs that are communicated) and the dofs local index on that rank.

https://jedbrown.org/files/StarForest.pdf

Definition at line 728 of file petsc_dm_wrapper.C.

References libMesh::ParallelObject::comm(), libMesh::DofMap::dof_owner(), libMesh::DofMapBase::first_dof(), TIMPI::Communicator::get(), libMesh::System::get_dof_map(), libMesh::DofMap::get_send_list(), libMesh::index_range(), and libMesh::DofMap::n_local_dofs().

Referenced by init_petscdm().

  {
    LOG_SCOPE ("build_sf()", "PetscDMWrapper");

    const DofMap & dof_map = system.get_dof_map();

    const std::vector<dof_id_type> & send_list = dof_map.get_send_list();

    // Number of ghost dofs that send information to this processor
    const PetscInt n_leaves = cast_int<PetscInt>(send_list.size());

    // Number of local dofs, including ghosts dofs
    const PetscInt n_roots = dof_map.n_local_dofs() + n_leaves;

    // This is the vector of dof indices coming from other processors
    // We need to give this to the PetscSF
    // We'll be extra paranoid about this ugly double cast
    static_assert(sizeof(PetscInt) == sizeof(dof_id_type),"PetscInt is not a dof_id_type!");
    PetscInt * local_dofs = reinterpret_cast<PetscInt *>(const_cast<dof_id_type *>(send_list.data()));

#ifdef DEBUG
    // PETSc 3.18 and above don't want duplicate entries here ... and
    // frankly we shouldn't have duplicates in the first place!
    {
      std::set<dof_id_type> send_set(send_list.begin(), send_list.end());
      libmesh_assert_equal_to(send_list.size(), send_set.size());
    }
#endif

    // This is the vector of PetscSFNode's for the local_dofs.
    // For each entry in local_dof, we have to supply the rank from which
    // that dof stems and its local index on that rank.
    // PETSc documentation here:
    // http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PetscSF/PetscSFNode.html
    std::vector<PetscSFNode> sf_nodes(send_list.size());

    for (auto i : index_range(send_list))
      {
        dof_id_type incoming_dof = send_list[i];

        const processor_id_type rank = dof_map.dof_owner(incoming_dof);

        // Dofs are sorted and continuous on the processor so local index
        // is counted up from the first dof on the processor.
        PetscInt index = incoming_dof - dof_map.first_dof(rank);

        sf_nodes[i].rank  = rank; /* Rank of owner */
        sf_nodes[i].index = index;/* Index of dof on rank */
      }

    PetscSFNode * remote_dofs = sf_nodes.data();

    LibmeshPetscCall2(system.comm(), PetscSFCreate(system.comm().get(), &star_forest));

    // TODO: We should create pointers to arrays so we don't have to copy
    //       and then can use PETSC_OWN_POINTER where PETSc will take ownership
    //       and delete the memory for us. But then we'd have to use PetscMalloc.
    LibmeshPetscCall2(system.comm(), PetscSFSetGraph(star_forest,
                                                     n_roots,
                                                     n_leaves,
                                                     local_dofs,
                                                     PETSC_COPY_VALUES,
                                                     remote_dofs,
                                                     PETSC_COPY_VALUES));
  }

check_section_n_dofs()

dof_id_type libMesh::PetscDMWrapper::check_section_n_dofs ( PetscSection &  section )

private

Helper function to sanity check PetscSection construction The PetscSection contains local dof information.

This helper function just facilitates sanity checking that in fact it only has n_local_dofs.

Definition at line 991 of file petsc_dm_wrapper.C.

Referenced by build_section().

  {
    PetscInt n_local_dofs = 0;

    // Grap the starting and ending points from the section
    PetscInt pstart, pend;
    PetscErrorCode ierr = PetscSectionGetChart(section, &pstart, &pend);
    if (ierr)
      libmesh_error();

    // Count up the n_dofs for each point from the section
    for( PetscInt p = pstart; p < pend; p++ )
      {
        PetscInt n_dofs;
        ierr = PetscSectionGetDof(section,p,&n_dofs);
        if (ierr)
          libmesh_error();
        n_local_dofs += n_dofs;
      }

    static_assert(sizeof(PetscInt) == sizeof(dof_id_type),"PetscInt is not a dof_id_type!");
    return n_local_dofs;
  }

clear()

void libMesh::PetscDMWrapper::clear

(

)

Destroys and clears all build DM-related data.

Definition at line 386 of file petsc_dm_wrapper.C.

References _ctx_vec, _dms, _mesh_dof_loc_sizes, _mesh_dof_sizes, _pmtx_vec, _sections, _star_forests, _subpmtx_vec, and _vec_vec.

Referenced by libMesh::PetscDiffSolver::clear().

  {
    // Calls custom deleters
    _dms.clear();
    _sections.clear();

    // We don't manage the lifetime of the PetscSF objects
    _star_forests.clear();

    // The relevant C++ destructors are called for these objects
    // automatically.
    _pmtx_vec.clear();
    _subpmtx_vec.clear();
    _vec_vec.clear();

    // These members are trivially clear()able.
    _ctx_vec.clear();
    _mesh_dof_sizes.clear();
    _mesh_dof_loc_sizes.clear();
  }

get_dm()

DM& libMesh::PetscDMWrapper::get_dm ( unsigned int  level )

inlineprivate

Get reference to DM for the given mesh level.

init_dm_data() should be called before this function.

Definition at line 184 of file petsc_dm_wrapper.h.

References _dms.

Referenced by init_and_attach_petscdm(), and init_petscdm().

      { libmesh_assert_less(level, _dms.size());
        return *_dms[level]; }

get_section()

PetscSection& libMesh::PetscDMWrapper::get_section ( unsigned int  level )

inlineprivate

Get reference to PetscSection for the given mesh level.

init_dm_data() should be called before this function.

Definition at line 192 of file petsc_dm_wrapper.h.

References _sections.

Referenced by init_petscdm().

      { libmesh_assert_less(level, _sections.size());
        return *_sections[level]; }

get_star_forest()

PetscSF& libMesh::PetscDMWrapper::get_star_forest ( unsigned int  level )

inlineprivate

Get reference to PetscSF for the given mesh level.

init_dm_data() should be called before this function.

Definition at line 200 of file petsc_dm_wrapper.h.

References _star_forests.

Referenced by init_petscdm().

      { libmesh_assert_less(level, _star_forests.size());
        return _star_forests[level]; }

init_and_attach_petscdm() [1/2]

void libMesh::PetscDMWrapper::init_and_attach_petscdm	(	System &	system,
		SNES	snes
	)

Definition at line 653 of file petsc_dm_wrapper.C.

References libMesh::ParallelObject::comm(), get_dm(), init_petscdm(), and libMesh::MeshTools::n_levels().

Referenced by libMesh::PetscDiffSolver::setup_petsc_data().

  {
    const auto n_levels = this->init_petscdm(system);

    // Lastly, give SNES the finest level DM
    DM & dm = this->get_dm(n_levels-1);
    LibmeshPetscCall2(system.comm(), SNESSetDM(snes, dm));
  }

init_and_attach_petscdm() [2/2]

void libMesh::PetscDMWrapper::init_and_attach_petscdm	(	System &	system,
		KSP	ksp
	)

Definition at line 662 of file petsc_dm_wrapper.C.

References libMesh::ParallelObject::comm(), get_dm(), init_petscdm(), and libMesh::MeshTools::n_levels().

  {
    const auto n_levels = this->init_petscdm(system);

    // Lastly, give KSP the finest level DM
    DM & dm = this->get_dm(n_levels-1);
    LibmeshPetscCall2(system.comm(), KSPSetDM(ksp, dm));
  }

init_dm_data()

void libMesh::PetscDMWrapper::init_dm_data ( unsigned int  n_levels, const Parallel::Communicator &  comm  )

private

Init all the n_mesh_level dependent data structures.

Definition at line 1015 of file petsc_dm_wrapper.C.

References _ctx_vec, _dms, _mesh_dof_loc_sizes, _mesh_dof_sizes, _pmtx_vec, _sections, _star_forests, _subpmtx_vec, _vec_vec, and libMesh::MeshTools::n_levels().

Referenced by init_petscdm().

  {
    _dms.resize(n_levels);
    _sections.resize(n_levels);
    _star_forests.resize(n_levels);
    _ctx_vec.resize(n_levels);
    _pmtx_vec.resize(n_levels);
    _subpmtx_vec.resize(n_levels);
    _vec_vec.resize(n_levels);
    _mesh_dof_sizes.resize(n_levels);
    _mesh_dof_loc_sizes.resize(n_levels);

    for( unsigned int i = 0; i < n_levels; i++ )
      {
        // Call C++ object constructors
        _pmtx_vec[i] = std::make_unique<PetscMatrix<Number>>(comm);
        _subpmtx_vec[i] = std::make_unique<PetscMatrix<Number>>(comm);
        _vec_vec[i] = std::make_unique<PetscVector<Number>>(comm);
      }
  }

init_petscdm()

unsigned int libMesh::PetscDMWrapper::init_petscdm ( System &  system )

private

Initialize the PETSc DM and return the number of geometric multigrid levels.

Definition at line 407 of file petsc_dm_wrapper.C.

References _ctx_vec, _dms, _mesh_dof_loc_sizes, _mesh_dof_sizes, _pmtx_vec, _subpmtx_vec, _vec_vec, build_section(), build_sf(), libMesh::ParallelObject::comm(), libMesh::command_line_next(), libMesh::DofMap::distribute_dofs(), libMesh::DofMapBase::end_old_dof(), libMesh::DofMapBase::first_old_dof(), TIMPI::Communicator::get(), get_dm(), libMesh::System::get_dof_map(), libMesh::System::get_mesh(), get_section(), get_star_forest(), init_dm_data(), libMesh::libmesh_assert(), libMesh::libmesh_petsc_DMCoarsen(), libMesh::libmesh_petsc_DMCreateInterpolation(), libMesh::libmesh_petsc_DMCreateRestriction(), libMesh::libmesh_petsc_DMCreateSubDM(), libMesh::libmesh_petsc_DMRefine(), libMesh::DofMap::local_variable_indices(), mesh, libMesh::DofMap::n_dofs(), libMesh::MeshTools::n_levels(), libMesh::DofMap::n_local_dofs(), libMesh::ParallelObject::n_processors(), n_vars, libMesh::System::n_vars(), libMesh::DofMap::prepare_send_list(), libMesh::System::projection_matrix(), libMesh::System::reinit_constraints(), libMesh::MeshRefinement::uniformly_coarsen(), and libMesh::MeshRefinement::uniformly_refine().

Referenced by init_and_attach_petscdm().

  {
    LOG_SCOPE ("init_and_attach_petscdm()", "PetscDMWrapper");

    MeshBase & mesh = system.get_mesh();   // Convenience
    MeshRefinement mesh_refinement(mesh); // Used for swapping between grids

    // There's no need for these code paths while traversing the hierarchy
    mesh.allow_renumbering(false);
    mesh.allow_remote_element_removal(false);
    mesh.partitioner() = nullptr;

    // First walk over the active local elements and see how many maximum MG levels we can construct
/*
    unsigned int n_levels = 0;
    for ( auto & elem : mesh.active_local_element_ptr_range() )
      {
        if ( elem->level() > n_levels )
          n_levels = elem->level();
      }
    // On coarse grids some processors may have no active local elements,
    // these processors shouldn't make projections
    if (n_levels >= 1)
      n_levels += 1;
*/

    unsigned int n_levels = MeshTools::n_levels(mesh);

    // How many MG levels did the user request?
    unsigned int usr_requested_mg_lvls = 0;
    usr_requested_mg_lvls = command_line_next("-pc_mg_levels", usr_requested_mg_lvls);

    // Only construct however many levels were requested if something was actually requested
    if ( usr_requested_mg_lvls != 0 )
      {
        // Dont request more than avail num levels on mesh, require at least 2 levels
        libmesh_assert_less_equal( usr_requested_mg_lvls, n_levels );
        libmesh_assert( usr_requested_mg_lvls > 1 );

        n_levels = usr_requested_mg_lvls;
      }
    else
      {
        // if -pc_mg_levels is not specified we just construct fieldsplit related
        // structures on the finest mesh.
        n_levels = 1;
      }


    // Init data structures: data[0] ~ coarse grid, data[n_levels-1] ~ fine grid
    this->init_dm_data(n_levels, system.comm());

    // Step 1.  contract : all active elements have no children
    mesh.contract();

    // Start on finest grid. Construct DM datas and stash some info for
    // later projection_matrix and vec sizing
    for(unsigned int level = n_levels; level >= 1; level--)
      {
        // Save the n_fine_dofs before coarsening for later projection matrix sizing
        _mesh_dof_sizes[level-1] = system.get_dof_map().n_dofs();
        _mesh_dof_loc_sizes[level-1] = system.get_dof_map().n_local_dofs();

        // Get refs to things we will fill
        DM & dm = this->get_dm(level-1);
        PetscSection & section = this->get_section(level-1);
        PetscSF & star_forest = this->get_star_forest(level-1);

        // The shell will contain other DM info
        LibmeshPetscCall2(system.comm(), DMShellCreate(system.comm().get(), &dm));

        // Set the DM embedding dimension to help PetscDS (Discrete System)
        LibmeshPetscCall2(system.comm(), DMSetCoordinateDim(dm, mesh.mesh_dimension()));

        // Build the PetscSection and attach it to the DM
        this->build_section(system, section);
#if PETSC_VERSION_LESS_THAN(3,10,0)
        LibmeshPetscCall2(system.comm(), DMSetDefaultSection(dm, section));
#elif PETSC_VERSION_LESS_THAN(3,23,0)
        LibmeshPetscCall2(system.comm(), DMSetSection(dm, section));
#else
        LibmeshPetscCall2(system.comm(), DMSetLocalSection(dm, section));
#endif

        // We only need to build the star forest if we're in a parallel environment
        if (system.n_processors() > 1)
          {
            // Build the PetscSF and attach it to the DM
            this->build_sf(system, star_forest);
#if PETSC_VERSION_LESS_THAN(3,12,0)
            LibmeshPetscCall2(system.comm(), DMSetDefaultSF(dm, star_forest));
#else
            LibmeshPetscCall2(system.comm(), DMSetSectionSF(dm, star_forest));
#endif
          }

        // Set PETSC's Restriction, Interpolation, Coarsen and Refine functions for the current DM
        LibmeshPetscCall2(system.comm(), DMShellSetCreateInterpolation ( dm, libmesh_petsc_DMCreateInterpolation ));

        // Not implemented. For now we rely on galerkin style restrictions
        bool supply_restriction = false;
        if (supply_restriction)
          LibmeshPetscCall2(system.comm(), DMShellSetCreateRestriction ( dm, libmesh_petsc_DMCreateRestriction  ));

        LibmeshPetscCall2(system.comm(), DMShellSetCoarsen ( dm, libmesh_petsc_DMCoarsen ));

        LibmeshPetscCall2(system.comm(), DMShellSetRefine ( dm, libmesh_petsc_DMRefine ));

        LibmeshPetscCall2(system.comm(), DMShellSetCreateSubDM(dm, libmesh_petsc_DMCreateSubDM));

        // Uniformly coarsen if not the coarsest grid and distribute dof info.
        if ( level != 1 )
          {
            LOG_CALL("PDM_coarsen", "PetscDMWrapper", mesh_refinement.uniformly_coarsen(1));
            LOG_CALL("PDM_dist_dof", "PetscDMWrapper", system.get_dof_map().distribute_dofs(mesh));
            LOG_CALL("PDM_prep_send", "PetscDMWrapper", system.get_dof_map().prepare_send_list());
          }
      } // End PETSc data structure creation

    // Now fill the corresponding internal PetscDMContext for each created DM
    for( unsigned int i=1; i <= n_levels; i++ )
      {
        // Set context dimension
        _ctx_vec[i-1].mesh_dim = mesh.mesh_dimension();

        // Create and attach a sized vector to the current ctx
        _vec_vec[i-1]->init( _mesh_dof_sizes[i-1] );
        _ctx_vec[i-1].current_vec = _vec_vec[i-1].get();

        // Set a global DM to be used as reference when using fieldsplit
        _ctx_vec[i-1].global_dm = &(this->get_dm(n_levels-1));

        if (n_levels > 1 )
          {
            // Set pointers to surrounding dm levels to help PETSc refine/coarsen
            if ( i == 1 ) // were at the coarsest mesh
              {
                _ctx_vec[i-1].coarser_dm = nullptr;
                _ctx_vec[i-1].finer_dm   = _dms[1].get();
              }
            else if( i == n_levels ) // were at the finest mesh
              {
                _ctx_vec[i-1].coarser_dm = _dms[_dms.size() - 2].get();
                _ctx_vec[i-1].finer_dm   = nullptr;
              }
            else // were in the middle of the hierarchy
              {
                _ctx_vec[i-1].coarser_dm = _dms[i-2].get();
                _ctx_vec[i-1].finer_dm   = _dms[i].get();
              }
          }
      } // End context creation

    // Attach a vector and context to each DM
    if ( n_levels >= 1 )
      {
        for ( unsigned int i = 1; i <= n_levels ; ++i)
          {
            DM & dm = this->get_dm(i-1);

            LibmeshPetscCall2(system.comm(), DMShellSetGlobalVector( dm, _ctx_vec[i-1].current_vec->vec() ));

            LibmeshPetscCall2(system.comm(), DMShellSetContext( dm, &_ctx_vec[i-1] ));
          }
      }

    // DM structures created, now we need projection matrixes if GMG is requested.
    // To prepare for projection creation go to second coarsest mesh so we can utilize
    // old_dof_indices information in the projection creation.
    if (n_levels > 1 )
      {
        // First, stash the coarse dof indices for FS purposes
        unsigned int n_vars  = system.n_vars();
        _ctx_vec[0].dof_vec.resize(n_vars);

        for( unsigned int v = 0; v < n_vars; v++ )
          {
            std::vector<numeric_index_type> di;
            system.get_dof_map().local_variable_indices(di, system.get_mesh(), v);
            _ctx_vec[0].dof_vec[v] = di;
          }

        LOG_CALL ("PDM_refine", "PetscDMWrapper", mesh_refinement.uniformly_refine(1));
        LOG_CALL ("PDM_dist_dof", "PetscDMWrapper", system.get_dof_map().distribute_dofs(mesh));
        LOG_CALL ("PDM_cnstrnts", "PetscDMWrapper", system.reinit_constraints());
        LOG_CALL ("PDM_prep_send", "PetscDMWrapper", system.get_dof_map().prepare_send_list());
      }

    // Create the Interpolation Matrices between adjacent mesh levels
    for ( unsigned int i = 1 ; i < n_levels ; ++i )
      {
        if ( i != n_levels )
          {
            // Stash the rest of the dof indices
            unsigned int n_vars  = system.n_vars();
            _ctx_vec[i].dof_vec.resize(n_vars);

            for( unsigned int v = 0; v < n_vars; v++ )
              {
                std::vector<numeric_index_type> di;
                system.get_dof_map().local_variable_indices(di, system.get_mesh(), v);
                _ctx_vec[i].dof_vec[v] = di;
              }

            unsigned int ndofs_c = _mesh_dof_sizes[i-1];
            unsigned int ndofs_f = _mesh_dof_sizes[i];

            // Create the Interpolation matrix and set its pointer
            _ctx_vec[i-1].K_interp_ptr = _pmtx_vec[i-1].get();
            _ctx_vec[i-1].K_sub_interp_ptr = _subpmtx_vec[i-1].get();

            unsigned int ndofs_local     = system.get_dof_map().n_local_dofs();
            unsigned int ndofs_old_first = system.get_dof_map().first_old_dof();
            unsigned int ndofs_old_end   = system.get_dof_map().end_old_dof();
            unsigned int ndofs_old_size  = ndofs_old_end - ndofs_old_first;

            // Init and zero the matrix
            _ctx_vec[i-1].K_interp_ptr->init(ndofs_f, ndofs_c, ndofs_local, ndofs_old_size, 30 , 20);

            // Disable Mat destruction since PETSc destroys these for us
            _ctx_vec[i-1].K_interp_ptr->set_destroy_mat_on_exit(false);
            _ctx_vec[i-1].K_sub_interp_ptr->set_destroy_mat_on_exit(false);

            // TODO: Projection matrix sparsity pattern?
            //MatSetOption(_ctx_vec[i-1].K_interp_ptr->mat(), MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

            // Compute the interpolation matrix and set K_interp_ptr
            LOG_CALL ("PDM_proj_mat", "PetscDMWrapper", system.projection_matrix(*_ctx_vec[i-1].K_interp_ptr));

            // Always close matrix that contains altered data
            _ctx_vec[i-1].K_interp_ptr->close();
          }

        // Move to next grid to make next projection
        if ( i != n_levels - 1 )
          {
            LOG_CALL ("PDM_refine", "PetscDMWrapper", mesh_refinement.uniformly_refine(1));
            LOG_CALL ("PDM_dist_dof", "PetscDMWrapper", system.get_dof_map().distribute_dofs(mesh));
            LOG_CALL ("PDM_cnstrnts", "PetscDMWrapper", system.reinit_constraints());
            LOG_CALL ("PDM_prep_send", "PetscDMWrapper", system.get_dof_map().prepare_send_list());
          }
      } // End create transfer operators. System back at the finest grid

    return n_levels;
  }

set_point_range_in_section()

void libMesh::PetscDMWrapper::set_point_range_in_section ( const System &  system, PetscSection &  section, 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  )

private

Helper function for build_section.

This function will count how many "points" on the current processor have DoFs associated with them and give that count to PETSc. We need to cache a mapping between the global node id and our local count that we do in this function because we will need the local number again in the add_dofs_to_section function.

Definition at line 794 of file petsc_dm_wrapper.C.

References libMesh::ParallelObject::comm(), libMesh::DofMap::dof_indices(), libMesh::FEType::family, libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::libmesh_assert(), libMesh::DofMap::local_index(), libMesh::make_range(), mesh, libMesh::DofMap::n_local_dofs(), libMesh::ParallelObject::n_processors(), libMesh::DofMap::n_SCALAR_dofs(), libMesh::System::n_vars(), libMesh::System::number(), libMesh::ParallelObject::processor_id(), libMesh::SCALAR, libMesh::Variable::type(), and libMesh::System::variable().

Referenced by build_section().

  {
    // We're expecting this to be empty coming in
    libmesh_assert(node_map.empty());

    // We need to count up the number of active "points" on this processor.
    // Nominally, a "point" in PETSc parlance is a geometric object that can
    // hold DoFs, i.e node, edge, face, elem. Since we handle the mesh and are only
    // interested in solvers, then the only thing PETSc needs is a unique *local* number
    // for each "point" that has active DoFs; note however this local numbering
    // we construct must be continuous.
    //
    // In libMesh, for most finite elements, we just associate those DoFs with the
    // geometric nodes. So can we loop over the nodes on this processor and check
    // if any of the fields are have active DoFs on that node.
    // If so, then we tell PETSc about that "point". At this stage, we just need
    // to count up how many active "points" we have and cache the local number to global id
    // mapping.


    // These will be our local counters. pstart should always be zero.
    // pend will track our local "point" count.
    // If we're on a processor who coarsened the mesh to have no local elements,
    // we should make an empty PetscSection. An empty PetscSection is specified
    // by passing [0,0) to the PetscSectionSetChart call at the end. So, if we
    // have nothing on this processor, these are the correct values to pass to
    // PETSc.
    dof_id_type pstart = 0;
    dof_id_type pend = 0;

    const MeshBase & mesh = system.get_mesh();

    const DofMap & dof_map = system.get_dof_map();

    // If we don't have any local dofs, then there's nothing to tell to the PetscSection
    if (dof_map.n_local_dofs() > 0)
      {
        // Conservative estimate of space needed so we don't thrash
        node_map.reserve(mesh.n_local_nodes());
        elem_map.reserve(mesh.n_active_local_elem());

        // We loop over active elements and then cache the global/local node mapping to make sure
        // we only count active nodes. For example, if we're calling this function and we're
        // not the finest level in the Mesh tree, we don't want to include nodes of child elements
        // that aren't active on this level.
        for (const auto & elem : mesh.active_local_element_ptr_range())
          {
            for (const Node & node : elem->node_ref_range())
              {
                // get the global id number of local node n

                // Only register nodes with the PetscSection if they have dofs that belong to
                // this processor. Even though we're active local elements, the dofs associated
                // with the node may belong to a different processor. The processor who owns
                // those dofs will register that node with the PetscSection on that processor.
                std::vector<dof_id_type> node_dof_indices;
                dof_map.dof_indices( &node, node_dof_indices );
                if( !node_dof_indices.empty() && dof_map.local_index(node_dof_indices[0]) )
                  {
#ifndef NDEBUG
                    // We're assuming that if the first dof for this node belongs to this processor,
                    // then all of them do.
                    for( auto dof : node_dof_indices )
                      libmesh_assert(dof_map.local_index(dof));
#endif
                    // Cache the global/local mapping if we haven't already
                    // Then increment our local count
                    dof_id_type node_id = node.id();
                    if( node_map.count(node_id) == 0 )
                      {
                        node_map.emplace(node_id, pend);
                        pend++;
                      }
                  }
              }

            // Some finite elements, e.g. Hierarchic, associate element interior DoFs with the element
            // rather than the node (since we ought to be able to use Hierachic elements on a QUAD4,
            // which has no interior node). Thus, we also need to check element interiors for DoFs
            // as well and, if the finite element has them, we also need to count the Elem in our
            // "point" accounting.
            if( elem->n_dofs(system.number()) > 0 )
              {
                dof_id_type elem_id = elem->id();
                elem_map.emplace(elem_id, pend);
                pend++;
              }
          }

        // SCALAR dofs live on the "last" processor, so only work there if there are any
        if( dof_map.n_SCALAR_dofs() > 0 && (system.processor_id() == (system.n_processors()-1)) )
          {
            // Loop through all the variables and cache the scalar ones. We cache the
            // SCALAR variable index along with the local point to make it easier when
            // we have to register dofs with the PetscSection
            for (auto v : make_range(system.n_vars()))
              {
                if( system.variable(v).type().family == SCALAR )
                  {
                    scalar_map.emplace(pend, v);
                    pend++;
                  }
              }
          }

      }

    LibmeshPetscCall2(system.comm(), PetscSectionSetChart(section, pstart, pend));
  }

Member Data Documentation

_ctx_vec

std::vector<PetscDMContext> libMesh::PetscDMWrapper::_ctx_vec

private

Vector of internal PetscDM context structs for all grid levels Pointers to these C++ objects are passed to DMShellSetContext(), they are cleaned up automatically by their destructors.

Definition at line 162 of file petsc_dm_wrapper.h.

Referenced by clear(), init_dm_data(), and init_petscdm().

_dms

std::vector<WrappedPetsc<DM> > libMesh::PetscDMWrapper::_dms

private

Vector of DMs for all grid levels.

These are PETSc objects created by calling DMShellCreate(), so we are responsible for cleaning them up.

Definition at line 122 of file petsc_dm_wrapper.h.

Referenced by clear(), get_dm(), init_dm_data(), and init_petscdm().

_mesh_dof_loc_sizes

std::vector<unsigned int> libMesh::PetscDMWrapper::_mesh_dof_loc_sizes

private

Stores n_local_dofs for each grid level, to be used for projection vector sizing.

Definition at line 175 of file petsc_dm_wrapper.h.

Referenced by clear(), init_dm_data(), and init_petscdm().

_mesh_dof_sizes

std::vector<unsigned int> libMesh::PetscDMWrapper::_mesh_dof_sizes

private

Stores n_dofs for each grid level, to be used for projection matrix sizing.

Definition at line 172 of file petsc_dm_wrapper.h.

Referenced by clear(), init_dm_data(), and init_petscdm().

_pmtx_vec

std::vector<std::unique_ptr<PetscMatrixBase<Number> > > libMesh::PetscDMWrapper::_pmtx_vec

private

Vector of projection matrixes for all grid levels.

These are C++ objects, they are cleaned up automatically by their destructors.

Definition at line 148 of file petsc_dm_wrapper.h.

Referenced by clear(), init_dm_data(), and init_petscdm().

_sections

std::vector<WrappedPetsc<PetscSection> > libMesh::PetscDMWrapper::_sections

private

Vector of PETScSections for all grid levels.

These are PETSc objects which are attached to the DM by calling DMSetLocalSection(). While the DM takes care of destroying existing PetscSections in calls to DMSetLocalSection(), it does not appear to clean up PetscSections it holds when it is destroyed, so we manage their lifetimes using WrappedPetsc objects.

Definition at line 132 of file petsc_dm_wrapper.h.

Referenced by clear(), get_section(), and init_dm_data().

_star_forests

std::vector<PetscSF> libMesh::PetscDMWrapper::_star_forests

private

Vector of star forests for all grid levels.

These are PETSc objects which are attached to the DM by calling DMSetSectionSF(). The DM seems to take care of cleaning these up itself as far as I can tell, so we do not try to manage their lifetime in any way.

Definition at line 141 of file petsc_dm_wrapper.h.

Referenced by clear(), get_star_forest(), and init_dm_data().

_subpmtx_vec

std::vector<std::unique_ptr<PetscMatrixBase<Number> > > libMesh::PetscDMWrapper::_subpmtx_vec

private

Vector of sub projection matrixes for all grid levels for fieldsplit.

These are C++ objects, they are cleaned up automatically by their destructors.

Definition at line 155 of file petsc_dm_wrapper.h.

Referenced by clear(), init_dm_data(), and init_petscdm().

_vec_vec

std::vector<std::unique_ptr<PetscVector<Number> > > libMesh::PetscDMWrapper::_vec_vec

private

Vector of solution vectors for all grid levels.

These are C++ objects, they are cleaned up automatically by their destructors.

Definition at line 169 of file petsc_dm_wrapper.h.

Referenced by clear(), init_dm_data(), and init_petscdm().

---

The documentation for this class was generated from the following files:

• include/solvers/petsc_dm_wrapper.h
• src/solvers/petsc_dm_wrapper.C