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.

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

partially address https://github.com/dealii/dealii/issues/14656 #14735

Merged
merged 3 commits into from
Jan 27, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
4 changes: 2 additions & 2 deletions cmake/configure/configure_20_petsc.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ macro(feature_petsc_find_external var)
endif()

#
# Petsc has to be configured with the same MPI configuration as
# PETSc has to be configured with the same MPI configuration as
# deal.II.
#
# petscconf.h should export PETSC_HAVE_MPIUNI 1 in case mpi support is
Expand All @@ -68,7 +68,7 @@ macro(feature_petsc_find_external var)
endif()

#
# Petsc has to be configured with the same number of bits for indices as
# PETSc has to be configured with the same number of bits for indices as
# deal.II.
#
# petscconf.h should export PETSC_WITH_64BIT_INDICES 1 in case 64bits
Expand Down
2 changes: 1 addition & 1 deletion examples/step-44/doc/results.dox
Original file line number Diff line number Diff line change
Expand Up @@ -278,7 +278,7 @@ There are a number of obvious extensions for this work:
- The program has been developed for solving problems on single-node
multicore machines. With a little effort, the program could be
extended to a large-scale computing environment through the use of
Petsc or Trilinos, using a similar technique to that demonstrated in
PETSc or Trilinos, using a similar technique to that demonstrated in
step-40. This would mostly involve changes to the setup, assembly,
<code>PointHistory</code> and linear solver routines.
- As this program assumes quasi-static equilibrium, extensions to
Expand Down
33 changes: 17 additions & 16 deletions include/deal.II/lac/petsc_block_sparse_matrix.h
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,11 @@ namespace PETScWrappers
BlockSparseMatrix() = default;

/**
* Create a BlockSparseMatrix with a PETSc Mat
* Create a BlockSparseMatrix with a PETSc Mat that describes the entire
* block matrix.
* It infers the number of blocks from the Mat if it is of type MATNEST,
* otherwise the block operator will only have a single block.
* Internally, we always store a MATNEST matrix.
*/
explicit BlockSparseMatrix(const Mat &);

Expand Down Expand Up @@ -173,6 +177,14 @@ namespace PETScWrappers
const MPI_Comm & com);


/**
* This method associates the PETSc Mat to the instance of the class.
* Infers the number of blocks from A if it is of type MATNEST, otherwise
* the block operator will only have a single block.
*/
void
reinit(Mat A);


/**
* Matrix-vector multiplication: let $dst = M*src$ with $M$ being this
Expand Down Expand Up @@ -292,29 +304,18 @@ namespace PETScWrappers
* are doing.
*
* The PETSc Mat object returned here describes the *entire* matrix,
* not just one of its blocks. Internally, this is done by creating
* a "nested" matrix using Petsc's MatCreateNest object whose individual
* not just one of its blocks. Internally, this is done using
* a "nested" matrix using PETSc's MATNEST object whose individual
* blocks are the blocks of this matrix.
*/
Mat &
petsc_matrix();

/**
* This method assigns the PETSc Mat to the instance of the class.
*
* Note that the matrix is not copied: instead, the instance of this class
* is initialized to use the given matrix. This is useful if you want to
* interpret a PETSc Mat object as a deal.II BlockMatrix, and you already
* have a BlockMatrix object that you want to use for this purpose.
*/
void
assign_petsc_matrix(Mat A);

private:
/**
* A PETSc Mat object that describes the entire block matrix.
* Internally, this is done by creating
* a "nested" matrix using Petsc's MatCreateNest object whose individual
* a "nested" matrix using PETSc's MATNEST object whose individual
* blocks are the blocks of this matrix.
*/
Mat petsc_nest_matrix = nullptr;
Expand All @@ -329,7 +330,7 @@ namespace PETScWrappers
inline BlockSparseMatrix::BlockSparseMatrix(const Mat &A)
: BlockSparseMatrix()
{
this->assign_petsc_matrix(A);
this->reinit(A);
}

inline BlockSparseMatrix &
Expand Down
17 changes: 8 additions & 9 deletions include/deal.II/lac/petsc_block_vector.h
Original file line number Diff line number Diff line change
Expand Up @@ -132,10 +132,12 @@ namespace PETScWrappers

/**
* Create a BlockVector with a PETSc Vec
* It infers the number of blocks from the Vec if it is of type VECNEST,
* otherwise the block vector will only have a single block.
* Internally, we always store a VECNEST vector.
*/
explicit BlockVector(Vec v);


/**
* Destructor. Clears memory
*/
Expand All @@ -155,15 +157,12 @@ namespace PETScWrappers
operator=(const BlockVector &V);

/**
* This method assigns the given PETSc Vec to the instance of the class.
*
* Note that the vector is not copied: instead, the instance of this class
* is initialized to use the given vector. This is useful if you want to
* interpret a PETSc vector as a deal.II vector, and you already have a
* BlockVector that you want to use for this purpose.
* This method associates the PETSc Vec to the instance of the class.
* Infers the number of blocks from v if it is of type VECNEST, otherwise
* the block vector will only have a single block.
*/
void
assign_petsc_vector(Vec v);
reinit(Vec v);

/**
* Reinitialize the BlockVector to contain @p n_blocks of size @p
Expand Down Expand Up @@ -399,7 +398,7 @@ namespace PETScWrappers
: BlockVectorBase<Vector>()
, petsc_nest_vector(nullptr)
{
this->assign_petsc_vector(v);
this->reinit(v);
}

inline BlockVector &
Expand Down
8 changes: 5 additions & 3 deletions include/deal.II/lac/petsc_matrix_base.h
Original file line number Diff line number Diff line change
Expand Up @@ -345,11 +345,13 @@ namespace PETScWrappers
virtual ~MatrixBase() override;

/**
* This method assigns the PETSc Mat to the instance of the class.
*
* This method associates the PETSc Mat to the instance of the class.
* This is particularly useful when performing PETSc to Deal.II operations
* since it allows to reuse the Deal.II MatrixBase and the PETSc Mat
* without incurring in memory copies.
*/
void
assign_petsc_matrix(Mat A);
reinit(Mat A);

/**
* This operator assigns a scalar to a matrix. Since this does usually not
Expand Down
8 changes: 6 additions & 2 deletions include/deal.II/lac/petsc_sparse_matrix.h
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,9 @@ namespace PETScWrappers

/**
* Initialize a Matrix from a PETSc Mat object. Note that we do not copy
* the matrix.
* the matrix. The Mat object is referenced by the newly created instance
* of the class using PetscObjectReference. This is in line with the PETSc
* approach to object ownership, which mimicks std::shared_ptr.
*/
explicit SparseMatrix(const Mat &);

Expand Down Expand Up @@ -396,7 +398,9 @@ namespace PETScWrappers

/**
* Initialize a SparseMatrix from a PETSc Mat object. Note that we do not
* copy the matrix.
* copy the matrix. The Mat object is referenced by the newly created
* instance of the class using PetscObjectReference. This is in line with
* the PETSc approach to object ownership, which mimicks std::shared_ptr.
*/
explicit SparseMatrix(const Mat &);

Expand Down
4 changes: 3 additions & 1 deletion include/deal.II/lac/petsc_vector.h
Original file line number Diff line number Diff line change
Expand Up @@ -289,6 +289,8 @@ namespace PETScWrappers
Vector &
operator=(const dealii::Vector<number> &v);

using VectorBase::reinit;

/**
* Change the dimension of the vector to @p N. It is unspecified how
* resizing the vector affects the memory allocation of this object;
Expand Down Expand Up @@ -476,7 +478,7 @@ namespace PETScWrappers
// in this case a) again takes up a whole lot of memory on the heap,
// and b) is totally dumb since its content would simply be the
// sequence 0,1,2,3,...,n. the best of all worlds would probably be a
// function in Petsc that would take a pointer to an array of
// function in PETSc that would take a pointer to an array of
// PetscScalar values and simply copy n elements verbatim into the
// vector...
for (size_type i = 0; i < v.size(); ++i)
Expand Down
8 changes: 4 additions & 4 deletions include/deal.II/lac/petsc_vector_base.h
Original file line number Diff line number Diff line change
Expand Up @@ -328,13 +328,13 @@ namespace PETScWrappers
operator=(const PetscScalar s);

/**
* This method assigns the PETSc Vec to the instance of the class.
* This is particularly useful when performing PETSc to deal.II operations
* since it allows one to reuse the VectorBase and the PETSc Vec
* This method associates the PETSc Vec to the instance of the class.
* This is particularly useful when performing PETSc to Deal.II operations
* since it allows to reuse the Deal.II VectorBase and the PETSc Vec
* without incurring in memory copies.
*/
void
assign_petsc_vector(Vec v);
reinit(Vec v);

/**
* Test for equality. This function assumes that the present vector and
Expand Down
4 changes: 2 additions & 2 deletions include/deal.II/matrix_free/vector_access_internal.h
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ namespace internal
// below we use type-traits from matrix-free/type_traits.h

// access to generic const vectors that have operator ().
// FIXME: this is wrong for Trilinos/Petsc MPI vectors
// FIXME: this is wrong for Trilinos/PETSc MPI vectors
// where we should first do Partitioner::local_to_global()
template <
typename VectorType,
Expand All @@ -51,7 +51,7 @@ namespace internal


// access to generic non-const vectors that have operator ().
// FIXME: this is wrong for Trilinos/Petsc MPI vectors
// FIXME: this is wrong for Trilinos/PETSc MPI vectors
// where we should first do Partitioner::local_to_global()
template <
typename VectorType,
Expand Down
2 changes: 1 addition & 1 deletion source/lac/petsc_matrix_base.cc
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ namespace PETScWrappers
}

void
MatrixBase::assign_petsc_matrix(Mat A)
MatrixBase::reinit(Mat A)
{
AssertThrow(last_action == ::dealii::VectorOperation::unknown,
ExcMessage("Cannot assign a new Mat."));
Expand Down
2 changes: 1 addition & 1 deletion source/lac/petsc_parallel_block_sparse_matrix.cc
Original file line number Diff line number Diff line change
Expand Up @@ -210,7 +210,7 @@ namespace PETScWrappers
}

void
BlockSparseMatrix::assign_petsc_matrix(Mat A)
BlockSparseMatrix::reinit(Mat A)
{
clear();

Expand Down
4 changes: 2 additions & 2 deletions source/lac/petsc_parallel_block_vector.cc
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ namespace PETScWrappers
}

void
BlockVector::assign_petsc_vector(Vec v)
BlockVector::reinit(Vec v)
{
PetscBool isnest;

Expand Down Expand Up @@ -81,7 +81,7 @@ namespace PETScWrappers
this->components.resize(nb);
for (unsigned int i = 0; i < nb; ++i)
{
this->components[i].assign_petsc_vector(sv[i]);
this->components[i].reinit(sv[i]);
}

this->collect_sizes();
Expand Down
2 changes: 1 addition & 1 deletion source/lac/petsc_vector_base.cc
Original file line number Diff line number Diff line change
Expand Up @@ -181,7 +181,7 @@ namespace PETScWrappers


void
VectorBase::assign_petsc_vector(Vec v)
VectorBase::reinit(Vec v)
{
/* TODO GHOSTED */
AssertThrow(last_action == ::dealii::VectorOperation::unknown,
Expand Down