31 Solver(curlcurl_op_.Height()),
32 curlcurl_op(curlcurl_op_),
35 pispacesolver(pispacesolver_),
36 gspacesolver(gspacesolver_),
38 ess_tdof_list(ess_tdof_list_)
46void GeneralAMS::FormResidual(
const Vector& rhs,
const Vector& x,
49 curlcurl_op.
Mult(x, residual);
72 MFEM_ASSERT(x.
Size() == y.
Size(),
"Sizes don't match!");
73 MFEM_ASSERT(curlcurl_op.
Height() == x.
Size(),
"Sizes don't match!");
83 FormResidual(x, y, residual);
87 gspacecorrection = 0.0;
88 gspacesolver.
Mult(gspacetemp, gspacecorrection);
89 gradient.
Mult(gspacecorrection, residual);
93 FormResidual(x, y, residual);
99 pispacecorrection = 0.0;
100 pispacesolver.
Mult(pispacetemp, pispacecorrection);
101 pi.
Mult(pispacecorrection, residual);
105 FormResidual(x, y, residual);
107 gspacecorrection = 0.0;
108 gspacesolver.
Mult(gspacetemp, gspacecorrection);
109 gradient.
Mult(gspacecorrection, residual);
113 FormResidual(x, y, residual);
115 smoother.
Mult(residual, temp);
129 comm(mesh_lor.GetComm()),
135 inner_aux_iterations(0)
152 if (alpha_coeff == NULL)
161 if (beta_mcoeff != NULL)
163 MFEM_VERIFY(beta_coeff == NULL,
"Only one beta coefficient should be defined.");
166 else if (beta_coeff != NULL)
184 if (cg_iterations > 0)
186 SetupCG(curlcurl_oper, pi, cg_iterations);
212 Solver(curlcurl_oper.Height()),
213 comm(mesh_lor.GetComm()),
219 inner_aux_iterations(0)
234 if (beta_mcoeff != NULL)
236 MFEM_VERIFY(beta_coeff == NULL,
"Only one beta coefficient should be defined.");
239 else if (beta_coeff != NULL)
266 if (cg_iterations > 0)
268 SetupCG(curlcurl_oper, g, cg_iterations);
278void MatrixFreeAuxiliarySpace::SetupCG(
280 int inner_cg_iterations)
282 MFEM_ASSERT(conn.
Height() == curlcurl_oper.
Width(),
283 "Operators don't match!");
284 matfree =
new RAPOperator(conn, curlcurl_oper, conn);
286 "Operators don't match!");
291 if (inner_cg_iterations > 99)
306void MatrixFreeAuxiliarySpace::SetupVCycle()
308 aspacewrapper = lor_pc;
311class ZeroWrapAMG :
public Solver
315 ZeroWrapAMG(HypreParMatrix& mat, Array<int>& ess_tdof_list_,
316 const bool useAmgX) :
318 ZeroWrapAMG(HypreParMatrix& mat, Array<int>& ess_tdof_list_) :
325 const bool amgx_verbose =
false;
326 AmgXSolver *amgx =
new AmgXSolver(mat.GetComm(),
329 amgx->SetOperator(mat);
335 HypreBoomerAMG *amg =
new HypreBoomerAMG(mat);
336 amg->SetPrintLevel(0);
341 void Mult(
const Vector& x, Vector& y)
const override
345 auto Y = y.HostReadWrite();
346 for (
int k : ess_tdof_list)
352 void SetOperator(
const Operator&)
override { }
354 ~ZeroWrapAMG()
override
361 Array<int>& ess_tdof_list;
364void MatrixFreeAuxiliarySpace::SetupAMG(
int system_dimension)
366 if (system_dimension == 0)
370 lor_pc =
new ZeroWrapAMG(*lor_matrix, ess_tdof_list, useAmgX);
372 lor_pc =
new ZeroWrapAMG(*lor_matrix, ess_tdof_list);
381 const bool amgx_verbose =
false;
382 AmgXSolver *amgx =
new AmgXSolver(lor_matrix->
GetComm(),
385 amgx->SetOperator(*lor_matrix);
391 HypreBoomerAMG* hpc =
new HypreBoomerAMG(*lor_matrix);
392 hpc->SetSystemsOptions(system_dimension);
393 hpc->SetPrintLevel(0);
402 MPI_Comm_rank(comm, &rank);
405 aspacewrapper->
Mult(x, y);
409 inner_aux_iterations += q;
418 if (lor_pc != aspacewrapper) {
delete aspacewrapper; }
419 if (cg != aspacewrapper) {
delete cg; }
433 int inner_pi_iterations,
int inner_g_iterations,
Solver * nd_smoother) :
445 smoother = nd_smoother;
449 const double scale = 0.25;
477 beta_mcoeff, ess_bdr, oper,
486 beta_coeff, beta_mcoeff,
491 inner_pi_iterations);
493 general_ams =
new GeneralAMS(oper, *Pi, *Gradient, *Pispacesolver,
494 *Gspacesolver, *smoother, ess_tdof_list);
505 delete Pispacesolver;
int Size() const
Return the logical size of the array.
@ GaussLobatto
Closed type.
Conjugate gradient method.
void SetOperator(const Operator &op) override
Set/update the solver for the given operator.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
void AddDomainInterpolator(DiscreteInterpolator *di)
Adds a domain interpolator. Assumes ownership of di.
virtual void Assemble(int skip_zeros=1)
Construct the internal matrix representation of the discrete linear operator.
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::FULL.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Mesh * GetMesh() const
Returns the mesh.
const FiniteElement * GetTypicalFE() const
Return GetFE(0) if the local mesh is not empty; otherwise return a typical FE based on the Geometry t...
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Perform AMS cycle with generic Operator objects.
void Mult(const Vector &x, Vector &y) const override
Operator application: y=A(x).
GeneralAMS(const Operator &curlcurl_op_, const Operator &pi_, const Operator &gradient_, const Operator &pispacesolver_, const Operator &gspacesolver_, const Operator &smoother_, const Array< int > &ess_tdof_list_)
Constructor.
Arbitrary order H1-conforming (continuous) finite elements.
MPI_Comm GetComm() const
MPI communicator.
void SetRelTol(real_t rtol)
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
void SetMaxIter(int max_it)
Base class for Matrix Coefficients that optionally depend on time and space.
MatrixFreeAMS(ParBilinearForm &aform, Operator &oper, ParFiniteElementSpace &nd_fespace, Coefficient *alpha_coeff, Coefficient *beta_coeff, MatrixCoefficient *beta_mcoeff, Array< int > &ess_bdr, #ifdef MFEM_USE_AMGX bool useAmgX=false, #endif int inner_pi_its=0, int inner_g_its=1, Solver *nd_smoother=NULL)
Construct matrix-free AMS preconditioner.
Auxiliary space solvers for MatrixFreeAMS preconditioner.
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
MatrixFreeAuxiliarySpace(ParMesh &mesh_lor, Coefficient *alpha_coeff, Coefficient *beta_coeff, MatrixCoefficient *beta_mcoeff, Array< int > &ess_bdr, Operator &curlcurl_oper, Operator &pi, #ifdef MFEM_USE_AMGX bool useAmgX_, #endif int cg_iterations=0)
Pi space constructor.
~MatrixFreeAuxiliarySpace()
int Dimension() const
Dimension of the reference space used within the elements.
Jacobi smoothing for a given bilinear form (no matrix necessary).
Operator(int s=0)
Construct a square Operator with given size s (default 0).
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
@ DIAG_ONE
Set the diagonal value to one.
@ DIAG_KEEP
Keep the diagonal value.
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
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 ...
virtual void FormRectangularSystemMatrix(OperatorHandle &A)
Return in A a parallel (on truedofs) version of this operator.
Abstract parallel finite element space.
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
ParMesh * GetParMesh() const
Class for parallel meshes.
static ParMesh MakeRefined(ParMesh &orig_mesh, int ref_factor, int ref_type)
Create a uniformly refined (by any factor) version of orig_mesh.
The operator x -> R*A*P*x constructed through the actions of R^T, A and P.
Solver(int s=0, bool iter_mode=false)
Initialize a square Solver with size s.
int Size() const
Returns the size of the vector.