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)

Private Member Functions
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. More...
void	add_dofs_helper (const System &system, const DofObject &dof_object, dof_id_type local_id, PetscSection &section)
	Helper function to reduce code duplication when setting dofs in section. More...

Private Attributes
std::vector< std::unique_ptr< DM > >	_dms
	Vector of DMs for all grid levels. More...
std::vector< std::unique_ptr< PetscSection > >	_sections
	Vector of PETScSections for all grid levels. More...
std::vector< std::unique_ptr< PetscSF > >	_star_forests
	Vector of star forests for all grid levels. More...
std::vector< std::unique_ptr< PetscMatrix< Number > > >	_pmtx_vec
	Vector of projection matrixes for all grid levels. More...
std::vector< std::unique_ptr< PetscMatrix< Number > > >	_subpmtx_vec
	Vector of sub projection matrixes for all grid levels for fieldsplit. More...
std::vector< std::unique_ptr< PetscDMContext > >	_ctx_vec
	Vector of internal PetscDM context structs for all grid levels. 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 95 of file petsc_dm_wrapper.h.

Constructor & Destructor Documentation

PetscDMWrapper()

libMesh::PetscDMWrapper::PetscDMWrapper ( )

default

~PetscDMWrapper()

libMesh::PetscDMWrapper::~PetscDMWrapper

(

)

Definition at line 425 of file petsc_dm_wrapper.C.

  {
    this->clear();
  }

References clear().

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

Helper function to reduce code duplication when setting dofs in section.

Definition at line 1039 of file petsc_dm_wrapper.C.

  {
    unsigned int total_n_dofs_at_dofobject = 0;
 
    // We are assuming variables are also numbered 0 to n_vars()-1
    for (auto v : IntRange<unsigned int>(0, system.n_vars()))
      {
        unsigned int n_dofs_at_dofobject = dof_object.n_dofs(system.number(), v);
 
        if( n_dofs_at_dofobject > 0 )
          {
            PetscErrorCode ierr = PetscSectionSetFieldDof( section,
                                                           local_id,
                                                           v,
                                                           n_dofs_at_dofobject );
 
            CHKERRABORT(system.comm().get(),ierr);
 
            total_n_dofs_at_dofobject += n_dofs_at_dofobject;
          }
      }
 
    libmesh_assert_equal_to(total_n_dofs_at_dofobject, dof_object.n_dofs(system.number()));
 
    PetscErrorCode ierr =
      PetscSectionSetDof( section, local_id, total_n_dofs_at_dofobject );
    CHKERRABORT(system.comm().get(),ierr);
  }

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

Referenced by add_dofs_to_section().

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 972 of file petsc_dm_wrapper.C.

  {
    const MeshBase & mesh = system.get_mesh();
 
    PetscErrorCode ierr;
 
    // 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 & nmap : node_map )
      {
        const dof_id_type global_node_id = nmap.first;
        const dof_id_type local_node_id = nmap.second;
 
        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 & emap : elem_map )
      {
        const dof_id_type global_elem_id = emap.first;
        const dof_id_type local_elem_id = emap.second;
 
        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 & smap : scalar_map )
          {
            const dof_id_type local_id = smap.first;
            const unsigned int scalar_var = smap.second;
 
            // The number of SCALAR dofs comes from the variable order
            const int n_dofs = system.variable(scalar_var).type().order.get_order();
 
            ierr = PetscSectionSetFieldDof( section, local_id, scalar_var, n_dofs );
            CHKERRABORT(system.comm().get(),ierr);
 
            // 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
            ierr = PetscSectionSetDof( section, local_id, n_dofs );
            CHKERRABORT(system.comm().get(),ierr);
          }
      }
 
  }

References add_dofs_helper(), libMesh::ParallelObject::comm(), libMesh::System::get_mesh(), libMesh::OrderWrapper::get_order(), libMesh::ierr, 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().

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 729 of file petsc_dm_wrapper.C.

  {
    START_LOG ("build_section()", "PetscDMWrapper");
 
    PetscErrorCode ierr;
    ierr = PetscSectionCreate(system.comm().get(),&section);
    CHKERRABORT(system.comm().get(),ierr);
 
    // Tell the PetscSection about all of our System variables
    ierr = PetscSectionSetNumFields(section,system.n_vars());
    CHKERRABORT(system.comm().get(),ierr);
 
    // Set the actual names of all the field variables
    for (auto v : IntRange<unsigned int>(0, system.n_vars()))
      {
        ierr = PetscSectionSetFieldName( section, v, system.variable_name(v).c_str() );
        CHKERRABORT(system.comm().get(),ierr);
      }
 
    // 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
    if (!libMesh::on_command_line("--node-major-dofs"))
      libmesh_error_msg("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);
    //  }
 
    ierr = PetscSectionSetUp(section);
    CHKERRABORT(system.comm().get(),ierr);
 
    // Sanity checking at least that local_n_dofs match
    libmesh_assert_equal_to(system.n_local_dofs(),this->check_section_n_dofs(section));
 
    STOP_LOG ("build_section()", "PetscDMWrapper");
  }

References add_dofs_to_section(), check_section_n_dofs(), libMesh::ParallelObject::comm(), libMesh::ierr, 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_and_attach_petscdm().

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 795 of file petsc_dm_wrapper.C.

  {
    START_LOG ("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()));
 
    // 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();
 
    PetscErrorCode ierr;
    ierr = PetscSFCreate(system.comm().get(), &star_forest);
    CHKERRABORT(system.comm().get(),ierr);
 
    // 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.
    ierr = PetscSFSetGraph(star_forest,
                           n_roots,
                           n_leaves,
                           local_dofs,
                           PETSC_COPY_VALUES,
                           remote_dofs,
                           PETSC_COPY_VALUES);
    CHKERRABORT(system.comm().get(),ierr);
 
    STOP_LOG ("build_sf()", "PetscDMWrapper");
  }

References libMesh::ParallelObject::comm(), libMesh::DofMap::dof_owner(), libMesh::DofMap::first_dof(), libMesh::System::get_dof_map(), libMesh::DofMap::get_send_list(), libMesh::ierr, libMesh::index_range(), libMesh::DofMap::n_local_dofs(), and PETSC_COPY_VALUES.

Referenced by init_and_attach_petscdm().

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 1072 of file petsc_dm_wrapper.C.

  {
    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;
  }

References libMesh::ierr.

Referenced by build_section().

clear()

void libMesh::PetscDMWrapper::clear

(

)

Destroys and clears all build DM-related data.

Definition at line 430 of file petsc_dm_wrapper.C.

  {
    // PETSc will destroy the attached PetscSection, PetscSF as well as
    // other relateds such as the Projections so we just tidy up the
    // containers here.
 
    _dms.clear();
    _sections.clear();
    _star_forests.clear();
    _pmtx_vec.clear();
    _vec_vec.clear();
    _ctx_vec.clear();
 
  }

References _ctx_vec, _dms, _pmtx_vec, _sections, _star_forests, and _vec_vec.

Referenced by libMesh::PetscDiffSolver::clear(), and ~PetscDMWrapper().

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 144 of file petsc_dm_wrapper.h.

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

References _dms.

Referenced by init_and_attach_petscdm().

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 152 of file petsc_dm_wrapper.h.

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

References _sections.

Referenced by init_and_attach_petscdm().

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 160 of file petsc_dm_wrapper.h.

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

References _star_forests.

Referenced by init_and_attach_petscdm().

init_and_attach_petscdm()

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

Definition at line 445 of file petsc_dm_wrapper.C.

  {
    START_LOG ("init_and_attach_petscdm()", "PetscDMWrapper");
 
    PetscErrorCode ierr;
 
    MeshBase & mesh = system.get_mesh();   // Convenience
    MeshRefinement mesh_refinement(mesh); // Used for swapping between grids
 
    // Theres 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 shouldnt make projections
    if (n_levels >= 1)
      n_levels += 1;
 
    // 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
        ierr = DMShellCreate(system.comm().get(), &dm);
        CHKERRABORT(system.comm().get(),ierr);
 
        // Set the DM embedding dimension to help PetscDS (Discrete System)
        ierr = DMSetCoordinateDim(dm, mesh.mesh_dimension());
        CHKERRABORT(system.comm().get(),ierr);
 
        // Build the PetscSection and attach it to the DM
        this->build_section(system, section);
#if PETSC_VERSION_LESS_THAN(3,12,0)
        ierr = DMSetDefaultSection(dm, section);
#else
        ierr = DMSetSection(dm, section);
#endif
        CHKERRABORT(system.comm().get(),ierr);
 
        // 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)
            ierr = DMSetDefaultSF(dm, star_forest);
#else
            ierr = DMSetSectionSF(dm, star_forest);
#endif
            CHKERRABORT(system.comm().get(),ierr);
          }
 
        // Set PETSC's Restriction, Interpolation, Coarsen and Refine functions for the current DM
        ierr = DMShellSetCreateInterpolation ( dm, libmesh_petsc_DMCreateInterpolation );
        CHKERRABORT(system.comm().get(),ierr);
 
        // Not implemented. For now we rely on galerkin style restrictions
        bool supply_restriction = false;
        if (supply_restriction)
          {
            ierr = DMShellSetCreateRestriction ( dm, libmesh_petsc_DMCreateRestriction  );
            CHKERRABORT(system.comm().get(),ierr);
          }
 
        ierr = DMShellSetCoarsen ( dm, libmesh_petsc_DMCoarsen );
        CHKERRABORT(system.comm().get(),ierr);
 
        ierr = DMShellSetRefine ( dm, libmesh_petsc_DMRefine );
        CHKERRABORT(system.comm().get(),ierr);
 
        ierr= DMShellSetCreateSubDM(dm, libmesh_petsc_DMCreateSubDM);
        CHKERRABORT(system.comm().get(), ierr);
 
        // Uniformly coarsen if not the coarsest grid and distribute dof info.
        if ( level != 1 )
          {
            START_LOG ("PDM_coarsen", "PetscDMWrapper");
            mesh_refinement.uniformly_coarsen(1);
            STOP_LOG  ("PDM_coarsen", "PetscDMWrapper");
 
            START_LOG ("PDM_dist_dof", "PetscDMWrapper");
            system.get_dof_map().distribute_dofs(mesh);
            STOP_LOG  ("PDM_dist_dof", "PetscDMWrapper");
          }
      } // 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);
 
            ierr = DMShellSetGlobalVector( dm, (*_ctx_vec[ i-1 ]).current_vec->vec() );
            CHKERRABORT(system.comm().get(),ierr);
 
            ierr = DMShellSetContext( dm, _ctx_vec[ i-1 ].get() );
            CHKERRABORT(system.comm().get(),ierr);
          }
      }
 
    // 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;
          }
 
        START_LOG ("PDM_refine", "PetscDMWrapper");
        mesh_refinement.uniformly_refine(1);
        STOP_LOG  ("PDM_refine", "PetscDMWrapper");
 
        START_LOG ("PDM_dist_dof", "PetscDMWrapper");
        system.get_dof_map().distribute_dofs(mesh);
        STOP_LOG  ("PDM_dist_dof", "PetscDMWrapper");
 
        START_LOG ("PDM_cnstrnts", "PetscDMWrapper");
        system.reinit_constraints();
        STOP_LOG ("PDM_cnstrnts", "PetscDMWrapper");
      }
 
    // 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_dofs_on_processor(system.processor_id());
            unsigned int ndofs_old_first = system.get_dof_map().first_old_dof(system.processor_id());
            unsigned int ndofs_old_end   = system.get_dof_map().end_old_dof(system.processor_id());
            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
            START_LOG ("PDM_proj_mat", "PetscDMWrapper");
            system.projection_matrix(*_ctx_vec[i-1]->K_interp_ptr);
            STOP_LOG  ("PDM_proj_mat", "PetscDMWrapper");
 
            // 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 )
          {
            START_LOG ("PDM_refine", "PetscDMWrapper");
            mesh_refinement.uniformly_refine(1);
            STOP_LOG  ("PDM_refine", "PetscDMWrapper");
 
            START_LOG ("PDM_dist_dof", "PetscDMWrapper");
            system.get_dof_map().distribute_dofs(mesh);
            STOP_LOG ("PDM_dist_dof", "PetscDMWrapper");
 
            START_LOG ("PDM_cnstrnts", "PetscDMWrapper");
            system.reinit_constraints();
            STOP_LOG  ("PDM_cnstrnts", "PetscDMWrapper");
 
          }
      } // End create transfer operators. System back at the finest grid
 
    // Lastly, give SNES the finest level DM
    DM & dm = this->get_dm(n_levels-1);
    ierr = SNESSetDM(snes, dm);
    CHKERRABORT(system.comm().get(),ierr);
 
    STOP_LOG ("init_and_attach_petscdm()", "PetscDMWrapper");
  }

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::DofMap::end_old_dof(), libMesh::DofMap::first_old_dof(), libMesh::ReferenceElem::get(), get_dm(), libMesh::System::get_dof_map(), libMesh::System::get_mesh(), get_section(), get_star_forest(), libMesh::ierr, 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::DofMap::n_dofs_on_processor(), libMesh::MeshTools::n_levels(), libMesh::DofMap::n_local_dofs(), libMesh::ParallelObject::n_processors(), n_vars, libMesh::System::n_vars(), libMesh::ParallelObject::processor_id(), libMesh::System::projection_matrix(), libMesh::System::reinit_constraints(), libMesh::MeshRefinement::uniformly_coarsen(), and libMesh::MeshRefinement::uniformly_refine().

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

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 1096 of file petsc_dm_wrapper.C.

  {
    _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++ )
      {
        _dms[i] = libmesh_make_unique<DM>();
        _sections[i] = libmesh_make_unique<PetscSection>();
        _star_forests[i] = libmesh_make_unique<PetscSF>();
        _ctx_vec[i] = libmesh_make_unique<PetscDMContext>();
        _pmtx_vec[i]= libmesh_make_unique<PetscMatrix<Number>>(comm);
        _subpmtx_vec[i]= libmesh_make_unique<PetscMatrix<Number>>(comm);
        _vec_vec[i] = libmesh_make_unique<PetscVector<Number>>(comm);
      }
  }

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_and_attach_petscdm().

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 857 of file petsc_dm_wrapper.C.

  {
    // 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.insert(std::make_pair(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.insert(std::make_pair(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 : IntRange<unsigned int>(0, system.n_vars()))
              {
                if( system.variable(v).type().family == SCALAR )
                  {
                    scalar_map.insert(std::make_pair(pend,v));
                    pend++;
                  }
              }
          }
 
      }
 
    PetscErrorCode ierr = PetscSectionSetChart(section, pstart, pend);
    CHKERRABORT(system.comm().get(),ierr);
  }

References libMesh::ParallelObject::comm(), libMesh::DofMap::dof_indices(), libMesh::FEType::family, libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::ierr, libMesh::libmesh_assert(), libMesh::DofMap::local_index(), 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().

Member Data Documentation

_ctx_vec

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

private

Vector of internal PetscDM context structs for all grid levels.

Definition at line 126 of file petsc_dm_wrapper.h.

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

_dms

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

private

Vector of DMs for all grid levels.

Definition at line 111 of file petsc_dm_wrapper.h.

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

_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 135 of file petsc_dm_wrapper.h.

Referenced by init_and_attach_petscdm(), and init_dm_data().

_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 132 of file petsc_dm_wrapper.h.

Referenced by init_and_attach_petscdm(), and init_dm_data().

_pmtx_vec

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

private

Vector of projection matrixes for all grid levels.

Definition at line 120 of file petsc_dm_wrapper.h.

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

_sections

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

private

Vector of PETScSections for all grid levels.

Definition at line 114 of file petsc_dm_wrapper.h.

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

_star_forests

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

private

Vector of star forests for all grid levels.

Definition at line 117 of file petsc_dm_wrapper.h.

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

_subpmtx_vec

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

private

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

Definition at line 123 of file petsc_dm_wrapper.h.

Referenced by init_and_attach_petscdm(), and init_dm_data().

_vec_vec

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

private

Vector of solution vectors for all grid levels.

Definition at line 129 of file petsc_dm_wrapper.h.

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

---

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

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