Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add 'requires' clauses to PETSc solvers. #15132

Merged
merged 3 commits into from
May 8, 2023
Merged

Add 'requires' clauses to PETSc solvers. #15132

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
81c5479
Add 'requires' clauses to PETSc solvers.
bangerth Apr 22, 2023
9927c90
Also allow vector types that are can just wrap Vec objects.
bangerth Apr 23, 2023
feb909a
Also annotate the implementation of functions with 'requires' clauses.
bangerth Apr 23, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
47 changes: 47 additions & 0 deletions include/deal.II/base/template_constraints.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -656,6 +656,7 @@ namespace LinearAlgebra
#ifdef DEAL_II_WITH_PETSC
namespace PETScWrappers
{
class VectorBase;
class Vector;
class BlockVector;

Expand Down Expand Up @@ -788,6 +789,40 @@ namespace concepts
dealii::LinearAlgebra::TpetraWrappers::Vector<Number>> = true;
# endif
# endif


/**
* A template variable that returns whether the template argument is
* a valid deal.II vector type that is internally built on PETSc
* functionality. Its general definition is `false`, with
* specializations dealing with actual vector types for which the
* predicate is `true`.
*/
template <typename T>
inline constexpr bool is_dealii_petsc_vector_type = false;

# ifdef DEAL_II_WITH_PETSC
template <>
inline constexpr bool
is_dealii_petsc_vector_type<dealii::PETScWrappers::VectorBase> = true;

template <>
inline constexpr bool
is_dealii_petsc_vector_type<dealii::PETScWrappers::Vector> = true;

template <>
inline constexpr bool
is_dealii_petsc_vector_type<dealii::PETScWrappers::BlockVector> = true;

template <>
inline constexpr bool
is_dealii_petsc_vector_type<dealii::PETScWrappers::MPI::Vector> = true;

template <>
inline constexpr bool
is_dealii_petsc_vector_type<dealii::PETScWrappers::MPI::BlockVector> =
true;
# endif
} // namespace internal


Expand Down Expand Up @@ -820,6 +855,18 @@ namespace concepts
(std::is_const_v<VectorType> ==
false);

/**
* A concept that tests whether a given template argument is a deal.II
* vector type that internally builds on PETSc functionality. This
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't be better to rename this to can_wrap_petsc_vector_type or something?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is to reflect the fact that the supported classes must have the ability to construct new instances by wrapping a Vec and can return the Vec via the petsc_vector() method.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same for the matrices

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A naive question: can this be done using traits? or with enable_if checking for explicit type::type(Vec) and Vec type::petsc_vector()?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, this can be done. The most explicit approach would be to do something like this:

template <typename VectorType>
concept can_wrap_petsc_vec =
  requires {
    VectorType vec_object (std::declval<Vec>());  // checks that an object can be constructed from a Vec object
  };

and then to use this concept in the same place as where I'm currently already using the other concept.

But it turns out that others have found such a concept useful in the past, and so it is also std::constructible_from<VectorType,Vec> :-)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks. So the whole point of std::concepts is to make the error message from the compilation stage readable? (and not a million lines long?)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's basically the idea. For example, if you have

  template <typename MeshType>
  int n_cells (const MeshType &m) {
    return m.size();
  }

then you will get awkward error messages if you call that function with an int argument -- something like error: 'int' variable 'm' does not have a member function 'size()'.

So you want to annotate the function with a "concept" like

  template <typename MeshType>
    requires is_mesh_type<MeshType>
  int n_cells (const MeshType &m) {
    return m.size();
  }

that encodes the kinds of properties the type MeshType has to have to make this function useful. If you call n_cells(1), the compiler will say something like error: can't call 'n_cells()' with argument 'MeshType=int' because 'int' does not satisfy the concept 'is_mesh_type'.

That is how I use concepts in this patch. There are other uses for concepts, but that's a different story. The idea for these other uses are of the kind where you want to use different algorithms based on properties of the involved types. Say something like this:

  template <typename MatrixType>
    requires (matrix_type_represents_symmetric_matrix<MatrixType> == false)
  std::pair<FullMatrix,FullMatrix>
  lu_decomposition (const MatrixType &m)
  {
    ...use LU algorithm...
  }

  template <typename MatrixType>
    requires (matrix_type_represents_symmetric_matrix<MatrixType> == true)
  std::pair<FullMatrix,FullMatrix>
  lu_decomposition (const MatrixType &m)
  {
    ...use Cholesky algorithm...
  }

Basically, concepts here serve the purpose of "dispatching" to overloaded versions of the same function name, based on whether or not the type satisfies certain properties or not. This was previously possible with things such as std::enable_if, but it was awkward at best.

* concept is used to constrain some classes that implement advanced
* functionality based on PETSc and that requires that the vector
* it works on are PETSc vectors. This includes, for example, the
* time stepping and nonlinear solver classes in namespace PETScWrappers.
*/
template <typename VectorType>
concept is_dealii_petsc_vector_type =
internal::is_dealii_petsc_vector_type<VectorType>;

#endif
} // namespace concepts

Expand Down
7 changes: 7 additions & 0 deletions include/deal.II/lac/petsc_snes.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,11 @@

# include <petscsnes.h>

# if defined(DEAL_II_HAVE_CXX20)
# include <concepts>
# endif


DEAL_II_NAMESPACE_OPEN

namespace PETScWrappers
Expand Down Expand Up @@ -224,6 +229,8 @@ namespace PETScWrappers
template <typename VectorType = PETScWrappers::VectorBase,
typename PMatrixType = PETScWrappers::MatrixBase,
typename AMatrixType = PMatrixType>
DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
std::constructible_from<VectorType, Vec>))
class NonlinearSolver
{
public:
Expand Down
62 changes: 40 additions & 22 deletions include/deal.II/lac/petsc_snes.templates.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,8 @@ DEAL_II_NAMESPACE_OPEN
namespace PETScWrappers
{
template <typename VectorType, typename PMatrixType, typename AMatrixType>
DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
std::constructible_from<VectorType, Vec>))
NonlinearSolver<VectorType, PMatrixType, AMatrixType>::NonlinearSolver(
const NonlinearSolverData &data,
const MPI_Comm & mpi_comm)
Expand All @@ -63,6 +65,8 @@ namespace PETScWrappers


template <typename VectorType, typename PMatrixType, typename AMatrixType>
DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
std::constructible_from<VectorType, Vec>))
NonlinearSolver<VectorType, PMatrixType, AMatrixType>::operator SNES() const
{
return snes;
Expand All @@ -71,25 +75,29 @@ namespace PETScWrappers


template <typename VectorType, typename PMatrixType, typename AMatrixType>
SNES
NonlinearSolver<VectorType, PMatrixType, AMatrixType>::petsc_snes()
DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
std::constructible_from<VectorType, Vec>))
SNES NonlinearSolver<VectorType, PMatrixType, AMatrixType>::petsc_snes()
{
return snes;
}



template <typename VectorType, typename PMatrixType, typename AMatrixType>
MPI_Comm
NonlinearSolver<VectorType, PMatrixType, AMatrixType>::get_mpi_communicator()
const
DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
std::constructible_from<VectorType, Vec>))
MPI_Comm NonlinearSolver<VectorType, PMatrixType, AMatrixType>::
get_mpi_communicator() const
{
return PetscObjectComm(reinterpret_cast<PetscObject>(snes));
}



template <typename VectorType, typename PMatrixType, typename AMatrixType>
DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
std::constructible_from<VectorType, Vec>))
NonlinearSolver<VectorType, PMatrixType, AMatrixType>::~NonlinearSolver()
{
AssertPETSc(SNESDestroy(&snes));
Expand All @@ -98,17 +106,19 @@ namespace PETScWrappers


template <typename VectorType, typename PMatrixType, typename AMatrixType>
void
NonlinearSolver<VectorType, PMatrixType, AMatrixType>::reinit()
DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
std::constructible_from<VectorType, Vec>))
void NonlinearSolver<VectorType, PMatrixType, AMatrixType>::reinit()
{
AssertPETSc(SNESReset(snes));
}



template <typename VectorType, typename PMatrixType, typename AMatrixType>
void
NonlinearSolver<VectorType, PMatrixType, AMatrixType>::reinit(
DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
std::constructible_from<VectorType, Vec>))
void NonlinearSolver<VectorType, PMatrixType, AMatrixType>::reinit(
const NonlinearSolverData &data)
{
reinit();
Expand Down Expand Up @@ -154,8 +164,9 @@ namespace PETScWrappers


template <typename VectorType, typename PMatrixType, typename AMatrixType>
void
NonlinearSolver<VectorType, PMatrixType, AMatrixType>::set_matrix(
DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
std::constructible_from<VectorType, Vec>))
void NonlinearSolver<VectorType, PMatrixType, AMatrixType>::set_matrix(
PMatrixType &P)
{
this->A = nullptr;
Expand All @@ -165,8 +176,9 @@ namespace PETScWrappers


template <typename VectorType, typename PMatrixType, typename AMatrixType>
void
NonlinearSolver<VectorType, PMatrixType, AMatrixType>::set_matrices(
DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
std::constructible_from<VectorType, Vec>))
void NonlinearSolver<VectorType, PMatrixType, AMatrixType>::set_matrices(
AMatrixType &A,
PMatrixType &P)
{
Expand All @@ -177,8 +189,10 @@ namespace PETScWrappers


template <typename VectorType, typename PMatrixType, typename AMatrixType>
unsigned int
NonlinearSolver<VectorType, PMatrixType, AMatrixType>::solve(VectorType &x)
DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
std::constructible_from<VectorType, Vec>))
unsigned int NonlinearSolver<VectorType, PMatrixType, AMatrixType>::solve(
VectorType &x)
{
const auto snes_function =
[](SNES, Vec x, Vec f, void *ctx) -> PetscErrorCode {
Expand Down Expand Up @@ -393,9 +407,11 @@ namespace PETScWrappers


template <typename VectorType, typename PMatrixType, typename AMatrixType>
unsigned int
NonlinearSolver<VectorType, PMatrixType, AMatrixType>::solve(VectorType & x,
PMatrixType &P)
DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
std::constructible_from<VectorType, Vec>))
unsigned int NonlinearSolver<VectorType, PMatrixType, AMatrixType>::solve(
VectorType & x,
PMatrixType &P)
{
set_matrix(P);
return solve(x);
Expand All @@ -404,10 +420,12 @@ namespace PETScWrappers


template <typename VectorType, typename PMatrixType, typename AMatrixType>
unsigned int
NonlinearSolver<VectorType, PMatrixType, AMatrixType>::solve(VectorType & x,
AMatrixType &A,
PMatrixType &P)
DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
std::constructible_from<VectorType, Vec>))
unsigned int NonlinearSolver<VectorType, PMatrixType, AMatrixType>::solve(
VectorType & x,
AMatrixType &A,
PMatrixType &P)
{
set_matrices(A, P);
return solve(x);
Expand Down
6 changes: 6 additions & 0 deletions include/deal.II/lac/petsc_ts.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,10 @@

# include <petscts.h>

# if defined(DEAL_II_HAVE_CXX20)
# include <concepts>
# endif

DEAL_II_NAMESPACE_OPEN

namespace PETScWrappers
Expand Down Expand Up @@ -289,6 +293,8 @@ namespace PETScWrappers
template <typename VectorType = PETScWrappers::VectorBase,
typename PMatrixType = PETScWrappers::MatrixBase,
typename AMatrixType = PMatrixType>
DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
std::constructible_from<VectorType, Vec>))
class TimeStepper
{
public:
Expand Down