17 #include "../config/config.hpp"
27 #include "../general/mem_manager.hpp"
29 #include "petscconf.h"
30 #if !defined(PETSC_USE_REAL_DOUBLE)
31 #error "MFEM does not work with PETSc compiled without double precision"
33 #if defined(PETSC_USE_COMPLEX)
34 #error "MFEM does not work with PETSc compiled with complex numbers support"
36 #if defined(PETSC_USE_64BIT_INDICES) && !defined(HYPRE_BIGINT)
37 #error "Mismatch between HYPRE (32bit) and PETSc (64bit) integer types"
39 #if !defined(PETSC_USE_64BIT_INDICES) && defined(HYPRE_BIGINT)
40 #error "Mismatch between HYPRE (64bit) and PETSc (32bit) integer types"
43 #include "petscversion.h"
44 #if PETSC_VERSION_GE(3,12,0)
45 #include "petscsystypes.h"
69 typedef struct ::_p_Vec *
Vec;
70 typedef struct ::_p_Mat *
Mat;
71 typedef struct ::_p_KSP *
KSP;
72 typedef struct ::_p_PC *
PC;
73 typedef struct ::_p_SNES *
SNES;
74 typedef struct ::_p_TS *
TS;
100 int size_,
bool usedev_)
110 bool read_,
bool write_,
bool usedev_)
150 class ParFiniteElementSpace;
151 class PetscParMatrix;
205 bool transpose =
false,
bool allocate =
true);
211 bool allocate =
true);
301 void Print(
const char *fname = NULL,
bool binary =
false)
const;
303 const double *
Read(
bool=
true)
const override;
304 const double *
HostRead()
const override;
305 double *
Write(
bool=
true)
override;
349 void BlockDiagonalConstructor(MPI_Comm comm,
PetscInt *row_starts,
353 void SetUpForDevice();
445 {
Mult(1.0, x, 0.0, y); }
502 void Print(
const char *fname = NULL,
bool binary =
false)
const;
556 PetscParMatrix *
ParMult(
const PetscParMatrix *A,
const PetscParMatrix *B);
559 PetscParMatrix *
RAP(PetscParMatrix *Rt, PetscParMatrix *A, PetscParMatrix *P);
566 PetscParMatrix *
RAP(PetscParMatrix *A, PetscParMatrix *P);
569 PetscParMatrix *
RAP(HypreParMatrix *A, PetscParMatrix *P);
576 void EliminateBC(PetscParMatrix &A, PetscParMatrix &Ae,
577 const Array<int> &ess_dof_list,
const Vector &X, Vector &B);
591 bctype(type_), setup(false), eval_t(0.0),
592 eval_t_cached(std::numeric_limits<double>::min()) {}
607 {
mfem_error(
"PetscBCHandler::Eval method not overloaded"); }
641 double eval_t_cached;
655 const char*
GetName() {
return name.c_str(); }
661 class PetscSolverMonitor;
709 void Customize(
bool customize =
true)
const;
747 bool wrap =
true,
bool iter_mode =
false);
749 const std::string &prefix = std::string(),
bool iter_mode =
false);
756 const std::string &prefix = std::string(),
bool iter_mode =
false);
783 PetscPCGSolver(MPI_Comm comm,
const std::string &prefix = std::string(),
784 bool iter_mode =
false);
786 bool iter_mode =
false);
788 const std::string &prefix = std::string(),
bool iter_mode =
false);
800 const std::string &prefix = std::string());
802 const std::string &prefix = std::string());
804 const std::string &prefix = std::string());
868 const std::string &prefix = std::string());
871 const std::string &prefix = std::string());
879 const std::string &prefix = std::string());
890 const std::string &prefix = std::string());
899 const std::string &prefix = std::string());
901 const std::string &prefix = std::string());
922 Vector &W,
bool &changed_y,
bool &changed_w));
949 PetscODESolver(MPI_Comm comm,
const std::string &prefix = std::string());
964 virtual void Step(
Vector &x,
double &
t,
double &dt);
965 virtual void Run(
Vector &x,
double &
t,
double &dt,
double t_final);
985 MFEM_ABORT(
"MonitorSolution() is not implemented!")
991 MFEM_ABORT(
"MonitorResidual() is not implemented!")
1001 #endif // MFEM_USE_MPI
1002 #endif // MFEM_USE_PETSC
void SetPreconditioner(Solver &precond)
void SetPrintLevel(int plev)
PetscParVector * B
Right-hand side and solution vector.
void CreatePrivateContext()
void Print(const char *fname=NULL, bool binary=false) const
Prints the matrix (to stdout if fname is NULL)
void trans(const Vector &u, Vector &x)
void EliminateBC(const HypreParMatrix &A, const HypreParMatrix &Ae, const Array< int > &ess_dof_list, const Vector &X, Vector &B)
Eliminate essential BC specified by ess_dof_list from the solution X to the r.h.s. B.
void SetUpdate(void(*update)(Operator *op, int it, const Vector &F, const Vector &X, const Vector &D, const Vector &P))
virtual void SetOperator(const Operator &op)
Specification of the nonlinear operator.
Abstract class for PETSc's preconditioners.
PetscParMatrix & operator+=(const PetscParMatrix &B)
petsc::Vec x
The actual PETSc object.
virtual void SetOperator(const Operator &op)
bool clcustom
Boolean to handle SetFromOptions calls.
void SetObjective(void(*obj)(Operator *op, const Vector &x, double *f))
Specification of an objective function to be used for line search.
virtual void MultTranspose(const Vector &b, Vector &x) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
void MakeAliasForSync(Memory< double > &base_, int offset_, int size_, bool read_, bool write_, bool usedev_)
petsc::Mat A
The actual PETSc object.
Wrapper for PETSc's matrix class.
PetscParVector & operator+=(const PetscParVector &y)
const Array< int > * nat_dof
void ApplyBC(const Vector &x, Vector &y)
y = x on ess_tdof_list_c and y = g (internally evaluated) on ess_tdof_list
void SetHostValid() const
Abstract class for PETSc's linear solvers.
Abstract class for PETSc's solvers.
PetscODESolver::Type GetType() const
virtual void Run(Vector &x, double &t, double &dt, double t_final)
Perform time integration from time t [in] to time tf [in].
void SetPostCheck(void(*post)(Operator *op, const Vector &X, Vector &Y, Vector &W, bool &changed_y, bool &changed_w))
PetscBCHandler(Type type_=ZERO)
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
ParFiniteElementSpace * fespace
Type GetType() const
Returns the type of boundary conditions.
Base abstract class for first order time dependent operators.
PetscInt GetNumCols() const
Returns the local number of columns.
void MakeAliasForSync(const Memory< double > &base_, int offset_, int size_, bool usedev_)
void SetType(PetscODESolver::Type)
void ConvertOperator(MPI_Comm comm, const Operator &op, petsc::Mat *B, Operator::Type tid)
virtual void MonitorSolver(PetscSolver *solver)
Generic monitor to take access to the solver.
Pointer to an Operator of a specified type.
void Mult(double a, const Vector &x, double b, Vector &y) const
Matvec: y = a A x + b y.
PetscBDDCSolver(MPI_Comm comm, Operator &op, const PetscBDDCSolverParams &opts=PetscBDDCSolverParams(), const std::string &prefix=std::string())
PetscParVector & AddValues(const Array< PetscInt > &, const Array< PetscScalar > &)
Add values in a vector.
virtual ~PetscSolver()
Destroy the PetscParVectors allocated (if any).
PetscNonlinearSolver(MPI_Comm comm, const std::string &prefix=std::string())
bool DeviceRequested() const
PetscInt GetGlobalNumRows() const
Returns the global number of rows.
void SetBlockSize(PetscInt rbs, PetscInt cbs=-1)
Set row and column block sizes of a matrix.
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
void Init()
Initialize with defaults. Does not initialize inherited members.
PetscPreconditioner(MPI_Comm comm, const std::string &prefix=std::string())
PetscClassId cid
The class id of the actual PETSc object.
Abstract parallel finite element space.
void SetRelTol(double tol)
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Wrapper for syncing PETSc's vector memory.
void Randomize(PetscInt seed=0)
Set random values.
void Destroy()
Delete all owned data. Does not perform re-initialization with defaults.
void SetMonitor(PetscSolverMonitor *ctx)
Sets user-defined monitoring routine.
void SetMaxIter(int max_iter)
Abstract class for PETSc's ODE solvers.
PetscParVector & operator*=(PetscScalar d)
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B, bool own_matrix)
void EliminateRowsCols(const Array< int > &rows_cols, const PetscParVector &X, PetscParVector &B, double diag=1.)
Eliminate rows and columns from the matrix, and rows from the vector B. Modify B with the BC values i...
const double * HostRead() const override
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
void Customize(bool customize=true) const
Customize object with options set.
PetscInt M() const
Returns the global number of rows.
unsigned flags
Bit flags defined from the FlagMask enum.
virtual void Mult(const Vector &b, Vector &x) const
Application of the solver.
virtual ~PetscLinearSolver()
PetscInt GetNumRows() const
Returns the local number of rows.
void SetHostInvalid() const
void Print(const char *fname=NULL, bool binary=false) const
Prints the vector (to stdout if fname is NULL)
void SetDeviceInvalid() const
void ResetArray()
Reset the PETSc Vec object to use its default data. Call this method after the use of PlaceArray()...
const double * GetHostPointer() const
PetscParVector & SetValues(const Array< PetscInt > &, const Array< PetscScalar > &)
Set values in a vector.
double * Write(bool=true) override
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
virtual void Init(TimeDependentOperator &f_, enum PetscODESolver::Type type)
Initialize the ODE solver.
PetscParMatrix * Transpose(bool action=false)
Returns the transpose of the PetscParMatrix.
double f(const Vector &xvec)
PetscObject obj
The actual PETSc object (KSP, PC, SNES or TS).
Abstract class for monitoring PETSc's solvers.
void SetDeviceValid() const
PetscLinearSolver(MPI_Comm comm, const std::string &prefix=std::string(), bool wrap=true, bool iter_mode=false)
struct s_NavierContext ctx
void SetJacobianType(Operator::Type type)
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Auxiliary class for BDDC customization.
virtual ~PetscPreconditioner()
MPI_Comm GetComm() const
Get the associated MPI communicator.
void FixResidualBC(const Vector &x, Vector &y)
y = x-g on ess_tdof_list, the rest of y is unchanged
bool UseDevice() const override
Return the device flag of the Memory object used by the Vector.
void SetFlagsFromMask_() const
void SetType(enum Type type_)
Sets the type of boundary conditions.
PetscParMatrix & operator=(const PetscParMatrix &B)
bool ReadRequested() const
PetscInt GetColStart() const
Returns the global index of the first local column.
virtual ~PetscParVector()
Calls PETSc's destroy function.
PetscFieldSplitSolver(MPI_Comm comm, Operator &op, const std::string &prefix=std::string())
PetscBCHandler * bchandler
Handler for boundary conditions.
PetscParMatrix & operator-=(const PetscParMatrix &B)
virtual Solver * NewPreconditioner(const OperatorHandle &oh)=0
double * ReadWrite(bool=true) override
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
void ScaleRows(const Vector &s)
Scale the local row i by s(i).
void Zero(Vector &x)
Replace boundary dofs with 0.
MPI_Comm GetComm() const
Get the associated MPI communicator.
bool WriteRequested() const
void Shift(double s)
Shift diagonal by a constant.
double * HostReadWrite() override
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
void UpdateVecFromFlags()
Update PETSc's Vec after having accessed its data via GetMemory()
void SetTime(double t)
Sets the current time.
ID for class PetscParMatrix, MATSHELL format.
void Reset()
Reset the memory to be empty, ensuring that Delete() will be a no-op.
void SetMat(petsc::Mat newA)
Replace the inner Mat Object. The reference count of newA is increased.
PetscParMatrix()
Create an empty matrix to be used as a reference to an existing matrix.
void MakeWrapper(MPI_Comm comm, const Operator *op, petsc::Mat *B)
Creates a wrapper around a mfem::Operator op using PETSc's MATSHELL object and returns the Mat in B...
Abstract class for PETSc's nonlinear solvers.
void SetBCHandler(PetscBCHandler *bch)
Sets the object to handle essential boundary conditions.
Wrapper for hypre's parallel vector class.
void SetBlockSize(PetscInt bs)
Set block size of a vector.
const double * Read(bool=true) const override
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
void ZeroBC(const Vector &x, Vector &y)
y = x on ess_tdof_list_c and y = 0 on ess_tdof_list
virtual void Eval(double t, Vector &g)
Boundary conditions evaluation.
Type
Enumeration defining IDs for some classes derived from Operator.
PetscInt NNZ() const
Returns the number of nonzeros.
PetscPreconditionerFactory(const std::string &name_="MFEM Factory")
virtual void Mult(const Vector &b, Vector &x) const
Application of the preconditioner.
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
void Sync(const Memory &other) const
Copy the host/device pointer validity flags from other to *this.
PetscInt N() const
Returns the global number of columns.
PetscParVector * GetY() const
Returns the inner vector in the range of A (it creates it if needed)
void SetJacobianType(Operator::Type type)
double * HostWrite() override
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
bool operatorset
Boolean to handle SetOperator calls.
PetscParVector * X
Auxiliary vectors for typecasting.
void MakeAlias(const Memory &base, int offset, int size)
Create a memory object that points inside the memory object base.
void MultTranspose(double a, const Vector &x, double b, Vector &y) const
Matvec transpose: y = a A^T x + b y.
bool IsAliasForSync() const
PetscParVector & operator-=(const PetscParVector &y)
virtual ~PetscPreconditionerFactory()
virtual ~PetscSolverMonitor()
virtual ~PetscNonlinearSolver()
void SetComputeNetFlux(bool net=true)
Setup BDDC with no-net-flux local solvers. Needs a ParFiniteElementSpace attached.
virtual ~PetscODESolver()
virtual ~PetscBCHandler()
PetscSolver()
Construct an empty PetscSolver. Initialize protected objects to NULL.
void SetPreconditionerFactory(PetscPreconditionerFactory *factory)
Sets the object for the creation of the preconditioner.
void MFEMInitializePetsc()
Convenience functions to initialize/finalize PETSc.
petsc::Mat ReleaseMat(bool dereference)
Release the PETSc Mat object. If dereference is true, decrement the refcount of the Mat object...
PetscParVector(MPI_Comm comm, PetscInt glob_size, PetscInt *col=NULL)
Creates vector with given global size and partitioning of the columns.
virtual void MonitorSolution(PetscInt it, PetscReal norm, const Vector &x)
Monitor the solution vector x.
ID for class PetscParMatrix, MATAIJ format.
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
void SetSpace(ParFiniteElementSpace *fe)
MPI_Comm GetComm() const
Get the associated MPI communicator.
PetscInt GlobalSize() const
Returns the global number of rows.
const Array< int > * ess_dof
PetscInt GetRowStart() const
Returns the global index of the first local row.
struct _p_PetscObject * PetscObject
PetscH2Solver(Operator &op, ParFiniteElementSpace *fes, const std::string &prefix=std::string())
void ResetMemory()
Completes the operation started with PlaceMemory.
void SetAbsTol(double tol)
Helper class for handling essential boundary conditions.
Class used by MFEM to store pointers to host and/or device memory.
PetscParVector * GetX() const
Returns the inner vector in the domain of A (it creates it if needed)
Array< int > & GetTDofs()
Gets essential dofs (local, per-process numbering)
PetscParVector & operator=(PetscScalar d)
Set constant values.
PetscSolverMonitor(bool monitor_sol=false, bool monitor_res=true)
virtual void MonitorResidual(PetscInt it, PetscReal norm, const Vector &r)
Monitor the residual vector r.
PetscParMatrix * TripleMatrixProduct(PetscParMatrix *R, PetscParMatrix *A, PetscParMatrix *P)
Returns the matrix R * A * P.
void PlaceArray(PetscScalar *temp_data)
Temporarily replace the data of the PETSc Vec object. To return to the original data array...
void FreePrivateContext()
Vector * GlobalVector() const
Returns the global vector in each processor.
PetscODESolver(MPI_Comm comm, const std::string &prefix=std::string())
void EliminateRows(const Array< int > &rows)
Eliminate only the rows from the matrix.
PetscInt GetGlobalNumCols() const
Returns the global number of columns.
PetscPCGSolver(MPI_Comm comm, const std::string &prefix=std::string(), bool iter_mode=false)
void SetNatBdrDofs(const Array< int > *natdofs, bool loc=false)
Specify dofs on the natural boundary.
Wrapper for hypre's ParCSR matrix class.
virtual void Mult(const Vector &b, Vector &x) const
Application of the solver.
void operator*=(double s)
Scale all entries by s: A_scaled = s*A.
void SetEssBdrDofs(const Array< int > *essdofs, bool loc=false)
Specify dofs on the essential boundary.
void SetTDofs(Array< int > &list)
Sets essential dofs (local, per-process numbering)
virtual void MultTranspose(const Vector &b, Vector &x) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
void PlaceMemory(Memory< double > &, bool=false)
This requests write access from where the memory is valid and temporarily replaces the corresponding ...
void * private_ctx
Private context for solver.
void SetUp(PetscInt n)
SetUp the helper object, where n is the size of the solution vector.
const double * GetDevicePointer() const
void MakeRef(const PetscParMatrix &master)
Makes this object a reference to another PetscParMatrix.
void ScaleCols(const Vector &s)
Scale the local col i by s(i).
virtual ~PetscParMatrix()
Calls PETSc's destroy function.