This class provides an interface to the suite of preconditioners available from PETSc. More...

#include <petsc_preconditioner.h>

Inheritance diagram for libMesh::PetscPreconditioner< T >:

Public Member Functions
	PetscPreconditioner (const libMesh::Parallel::Communicator &comm_in)
	Constructor.

virtual	~PetscPreconditioner ()=default

virtual void	apply (const NumericVector< T > &x, NumericVector< T > &y) override
	Computes the preconditioned vector `y` based on input vector `x`.

virtual void	clear () override
	Release all memory and clear data structures.

virtual void	init () override
	Initialize data structures if not done so already.

PC	pc ()

virtual bool	initialized () const

virtual void	setup ()
	This is called every time the "operator might have changed".

virtual void	zero ()
	Can be used to zero items relevant to the preconditioner.

void	set_matrix (SparseMatrix< Number > &mat)
	Sets the matrix to be preconditioned.

PreconditionerType	type () const

void	set_type (const PreconditionerType pct)
	Sets the type of preconditioner to use.

const Parallel::Communicator &	comm () const

processor_id_type	n_processors () const

processor_id_type	processor_id () const

Static Public Member Functions
static void	set_petsc_preconditioner_type (const PreconditionerType &preconditioner_type, PC &pc)
	Tells PETSc to use the user-specified preconditioner.

static void	set_petsc_aux_data (PC &pc, System &sys, const unsigned v=0)
	Builds PETSc auxiliary data needed by preconditioners such as hypre ams/ads.

static void	set_hypre_ams_data (PC &pc, System &sys, const unsigned v)

static void	set_hypre_ads_data (PC &pc, System &sys, const unsigned v)

static void	set_petsc_aux_data (PC &, System &, const unsigned=0)

static std::unique_ptr< Preconditioner< T > >	build_preconditioner (const libMesh::Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
	Builds a `Preconditioner` using the linear solver package specified by `solver_package`, returning the result wrapped in a std::unique_ptr for safety.

static std::string	get_info ()
	Gets a string containing the reference information.

static void	print_info (std::ostream &out_stream=libMesh::out)
	Prints the reference information, by default to `libMesh::out`.

static unsigned int	n_objects ()
	Prints the number of outstanding (created, but not yet destroyed) objects.

static void	enable_print_counter_info ()
	Methods to enable/disable the reference counter output from print_info().

static void	disable_print_counter_info ()

Protected Types
typedef std::map< std::string, std::pair< unsigned int, unsigned int > >	Counts
	Data structure to log the information.

Protected Member Functions
void	increment_constructor_count (const std::string &name) noexcept
	Increments the construction counter.

void	increment_destructor_count (const std::string &name) noexcept
	Increments the destruction counter.

Protected Attributes
WrappedPetsc< PC >	_pc
	Preconditioner context.

Mat	_mat
	PETSc Mat pulled out of the _matrix object during init().

SparseMatrix< T > *	_matrix
	The matrix P... ie the matrix to be preconditioned.

PreconditionerType	_preconditioner_type
	Enum stating with type of preconditioner to use.

bool	_is_initialized
	Flag indicating if the data structures have been initialized.

const Parallel::Communicator &	_communicator

Static Protected Attributes
static Counts	_counts
	Actually holds the data.

static Threads::atomic< unsigned int >	_n_objects
	The number of objects.

static Threads::spin_mutex	_mutex
	Mutual exclusion object to enable thread-safe reference counting.

static bool	_enable_print_counter = true
	Flag to control whether reference count information is printed when print_info is called.

Detailed Description

template<typename T>
class libMesh::PetscPreconditioner< T >

This class provides an interface to the suite of preconditioners available from PETSc.

All overridden virtual functions are documented in preconditioner.h.

Author: Derek Gaston

Date: 2009

Definition at line 57 of file petsc_preconditioner.h.

Member Typedef Documentation

◆ Counts

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > libMesh::ReferenceCounter::Counts

protectedinherited

Data structure to log the information.

The log is identified by the class name.

Definition at line 119 of file reference_counter.h.

Constructor & Destructor Documentation

◆ PetscPreconditioner()

template<typename T >

libMesh::PetscPreconditioner< T >::PetscPreconditioner ( const libMesh::Parallel::Communicator & comm_in )

Constructor.

Initializes PetscPreconditioner data structures

Definition at line 37 of file petsc_preconditioner.C.

                                                                                        :
  Preconditioner<T>(comm_in)
{}

◆ ~PetscPreconditioner()

template<typename T >

virtual libMesh::PetscPreconditioner< T >::~PetscPreconditioner ( )

virtualdefault

Member Function Documentation

◆ apply()

template<typename T >

void libMesh::PetscPreconditioner< T >::apply	(	const NumericVector< T > &	x,
		NumericVector< T > &	y
	)

overridevirtual

Computes the preconditioned vector y based on input vector x.

This is usually done by solving \( Py=x \) to get the action of \( P^-1 x \).

Implements libMesh::Preconditioner< T >.

Definition at line 44 of file petsc_preconditioner.C.

{
  PetscVector<T> & x_pvec = cast_ref<PetscVector<T> &>(const_cast<NumericVector<T> &>(x));
  PetscVector<T> & y_pvec = cast_ref<PetscVector<T> &>(const_cast<NumericVector<T> &>(y));
 
  Vec x_vec = x_pvec.vec();
  Vec y_vec = y_pvec.vec();
 
  LibmeshPetscCall(PCApply(_pc, x_vec, y_vec));
}

References libMesh::PetscVector< T >::vec().

◆ build_preconditioner()

template<typename T >

std::unique_ptr< Preconditioner< T > > libMesh::Preconditioner< T >::build_preconditioner	(	const libMesh::Parallel::Communicator &	comm,
		const SolverPackage	solver_package = `libMesh::default_solver_package()`
	)

staticinherited

Builds a Preconditioner using the linear solver package specified by solver_package, returning the result wrapped in a std::unique_ptr for safety.

Definition at line 50 of file preconditioner.C.

{
  // Avoid unused parameter warnings when no solver packages are enabled.
  libmesh_ignore(comm);
 
  // Build and return the appropriate Preconditioner object.
  switch (solver_package)
    {
 
#ifdef LIBMESH_HAVE_PETSC
    case PETSC_SOLVERS:
      {
        return std::make_unique<PetscPreconditioner<T>>(comm);
      }
#endif
 
#ifdef LIBMESH_TRILINOS_HAVE_EPETRA
    case TRILINOS_SOLVERS:
      return std::make_unique<TrilinosPreconditioner<T>>(comm);
#endif
 
#ifdef LIBMESH_HAVE_EIGEN
    case EIGEN_SOLVERS:
      return std::make_unique<EigenPreconditioner<T>>(comm);
#endif
 
    default:
      libmesh_error_msg("ERROR:  Unrecognized solver package: " << solver_package);
    }
}

References libMesh::libmesh_ignore(), and libMesh::PETSC_SOLVERS.

◆ clear()

template<typename T >

void libMesh::PetscPreconditioner< T >::clear ( )

overridevirtual

Release all memory and clear data structures.

Reimplemented from libMesh::Preconditioner< T >.

Definition at line 93 of file petsc_preconditioner.C.

{
  // Calls custom deleter
  _pc.destroy();
}

◆ comm()

const Parallel::Communicator & libMesh::ParallelObject::comm ( ) const

inlineinherited

Returns: A reference to the Parallel::Communicator object used by this mesh.

Definition at line 97 of file parallel_object.h.

98 { return _communicator; }

libMesh::ParallelObject::_communicator

const Parallel::Communicator & _communicator

Definition parallel_object.h:120

References libMesh::ParallelObject::_communicator.

◆ disable_print_counter_info()

void libMesh::ReferenceCounter::disable_print_counter_info ( )

staticinherited

Definition at line 100 of file reference_counter.C.

{
  _enable_print_counter = false;
  return;
}

References libMesh::ReferenceCounter::_enable_print_counter.

◆ enable_print_counter_info()

void libMesh::ReferenceCounter::enable_print_counter_info ( )

staticinherited

Methods to enable/disable the reference counter output from print_info().

Enabled by default.

Definition at line 94 of file reference_counter.C.

{
  _enable_print_counter = true;
  return;
}

References libMesh::ReferenceCounter::_enable_print_counter.

Referenced by libMesh::LibMeshInit::~LibMeshInit().

◆ get_info()

std::string libMesh::ReferenceCounter::get_info ( )

staticinherited

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
 
  std::ostringstream oss;
 
  oss << '\n'
      << " ---------------------------------------------------------------------------- \n"
      << "| Reference count information                                                |\n"
      << " ---------------------------------------------------------------------------- \n";
 
  for (const auto & [name, cd] : _counts)
    oss << "| " << name << " reference count information:\n"
        << "|  Creations:    " << cd.first    << '\n'
        << "|  Destructions: " << cd.second << '\n';
 
  oss << " ---------------------------------------------------------------------------- \n";
 
  return oss.str();
 
#else
 
  return "";
 
#endif
}

References libMesh::ReferenceCounter::_counts.

Referenced by libMesh::ReferenceCounter::print_info().

◆ increment_constructor_count()

void libMesh::ReferenceCounter::increment_constructor_count ( const std::string & name )

inlineprotectednoexceptinherited

Increments the construction counter.

Should be called in the constructor of any derived class that will be reference counted.

Definition at line 183 of file reference_counter.h.

{
  libmesh_try
  {
    Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
    std::pair<unsigned int, unsigned int> & p = _counts[name];
    p.first++;
  }
  libmesh_catch (...)
  {
    auto stream = libMesh::err.get();
    stream->exceptions(stream->goodbit); // stream must not throw
    libMesh::err << "Encountered unrecoverable error while calling "
                 << "ReferenceCounter::increment_constructor_count() "
                 << "for a(n) " << name << " object." << std::endl;
    std::terminate();
  }
}

References libMesh::err, libMesh::BasicOStreamProxy< charT, traits >::get(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< T >::ReferenceCountedObject(), libMesh::ReferenceCountedObject< T >::ReferenceCountedObject(), and libMesh::ReferenceCountedObject< T >::ReferenceCountedObject().

◆ increment_destructor_count()

void libMesh::ReferenceCounter::increment_destructor_count ( const std::string & name )

inlineprotectednoexceptinherited

Increments the destruction counter.

Should be called in the destructor of any derived class that will be reference counted.

Definition at line 207 of file reference_counter.h.

{
  libmesh_try
  {
    Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
    std::pair<unsigned int, unsigned int> & p = _counts[name];
    p.second++;
  }
  libmesh_catch (...)
  {
    auto stream = libMesh::err.get();
    stream->exceptions(stream->goodbit); // stream must not throw
    libMesh::err << "Encountered unrecoverable error while calling "
                 << "ReferenceCounter::increment_destructor_count() "
                 << "for a(n) " << name << " object." << std::endl;
    std::terminate();
  }
}

References libMesh::err, libMesh::BasicOStreamProxy< charT, traits >::get(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< T >::~ReferenceCountedObject().

◆ init()

template<typename T >

void libMesh::PetscPreconditioner< T >::init ( )

overridevirtual

Initialize data structures if not done so already.

Note: This MUST be called before the preconditioning object is used.

Reimplemented from libMesh::Preconditioner< T >.

Definition at line 59 of file petsc_preconditioner.C.

{
  libmesh_error_msg_if(!this->_matrix, "ERROR: No matrix set for PetscPreconditioner, but init() called");
 
  // Clear the preconditioner in case it has been created in the past
  if (!this->_is_initialized)
    {
      // Should probably use PCReset(), but it's not working at the moment so we'll destroy instead
      if (_pc)
        _pc.destroy();
 
      LibmeshPetscCall(PCCreate(this->comm().get(), _pc.get()));
 
      auto pmatrix = cast_ptr<PetscMatrixBase<T> *>(this->_matrix);
      _mat = pmatrix->mat();
    }
 
  LibmeshPetscCall(PCSetOperators(_pc, _mat, _mat));
 
  // Set the PCType.  Note: this used to be done *before* the call to
  // PCSetOperators(), and only when !_is_initialized, but
  // 1.) Some preconditioners (those employing sub-preconditioners,
  // for example) have to call PCSetUp(), and can only do this after
  // the operators have been set.
  // 2.) It should be safe to call set_petsc_preconditioner_type()
  // multiple times.
  set_petsc_preconditioner_type(this->_preconditioner_type, *_pc);
 
  this->_is_initialized = true;
}

◆ initialized()

template<typename T >

virtual bool libMesh::Preconditioner< T >::initialized ( ) const

inlinevirtualinherited

Returns: true if the data structures are initialized, false otherwise.

Reimplemented in libMesh::StaticCondensationPreconditioner.

Definition at line 84 of file preconditioner.h.

84{ return _is_initialized; }

References libMesh::Preconditioner< T >::_is_initialized.

Referenced by libMesh::libmesh_petsc_preconditioner_setup().

◆ n_objects()

static unsigned int libMesh::ReferenceCounter::n_objects ( )

inlinestaticinherited

Prints the number of outstanding (created, but not yet destroyed) objects.

Definition at line 85 of file reference_counter.h.

86 { return _n_objects; }

libMesh::ReferenceCounter::_n_objects

static Threads::atomic< unsigned int > _n_objects

The number of objects.

Definition reference_counter.h:132

References libMesh::ReferenceCounter::_n_objects.

Referenced by libMesh::LibMeshInit::~LibMeshInit().

◆ n_processors()

processor_id_type libMesh::ParallelObject::n_processors ( ) const

inlineinherited

Returns: The number of processors in the group.

Definition at line 103 of file parallel_object.h.

  {
    processor_id_type returnval =
      cast_int<processor_id_type>(_communicator.size());
    libmesh_assert(returnval); // We never have an empty comm
    return returnval;
  }

References libMesh::ParallelObject::_communicator, libMesh::libmesh_assert(), and libMesh::Parallel::Communicator::size().

◆ pc()

template<typename T >

PC libMesh::PetscPreconditioner< T >::pc ( )

Returns: The PETSc PC object. Can be useful for implementing more advanced algorithms.

Definition at line 102 of file petsc_preconditioner.C.

{
  return _pc;
}

◆ print_info()

void libMesh::ReferenceCounter::print_info ( std::ostream & out_stream = libMesh::out )

staticinherited

Prints the reference information, by default to libMesh::out.

Definition at line 81 of file reference_counter.C.

{
  if (_enable_print_counter)
    out_stream << ReferenceCounter::get_info();
}

References libMesh::ReferenceCounter::_enable_print_counter, and libMesh::ReferenceCounter::get_info().

Referenced by libMesh::LibMeshInit::~LibMeshInit().

◆ processor_id()

processor_id_type libMesh::ParallelObject::processor_id ( ) const

inlineinherited

Returns: The rank of this processor in the group.

Definition at line 114 of file parallel_object.h.

115 { return cast_int<processor_id_type>(_communicator.rank()); }

libMesh::Parallel::Communicator::rank

processor_id_type rank() const

References libMesh::ParallelObject::_communicator, and libMesh::Parallel::Communicator::rank().

◆ set_hypre_ads_data()

template<typename T >

void libMesh::PetscPreconditioner< T >::set_hypre_ads_data	(	PC &	pc,
		System &	sys,
		const unsigned	v
	)

static

Definition at line 355 of file petsc_preconditioner.C.

{
  // Get the communicator from the PETSc object
  Parallel::communicator comm;
  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)pc, &comm);
  libmesh_error_msg_if(ierr != LIBMESH_PETSC_SUCCESS,
                       "Error retrieving communicator");
  Parallel::Communicator Comm(comm);
 
  // The mesh/problem dimension
  const unsigned dim = sys.get_mesh().mesh_dimension();
 
  // Dummy Lagrange and Nédélec systems defined over the same mesh so we can
  // enumerate the vertices and edges, respectively
  System & lagrange_sys = sys.get_equation_systems().add_system<System>("__hypre_ads_vertices");
  lagrange_sys.hide_output() = true;
  lagrange_sys.add_variable("__lagrange");
  lagrange_sys.reinit_mesh();
  System & nedelec_sys = sys.get_equation_systems().add_system<System>("__hypre_ads_edges");
  nedelec_sys.hide_output() = true;
  nedelec_sys.add_variable("__nedelec", FIRST, NEDELEC_ONE);
  nedelec_sys.reinit_mesh();
 
  // Global (i.e. total) and local (i.e. to this processor) number of faces,
  // edges and vertices
  const std::vector<dof_id_type> n_all_faces = sys.get_dof_map().n_dofs_per_processor(v);
  const dof_id_type n_glb_faces = std::accumulate(n_all_faces.begin(), n_all_faces.end(), 0);
  const dof_id_type n_loc_faces = n_all_faces[global_processor_id()];
  const dof_id_type n_glb_edges = nedelec_sys.n_dofs();
  const dof_id_type n_loc_edges = nedelec_sys.n_local_dofs();
  const dof_id_type n_glb_verts = lagrange_sys.n_dofs();
  const dof_id_type n_loc_verts = lagrange_sys.n_local_dofs();
 
  // We require local indexing through the vertices and contiguous indexing
  // through the edges and faces: since the vertices are enumerated on their own
  // system, a local index can simply be obtained by subtracting those vertices
  // in all preceding processors; since the edges are enumerated on their own
  // system, we already have contiguous indexing; for the faces, if the system
  // we are given by the user has multiple variables/splits, we effectively need
  // to add local face numbers to those faces in all preceding processors
 
  // The number of vertices and faces in all preceding processors
  dof_id_signed_type vert_offset = -lagrange_sys.get_dof_map().first_dof();
  dof_id_signed_type face_offset =
    std::accumulate(n_all_faces.begin(), n_all_faces.begin() + global_processor_id(), 0);
 
  // Whether dofs are in variable-major or node-major order
  bool var_major = !libMesh::on_command_line("--node-major-dofs");
 
  // If variable-major order, we need only subtract all the preceding dofs; but
  // if node-major order, we need to enumerate the dofs on this processor
  std::vector<dof_id_type> idx_faces;
  if (var_major)
  {
    face_offset -= sys.get_dof_map().first_dof();
    for (auto i : make_range(v))
      face_offset -= sys.get_dof_map().n_local_dofs(i);
  }
  else
    sys.get_dof_map().local_variable_indices(idx_faces, sys.get_mesh(), v);
 
  // Create the discrete grandient matrix, representing the edges in terms of its vertices
  // Preallocate 2 diagonal + 2 off-diagonal nonzeros as the vertices could fall on either
  PetscMatrix<Real> G(Comm, n_glb_edges, n_glb_verts, n_loc_edges, n_loc_verts, 2, 2);
 
  // Create the discrete curl matrix, representing the faces in terms of its edges
  // Preallocate 4 diagonal + 4 off-diagonal nonzeros as the edges could fall on either
  PetscMatrix<Real> C(Comm, n_glb_faces, n_glb_edges, n_loc_faces, n_loc_edges, 4, 4);
 
  // Create array for the coordinates of the local vertices
  PetscReal * coords;
  LibmeshPetscCallA(comm, PetscMalloc1(dim * n_loc_verts, &coords));
 
  // Populate the discrete gradient matrix and the coordinates array
  for (const auto & elem : sys.get_mesh().active_local_element_ptr_range())
    for (auto face : make_range(elem->n_faces()))
    {
      // The number of edges on this face
      const unsigned n_face_edges = Elem::type_to_n_sides_map[elem->side_type(face)];
 
      // The faces's three/four edges
      std::vector<dof_id_type> edge_dofs(n_face_edges);
      std::vector<bool> edge_orients(n_face_edges);
      for (auto face_edge : make_range(n_face_edges))
      {
        // Convert from face-wise to element-wise edge id
        const unsigned edge = elem->local_side_node(face, n_face_edges + face_edge)
                            - elem->n_vertices();
        libmesh_assert(elem->is_edge_on_side(edge, face));
 
        // The edge's two vertices
        dof_id_type vert_dofs[2];
        for (auto vert : make_range(2))
        {
          const unsigned loc_vert_node_id = elem->local_edge_node(edge, vert);
          libmesh_assert(elem->is_vertex(loc_vert_node_id));
 
          const Node & vert_node = elem->node_ref(loc_vert_node_id);
          vert_dofs[vert] = vert_node.dof_number(lagrange_sys.number(), 0, 0);
 
          // If owned, populate coordinates array
          if (vert_node.processor_id() == global_processor_id())
          {
            const dof_id_type loc_vert_dof = vert_offset + vert_dofs[vert];
            libmesh_assert(loc_vert_dof <  n_loc_verts);
 
            for (auto d : make_range(dim))
              coords[dim * loc_vert_dof + d] = vert_node(d);
          }
        }
 
        // The edge's (middle) node
        const unsigned loc_edge_node_id = elem->local_edge_node(edge, 2);
        libmesh_assert(elem->is_edge(loc_edge_node_id));
 
        const Node & edge_node = elem->node_ref(loc_edge_node_id);
        edge_dofs[face_edge] = edge_node.dof_number(nedelec_sys.number(), 0, 0);
        edge_orients[face_edge] = elem->positive_edge_orientation(edge) ^
                                  elem->relative_edge_face_order(edge, face);
 
        // If owned, populate discrete gradient matrix
        if (edge_node.processor_id() == global_processor_id())
        {
          const Real sign = elem->positive_edge_orientation(edge) ? 1 : -1;
 
          libmesh_assert(edge_dofs[face_edge] >= G.row_start());
          libmesh_assert(edge_dofs[face_edge] <  G.row_stop());
          G.set(edge_dofs[face_edge], vert_dofs[0],  sign);
          G.set(edge_dofs[face_edge], vert_dofs[1], -sign);
        }
      }
 
      // The faces's (middle) node
      const unsigned loc_face_node_id = elem->local_side_node(face, 2 * n_face_edges);
      libmesh_assert(elem->is_face(loc_face_node_id));
 
      const Node & face_node = elem->node_ref(loc_face_node_id);
 
      // If owned, populate discrete curl matrix
      if (face_node.processor_id() == global_processor_id())
      {
        const dof_id_type face_dof = face_node.dof_number(sys.number(), v, 0);
 
        const dof_id_type cont_face_dof = face_offset +
          (var_major ? face_dof
                     : std::distance(idx_faces.begin(),
                       std::find(idx_faces.begin(), idx_faces.end(), face_dof)));
 
        const bool face_orient = elem->positive_face_orientation(face);
 
        for (auto face_edge : make_range(n_face_edges))
        {
          const Real sign = face_orient ^ edge_orients[face_edge] ? 1 : -1;
 
          libmesh_assert(cont_face_dof >= C.row_start());
          libmesh_assert(cont_face_dof <  C.row_stop());
          C.set(cont_face_dof, edge_dofs[face_edge], sign);
        }
      }
    }
 
  // Assemble the discrete gradient and discrete curl matrices
  G.close();
  C.close();
 
#ifndef NDEBUG
  // The product CG of the two matrices should be the zero matrix
  PetscMatrix<Real> CG(Comm);
  C.matrix_matrix_mult(G, CG);
  libmesh_assert(CG.linfty_norm() == 0);
#endif
 
  // Hand over the matrices and coordinates array
  LibmeshPetscCallA(comm, PCHYPRESetDiscreteGradient(pc, G.mat()));
  LibmeshPetscCallA(comm, PCHYPRESetDiscreteCurl(pc, C.mat()));
  LibmeshPetscCallA(comm, PCSetCoordinates(pc, dim, n_loc_verts, coords));
 
  // Free allocated memory for the coordinates array
  LibmeshPetscCallA(comm, PetscFree(coords));
}

References libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), libMesh::CG, libMesh::PetscMatrixBase< T >::close(), TIMPI::communicator, dim, libMesh::DofObject::dof_number(), libMesh::FIRST, libMesh::DofMapBase::first_dof(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::System::get_mesh(), libMesh::global_processor_id(), libMesh::System::hide_output(), libMesh::libmesh_assert(), libMesh::DofMap::local_variable_indices(), libMesh::make_range(), libMesh::PetscMatrixBase< T >::mat(), libMesh::PetscMatrix< T >::matrix_matrix_mult(), libMesh::MeshBase::mesh_dimension(), libMesh::System::n_dofs(), libMesh::DofMap::n_dofs_per_processor(), libMesh::System::n_local_dofs(), libMesh::DofMap::n_local_dofs(), libMesh::NEDELEC_ONE, libMesh::System::number(), libMesh::on_command_line(), libMesh::DofObject::processor_id(), libMesh::Real, libMesh::System::reinit_mesh(), libMesh::PetscMatrixBase< T >::row_start(), libMesh::PetscMatrixBase< T >::row_stop(), libMesh::PetscMatrix< T >::set(), and libMesh::Elem::type_to_n_sides_map.

◆ set_hypre_ams_data()

template<typename T >

void libMesh::PetscPreconditioner< T >::set_hypre_ams_data	(	PC &	pc,
		System &	sys,
		const unsigned	v
	)

static

Definition at line 231 of file petsc_preconditioner.C.

{
  // Get the communicator from the PETSc object
  Parallel::communicator comm;
  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)pc, &comm);
  libmesh_error_msg_if(ierr != LIBMESH_PETSC_SUCCESS,
                       "Error retrieving communicator");
  Parallel::Communicator Comm(comm);
 
  // The mesh/problem dimension
  const unsigned dim = sys.get_mesh().mesh_dimension();
 
  // Dummy Lagrange system defined over the same mesh so we can enumerate the vertices
  System & lagrange_sys = sys.get_equation_systems().add_system<System>("__hypre_ams_vertices");
  lagrange_sys.hide_output() = true;
  lagrange_sys.add_variable("__lagrange");
  lagrange_sys.reinit_mesh();
 
  // Global (i.e. total) and local (i.e. to this processor) number of edges and vertices
  const std::vector<dof_id_type> n_all_edges = sys.get_dof_map().n_dofs_per_processor(v);
  const dof_id_type n_glb_edges = std::accumulate(n_all_edges.begin(), n_all_edges.end(), 0);
  const dof_id_type n_loc_edges = n_all_edges[global_processor_id()];
  const dof_id_type n_glb_verts = lagrange_sys.n_dofs();
  const dof_id_type n_loc_verts = lagrange_sys.n_local_dofs();
 
  // We require local indexing through the vertices and contiguous indexing
  // through the edges: since the vertices are enumerated on their own system,
  // a local index can simply be obtained by subtracting those vertices in all
  // preceding processors; for the edges, if the system we are given by the
  // user has multiple variables/splits, we effectively need to add local edge
  // numbers to those edges in all preceding processors
 
  // The number of vertices and edges in all preceding processors
  dof_id_signed_type vert_offset = -lagrange_sys.get_dof_map().first_dof();
  dof_id_signed_type edge_offset =
    std::accumulate(n_all_edges.begin(), n_all_edges.begin() + global_processor_id(), 0);
 
  // Whether dofs are in variable-major or node-major order
  bool var_major = !libMesh::on_command_line("--node-major-dofs");
 
  // If variable-major order, we need only subtract all the preceding dofs; but
  // if node-major order, we need to enumerate the dofs on this processor
  std::vector<dof_id_type> idx_edges;
  if (var_major)
  {
    edge_offset -= sys.get_dof_map().first_dof();
    for (auto i : make_range(v))
      edge_offset -= sys.get_dof_map().n_local_dofs(i);
  }
  else
    sys.get_dof_map().local_variable_indices(idx_edges, sys.get_mesh(), v);
 
  // Create the discrete grandient matrix, representing the edges in terms of its vertices
  // Preallocate 2 diagonal + 2 off-diagonal nonzeros as the vertices could fall on either
  PetscMatrix<Real> G(Comm, n_glb_edges, n_glb_verts, n_loc_edges, n_loc_verts, 2, 2);
 
  // Create array for the coordinates of the local vertices
  PetscReal * coords;
  LibmeshPetscCallA(comm, PetscMalloc1(dim * n_loc_verts, &coords));
 
  // Populate the discrete gradient matrix and the coordinates array
  for (const auto & elem : sys.get_mesh().active_local_element_ptr_range())
    for (auto edge : make_range(elem->n_edges()))
    {
      // The edge's two vertices
      dof_id_type vert_dofs[2];
      for (auto vert : make_range(2))
      {
        const unsigned loc_vert_node_id = elem->local_edge_node(edge, vert);
        libmesh_assert(elem->is_vertex(loc_vert_node_id));
 
        const Node & vert_node = elem->node_ref(loc_vert_node_id);
        vert_dofs[vert] = vert_node.dof_number(lagrange_sys.number(), 0, 0);
 
        // If owned, populate coordinates array
        if (vert_node.processor_id() == global_processor_id())
        {
          const dof_id_type loc_vert_dof = vert_offset + vert_dofs[vert];
          libmesh_assert(loc_vert_dof <  n_loc_verts);
 
          for (auto d : make_range(dim))
            coords[dim * loc_vert_dof + d] = vert_node(d);
        }
      }
 
      // The edge's (middle) node
      const unsigned loc_edge_node_id = elem->local_edge_node(edge, 2);
      libmesh_assert(elem->is_edge(loc_edge_node_id));
 
      const Node & edge_node = elem->node_ref(loc_edge_node_id);
 
      // If owned, populate discrete gradient matrix
      if (edge_node.processor_id() == global_processor_id())
      {
        const dof_id_type edge_dof = edge_node.dof_number(sys.number(), v, 0);
 
        const dof_id_type cont_edge_dof = edge_offset +
          (var_major ? edge_dof
                     : std::distance(idx_edges.begin(),
                       std::find(idx_edges.begin(), idx_edges.end(), edge_dof)));
 
        const Real sign = elem->positive_edge_orientation(edge) ? 1 : -1;
 
        libmesh_assert(cont_edge_dof >= G.row_start());
        libmesh_assert(cont_edge_dof <  G.row_stop());
        G.set(cont_edge_dof, vert_dofs[0],  sign);
        G.set(cont_edge_dof, vert_dofs[1], -sign);
      }
    }
 
  // Assemble the discrete gradient matrix
  G.close();
 
  // Hand over the matrix and coordinates array
  LibmeshPetscCallA(comm, PCHYPRESetDiscreteGradient(pc, G.mat()));
  LibmeshPetscCallA(comm, PCSetCoordinates(pc, dim, n_loc_verts, coords));
 
  // Free allocated memory for the coordinates array
  LibmeshPetscCallA(comm, PetscFree(coords));
}

References libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), libMesh::PetscMatrixBase< T >::close(), TIMPI::communicator, dim, libMesh::DofObject::dof_number(), libMesh::DofMapBase::first_dof(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::System::get_mesh(), libMesh::global_processor_id(), libMesh::System::hide_output(), libMesh::libmesh_assert(), libMesh::DofMap::local_variable_indices(), libMesh::make_range(), libMesh::PetscMatrixBase< T >::mat(), libMesh::MeshBase::mesh_dimension(), libMesh::System::n_dofs(), libMesh::DofMap::n_dofs_per_processor(), libMesh::System::n_local_dofs(), libMesh::DofMap::n_local_dofs(), libMesh::System::number(), libMesh::on_command_line(), libMesh::DofObject::processor_id(), libMesh::Real, libMesh::System::reinit_mesh(), libMesh::PetscMatrixBase< T >::row_start(), libMesh::PetscMatrixBase< T >::row_stop(), and libMesh::PetscMatrix< T >::set().

◆ set_matrix()

template<typename T >

void libMesh::Preconditioner< T >::set_matrix ( SparseMatrix< Number > & mat )

inherited

Sets the matrix to be preconditioned.

Definition at line 164 of file preconditioner.h.

{
  //If the matrix is changing then we (probably) need to reinitialize.
  _is_initialized = false;
  _matrix = &mat;
}

◆ set_petsc_aux_data() [1/2]

template<typename T >

static void libMesh::PetscPreconditioner< T >::set_petsc_aux_data	(	PC &	,
		System &	,
		const unsigned	= `0`
	)

inlinestatic

Definition at line 93 of file petsc_preconditioner.h.

93{}

◆ set_petsc_aux_data() [2/2]

template<typename T >

void libMesh::PetscPreconditioner< T >::set_petsc_aux_data	(	PC &	pc,
		System &	sys,
		const unsigned	v = `0`
	)

static

Builds PETSc auxiliary data needed by preconditioners such as hypre ams/ads.

Definition at line 161 of file petsc_preconditioner.C.

{
  // Get the communicator from the PETSc object
  Parallel::communicator comm;
  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)pc, &comm);
  libmesh_error_msg_if(ierr != LIBMESH_PETSC_SUCCESS,
                       "Error retrieving communicator");
 
  // Make sure the preconditioner options are set
  LibmeshPetscCallA(comm, PCSetFromOptions(pc));
 
  // Get the type of preconditioner we are using
  PCType pc_type = nullptr;
  LibmeshPetscCallA(comm, PCGetType(pc, &pc_type));
 
  // Check if hypre ams/ads, otherwise we quit with nothing to do
  if (pc_type && std::string(pc_type) == PCHYPRE)
  {
    // Get the hypre preconditioner we are using
    PCType hypre_type = nullptr;
    LibmeshPetscCallA(comm, PCHYPREGetType(pc, &hypre_type));
 
    // If not ams/ads, we quit with nothing to do
    if (std::string(hypre_type) == "ams")
    {
      // If multiple variables, we error out as senseless
      libmesh_error_msg_if(sys.n_vars() > 1,
                           "Error applying hypre AMS to a system with multiple "
                           "variables");
      // If not a 1st order Nédélec or a 2d 1st order Raviart-Thomas system, we
      // error out as we do not support anything else at the moment
      const FEType & var_type = sys.variable(v).type();
      const bool first_order_nedelec =
        var_type.order == FIRST && var_type.family == NEDELEC_ONE;
      const bool first_order_raviart_thomas =
        var_type.order == FIRST && var_type.family == RAVIART_THOMAS;
 
      libmesh_error_msg_if(!first_order_nedelec &&
                           (!first_order_raviart_thomas ||
                            sys.get_mesh().mesh_dimension() != 2),
                           "Error applying hypre AMS to a system "
                           "whose variable is not 1st order Nedelec or 1st "
                           "order Raviart-Thomas on a 2d mesh");
      set_hypre_ams_data(pc, sys, v);
    }
    else if (std::string(hypre_type) == "ads")
    {
      // If multiple variables, we error out as senseless
      libmesh_error_msg_if(sys.n_vars() > 1,
                           "Error applying hypre ADS to a system with multiple "
                           "variables");
      // If not a 3d 1st order Raviart-Thomas system, we error out as we do not
      // support anything else at the moment
      const FEType & var_type = sys.variable(v).type();
      const bool first_order_raviart_thomas =
        var_type.order == FIRST && var_type.family == RAVIART_THOMAS;
 
      libmesh_error_msg_if(!first_order_raviart_thomas ||
                           sys.get_mesh().mesh_dimension() != 3,
                           "Error applying hypre ADS to a system "
                           "whose variable is not 1st "
                           "order Raviart-Thomas on a 3d mesh");
      set_hypre_ads_data(pc, sys, v);
    }
  }
}

References TIMPI::communicator, libMesh::FEType::family, libMesh::FIRST, libMesh::System::get_mesh(), libMesh::MeshBase::mesh_dimension(), libMesh::System::n_vars(), libMesh::NEDELEC_ONE, libMesh::FEType::order, libMesh::RAVIART_THOMAS, libMesh::Variable::type(), and libMesh::System::variable().

Referenced by libMesh::PetscLinearSolver< Number >::ksp(), and libMesh::PetscNonlinearSolver< T >::solve().

◆ set_petsc_preconditioner_type()

template<typename T >

void libMesh::PetscPreconditioner< T >::set_petsc_preconditioner_type	(	const PreconditionerType &	preconditioner_type,
		PC &	pc
	)

static

Tells PETSc to use the user-specified preconditioner.

Definition at line 110 of file petsc_preconditioner.C.

{
  // get the communicator from the PETSc object
  Parallel::communicator comm;
  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)pc, & comm);
  if (ierr != LIBMESH_PETSC_SUCCESS)
    libmesh_error_msg("Error retrieving communicator");
 
  #define CasePCSetType(PreconditionerType, PCType)                       \
  case PreconditionerType:                                                \
    LibmeshPetscCallA(comm, PCSetType (pc, const_cast<KSPType>(PCType))); \
    break;
 
  switch (preconditioner_type)
    {
    CasePCSetType(IDENTITY_PRECOND,     PCNONE)
    CasePCSetType(CHOLESKY_PRECOND,     PCCHOLESKY)
    CasePCSetType(ICC_PRECOND,          PCICC)
    CasePCSetType(ILU_PRECOND,          PCILU)
    CasePCSetType(LU_PRECOND,           PCLU)
    CasePCSetType(ASM_PRECOND,          PCASM)
    CasePCSetType(JACOBI_PRECOND,       PCJACOBI)
    CasePCSetType(BLOCK_JACOBI_PRECOND, PCBJACOBI)
    CasePCSetType(SOR_PRECOND,          PCSOR)
    CasePCSetType(EISENSTAT_PRECOND,    PCEISENSTAT)
    CasePCSetType(AMG_PRECOND,          PCHYPRE)
    CasePCSetType(SVD_PRECOND,          PCSVD)
    CasePCSetType(USER_PRECOND,         PCMAT)
    CasePCSetType(SHELL_PRECOND,        PCSHELL)
 
    default:
      libMesh::err << "ERROR:  Unsupported PETSC Preconditioner: "
                   << Utility::enum_to_string(preconditioner_type) << std::endl
                   << "Continuing with PETSC defaults" << std::endl;
    }
 
  // Set additional options if we are doing AMG and
  // HYPRE is available
#ifdef LIBMESH_HAVE_PETSC_HYPRE
  if (preconditioner_type == AMG_PRECOND)
    LibmeshPetscCallA(comm, PCHYPRESetType(pc, "boomeramg"));
#endif
 
  // Let the commandline override stuff
  LibmeshPetscCallA(comm, PCSetFromOptions(pc));
}

References libMesh::AMG_PRECOND, libMesh::ASM_PRECOND, libMesh::BLOCK_JACOBI_PRECOND, libMesh::CHOLESKY_PRECOND, TIMPI::communicator, libMesh::EISENSTAT_PRECOND, libMesh::Utility::enum_to_string(), libMesh::err, libMesh::ICC_PRECOND, libMesh::IDENTITY_PRECOND, libMesh::ILU_PRECOND, libMesh::JACOBI_PRECOND, libMesh::LU_PRECOND, libMesh::SHELL_PRECOND, libMesh::SOR_PRECOND, libMesh::SVD_PRECOND, and libMesh::USER_PRECOND.

◆ set_type()

template<typename T >

void libMesh::Preconditioner< T >::set_type ( const PreconditionerType pct )

inherited

Sets the type of preconditioner to use.

Definition at line 173 of file preconditioner.h.

{
  //If the preconditioner type changes we (probably) need to reinitialize.
  _is_initialized = false;
  _preconditioner_type = pct;
}

◆ setup()

template<typename T >

virtual void libMesh::Preconditioner< T >::setup ( )

inlinevirtualinherited

This is called every time the "operator might have changed".

This is where you need to fill in your preconditioning matrix.

Reimplemented in libMesh::StaticCondensationPreconditioner.

Definition at line 110 of file preconditioner.h.

110{}

Referenced by libMesh::libmesh_petsc_preconditioner_setup().

◆ type()

template<typename T >

PreconditionerType libMesh::Preconditioner< T >::type ( ) const

inlineinherited

Returns: The type of preconditioner to use.

Definition at line 125 of file preconditioner.h.

125{ return _preconditioner_type; }

References libMesh::Preconditioner< T >::_preconditioner_type.

◆ zero()

template<typename T >

virtual void libMesh::Preconditioner< T >::zero ( )

inlinevirtualinherited

Can be used to zero items relevant to the preconditioner.

Reimplemented in libMesh::StaticCondensationPreconditioner.

Definition at line 115 of file preconditioner.h.

115{}

Member Data Documentation

◆ _communicator

const Parallel::Communicator& libMesh::ParallelObject::_communicator

protectedinherited

Definition at line 120 of file parallel_object.h.

Referenced by libMesh::EquationSystems::build_parallel_elemental_solution_vector(), libMesh::EquationSystems::build_parallel_solution_vector(), libMesh::StaticCondensation::close(), libMesh::ParallelObject::comm(), libMesh::CondensedEigenSystem::initialize_condensed_matrices(), libMesh::ParallelObject::n_processors(), libMesh::ParallelObject::operator=(), libMesh::ParallelObject::processor_id(), libMesh::BoundaryInfo::regenerate_id_sets(), libMesh::StaticCondensationDofMap::reinit(), and libMesh::BoundaryInfo::synchronize_global_id_set().

◆ _counts

ReferenceCounter::Counts libMesh::ReferenceCounter::_counts

staticprotectedinherited

Actually holds the data.

Definition at line 124 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::get_info().

◆ _enable_print_counter

bool libMesh::ReferenceCounter::_enable_print_counter = true

staticprotectedinherited

Flag to control whether reference count information is printed when print_info is called.

Definition at line 143 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::disable_print_counter_info(), libMesh::ReferenceCounter::enable_print_counter_info(), and libMesh::ReferenceCounter::print_info().

◆ _is_initialized

template<typename T >

bool libMesh::Preconditioner< T >::_is_initialized

protectedinherited

Flag indicating if the data structures have been initialized.

Definition at line 148 of file preconditioner.h.

Referenced by libMesh::Preconditioner< T >::initialized().

◆ _mat

template<typename T >

Mat libMesh::PetscPreconditioner< T >::_mat

protected

PETSc Mat pulled out of the _matrix object during init().

We aren't responsible for cleaning up this one.

Definition at line 107 of file petsc_preconditioner.h.

◆ _matrix

template<typename T >

SparseMatrix<T>* libMesh::Preconditioner< T >::_matrix

protectedinherited

The matrix P... ie the matrix to be preconditioned.

This is often the actual system matrix of a linear system.

Definition at line 138 of file preconditioner.h.

◆ _mutex

Threads::spin_mutex libMesh::ReferenceCounter::_mutex

staticprotectedinherited

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 137 of file reference_counter.h.

◆ _n_objects

Threads::atomic< unsigned int > libMesh::ReferenceCounter::_n_objects

staticprotectedinherited

The number of objects.

Print the reference count information when the number returns to 0.

Definition at line 132 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::n_objects(), libMesh::ReferenceCounter::ReferenceCounter(), libMesh::ReferenceCounter::ReferenceCounter(), libMesh::ReferenceCounter::ReferenceCounter(), and libMesh::ReferenceCounter::~ReferenceCounter().

◆ _pc

template<typename T >

WrappedPetsc<PC> libMesh::PetscPreconditioner< T >::_pc

protected

Preconditioner context.

Definition at line 101 of file petsc_preconditioner.h.

◆ _preconditioner_type

template<typename T >

PreconditionerType libMesh::Preconditioner< T >::_preconditioner_type

protectedinherited

Enum stating with type of preconditioner to use.

Definition at line 143 of file preconditioner.h.

Referenced by libMesh::Preconditioner< T >::type().

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

include/numerics/petsc_preconditioner.h
src/numerics/petsc_preconditioner.C

Public Member Functions

Static Public Member Functions

Protected Types

Protected Member Functions

Protected Attributes

Static Protected Attributes

Detailed Description

Member Typedef Documentation

◆ Counts

Constructor & Destructor Documentation

◆ PetscPreconditioner()

◆ ~PetscPreconditioner()

Member Function Documentation

◆ apply()

◆ build_preconditioner()

◆ clear()

◆ comm()

◆ disable_print_counter_info()

◆ enable_print_counter_info()

◆ get_info()

◆ increment_constructor_count()

◆ increment_destructor_count()

◆ init()

◆ initialized()

◆ n_objects()

◆ n_processors()

◆ pc()

◆ print_info()

◆ processor_id()

◆ set_hypre_ads_data()

◆ set_hypre_ams_data()

◆ set_matrix()

◆ set_petsc_aux_data() [1/2]

◆ set_petsc_aux_data() [2/2]

◆ set_petsc_preconditioner_type()

◆ set_type()

◆ setup()

◆ type()

◆ zero()

Member Data Documentation

◆ _communicator

◆ _counts

◆ _enable_print_counter

◆ _is_initialized

◆ _mat

◆ _matrix

◆ _mutex

◆ _n_objects

◆ _pc

◆ _preconditioner_type