Merge pull request #58 from aurora-multiphysics/k-collie/grad-div-pro…

…blem

Add grad-div problem and required classes

doc/content/examples/GradDiv.md

@@ -0,0 +1,50 @@
+# Grad-div Problem
+
+## Summary
+
+Solves a diffusion problem for a vector field on a cuboid domain, discretized
+using $H(\mathrm{div})$ conforming Raviart-Thomas elements. This example is
+based on [MFEM Example 4](https://mfem.org/examples/) and is relevant for
+solving Maxwell's equations using potentials without the Coulomb gauge.
+
+
+## Description
+
+This problem solves a grad-div equation with strong form:
+
+\begin{equation}
+\begin{split}
+\vec\nabla \left( \alpha \vec\nabla \cdot \vec F \right) + \beta \vec F = \vec f \,\,\,&\mathrm{on}\,\, \Omega \\
+\vec F \cdot \hat n= \vec g \,\,\, &\mathrm{on}\,\, \partial\Omega
+\end{split}
+\end{equation}
+
+where
+
+\begin{equation}
+\vec g = \begin{pmatrix}
+    \cos(k x)\sin(k y)\\
+    \cos(k y)\sin(k x)\\
+    0
+\end{pmatrix},\,\,\,
+\vec f = \left( \beta + 2 \alpha k^2 \right) \vec g
+\end{equation}
+
+In this example, we solve this using the weak form
+
+!equation
+(\alpha \vec\nabla \cdot \vec F , \vec\nabla \cdot \vec v)_\Omega + (\beta \vec F, \vec F)_\Omega
+= (\vec f, \vec v)_\Omega \,\,\, \forall \vec v \in V
+
+where
+
+\begin{equation}
+\begin{split}
+\vec F \in H(\mathrm{div})(\Omega) &: \vec F \cdot \hat n= \vec g \,\,\, \mathrm{on}\,\, \partial\Omega \\
+\vec v \in H(\mathrm{div})(\Omega) &: \vec v \cdot \hat n= \vec 0 \,\,\, \mathrm{on}\,\, \partial\Omega
+\end{split}
+\end{equation}
+
+## Example File
+
+!listing kernels/graddiv.i

doc/content/examples/index.md

@@ -25,3 +25,8 @@ starting point for users to adapt:
 - [Definite Maxwell](examples/DefiniteMaxwell.md): Solves a 3D electromagnetic diffusion problem for
   the electric field on a cube missing an octant, discretized using $H(\mathrm{curl})$ conforming
   Nédélec elements. This example is based on [MFEM Example 3](https://mfem.org/examples/).
+
+- [GradDiv](examples/GradDiv.md): Solves a diffusion problem for a vector field
+  on a cuboid domain, discretized using $H(\mathrm{div})$ conforming
+  Raviart-Thomas elements. This example is based on [MFEM Example 4](https://mfem.org/examples/)
+  and is relevant for solving Maxwell's equations using potentials without the Coulomb gauge.

doc/content/source/bcs/MFEMVectorNormalDirichletBC.md

@@ -0,0 +1,20 @@
+# MFEMVectorNormalDirichletBC
+
+## Summary
+
+!syntax description /BCs/MFEMVectorNormalDirichletBC
+
+## Overview
+
+Boundary condition for enforcing an essential (Dirichlet) boundary condition on the normal
+components of a $H(\mathrm{div})$ conforming vector FE at a boundary.
+
+## Example Input File Syntax
+
+!listing test/tests/kernels/graddiv.i block=BCs
+
+!syntax parameters /BCs/MFEMVectorNormalDirichletBC
+
+!syntax inputs /BCs/MFEMVectorNormalDirichletBC
+
+!syntax children /BCs/MFEMVectorNormalDirichletBC

doc/content/source/kernels/MFEMDivDivKernel.md

@@ -0,0 +1,29 @@
+# MFEMDivDivKernel
+
+## Summary
+
+!syntax description /Kernels/MFEMDivDivKernel
+
+## Overview
+
+Adds the domain integrator for integrating the bilinear form
+
+!equation
+(k\vec\nabla \cdot \vec u, \vec\nabla \cdot \vec v)_\Omega \,\,\, \forall \vec v \in V
+
+where $\vec u, \vec v \in H(\mathrm{div})$ and $k$ is a scalar coefficient.
+
+This term arises from the weak form of the grad div operator
+
+!equation
+-\vec\nabla \left(k \vec\nabla \cdot \vec u\right)
+
+## Example Input File Syntax
+
+!listing graddiv.i
+
+!syntax parameters /Kernels/MFEMDivDivKernel
+
+!syntax inputs /Kernels/MFEMDivDivKernel
+
+!syntax children /Kernels/MFEMDivDivKernel

doc/content/source/solvers/MFEMHypreADS.md

@@ -0,0 +1,22 @@
+# MFEMHypreADS
+
+## Summary
+
+!syntax description /Solver/MFEMHypreADS
+
+## Overview
+
+Defines and builds an `mfem::HypreADS` solver to use as a preconditioner or solver to solve the
+MFEM equation system. Most effective for preconditioning and solving a grad-divergence problem when using
+Raviart-Thomas elements, in which case the $H(\mathrm{div})$ FE space should be passed to the
+`mfem::HypreADS` solver during construction.
+
+## Example Input File Syntax
+
+!listing test/tests/kernels/graddiv.i block=FESpace Preconditioner Solver
+
+!syntax parameters /Solver/MFEMHypreADS
+
+!syntax inputs /Solver/MFEMHypreADS
+
+!syntax children /Solver/MFEMHypreADS

include/bcs/MFEMVectorNormalDirichletBC.h

@@ -0,0 +1,12 @@
+#pragma once
+
+#include "MFEMVectorDirichletBCBase.h"
+
+class MFEMVectorNormalDirichletBC : public MFEMVectorDirichletBCBase
+{
+public:
+  static InputParameters validParams();
+  MFEMVectorNormalDirichletBC(const InputParameters & parameters);
+  ~MFEMVectorNormalDirichletBC() override = default;
+  void ApplyBC(mfem::GridFunction & gridfunc, mfem::Mesh * mesh_) override;
+};

include/kernels/MFEMDivDivKernel.h

@@ -0,0 +1,20 @@
+#pragma once
+#include "MFEMKernel.h"
+
+/*
+(α∇.u, ∇.u')
+*/
+class MFEMDivDivKernel : public MFEMKernel<mfem::BilinearFormIntegrator>
+{
+public:
+  static InputParameters validParams();
+
+  MFEMDivDivKernel(const InputParameters & parameters);
+  ~MFEMDivDivKernel() override {}
+
+  virtual mfem::BilinearFormIntegrator * createIntegrator() override;
+
+protected:
+  std::string _coef_name;
+  mfem::Coefficient & _coef;
+};

include/solvers/MFEMHypreADS.h

@@ -0,0 +1,24 @@
+#pragma once
+#include "MFEMSolverBase.h"
+#include "MFEMFESpace.h"
+
+/**
+ * Wrapper for mfem::HypreADS solver.
+ */
+class MFEMHypreADS : public MFEMSolverBase
+{
+public:
+  static InputParameters validParams();
+
+  MFEMHypreADS(const InputParameters &);
+
+  /// Returns a shared pointer to the instance of the Solver derived-class.
+  std::shared_ptr<mfem::Solver> getSolver() const override { return _preconditioner; }
+
+protected:
+  void constructSolver(const InputParameters & parameters) override;
+
+private:
+  const MFEMFESpace & _mfem_fespace;
+  std::shared_ptr<mfem::HypreADS> _preconditioner{nullptr};
+};

src/bcs/MFEMVectorNormalDirichletBC.C

@@ -0,0 +1,25 @@
+#include "MFEMVectorNormalDirichletBC.h"
+
+registerMooseObject("PlatypusApp", MFEMVectorNormalDirichletBC);
+
+InputParameters
+MFEMVectorNormalDirichletBC::validParams()
+{
+  InputParameters params = MFEMVectorDirichletBCBase::validParams();
+  params.addClassDescription(
+      "Applies a Dirichlet condition to the normal components of a vector variable.");
+  return params;
+}
+
+MFEMVectorNormalDirichletBC::MFEMVectorNormalDirichletBC(const InputParameters & parameters)
+  : MFEMVectorDirichletBCBase(parameters)
+{
+}
+
+void
+MFEMVectorNormalDirichletBC::ApplyBC(mfem::GridFunction & gridfunc, mfem::Mesh * mesh_)
+{
+  mfem::Array<int> ess_bdrs(mesh_->bdr_attributes.Max());
+  ess_bdrs = GetMarkers(*mesh_);
+  gridfunc.ProjectBdrCoefficientNormal(*_vec_coef->getVectorCoefficient(), ess_bdrs);
+}

src/kernels/MFEMDivDivKernel.C

@@ -0,0 +1,34 @@
+#include "MFEMDivDivKernel.h"
+#include "MFEMProblem.h"
+
+registerMooseObject("PlatypusApp", MFEMDivDivKernel);
+
+InputParameters
+MFEMDivDivKernel::validParams()
+{
+  InputParameters params = MFEMKernel::validParams();
+  params.addClassDescription(
+      "Adds the domain integrator to an MFEM problem for the bilinear form "
+      "$(k\\vec\\nabla \\cdot \\vec u, \\vec\\nabla \\cdot \\vec v)_\\Omega$ "
+      "arising from the weak form of the grad-div operator "
+      "$-\\vec\\nabla \\left( k \\vec\\nabla \\cdot \\vec u \\right)$.");
+
+  params.addParam<std::string>("coefficient", "Name of property k to multiply the Laplacian by");
+
+  return params;
+}
+
+MFEMDivDivKernel::MFEMDivDivKernel(const InputParameters & parameters)
+  : MFEMKernel(parameters),
+    _coef_name(getParam<std::string>("coefficient")),
+    // FIXME: The MFEM bilinear form can also handle vector and matrix
+    // coefficients, so ideally we'd handle all three too.
+    _coef(getMFEMProblem().getProperties().getScalarProperty(_coef_name))
+{
+}
+
+mfem::BilinearFormIntegrator *
+MFEMDivDivKernel::createIntegrator()
+{
+  return new mfem::DivDivIntegrator(_coef);
+}

src/solvers/MFEMHypreADS.C

@@ -0,0 +1,28 @@
+#pragma once
+#include "MFEMHypreADS.h"
+
+registerMooseObject("PlatypusApp", MFEMHypreADS);
+
+InputParameters
+MFEMHypreADS::validParams()
+{
+  InputParameters params = MFEMSolverBase::validParams();
+  params.addClassDescription("Hypre auxiliary-space divergence solver and preconditioner for the "
+                             "iterative solution of MFEM equation systems.");
+  params.addParam<UserObjectName>("fespace", "H(div) FESpace to use in HypreADS setup.");
+  params.addParam<int>("print_level", 2, "Set the solver verbosity.");
+  return params;
+}
+
+MFEMHypreADS::MFEMHypreADS(const InputParameters & parameters)
+  : MFEMSolverBase(parameters), _mfem_fespace(getUserObject<MFEMFESpace>("fespace"))
+{
+  constructSolver(parameters);
+}
+
+void
+MFEMHypreADS::constructSolver(const InputParameters & parameters)
+{
+  _preconditioner = std::make_shared<mfem::HypreADS>(_mfem_fespace.getFESpace().get());
+  _preconditioner->SetPrintLevel(getParam<int>("print_level"));
+}

test/tests/kernels/gold/OutputData/GradDiv/Run0/Cycle000000/data.pvtu

@@ -0,0 +1,20 @@
+<?xml version="1.0"?>
+<VTKFile type="PUnstructuredGrid" version ="0.1" byte_order="LittleEndian">
+<PUnstructuredGrid GhostLevel="0">
+<PPoints>
+	<PDataArray type="Float64"  Name="Points" NumberOfComponents="3" format="ascii"/>
+</PPoints>
+<PCells>
+	<PDataArray type="Int32"  Name="connectivity" NumberOfComponents="1" format="ascii"/>
+	<PDataArray type="Int32"  Name="offsets"      NumberOfComponents="1" format="ascii"/>
+	<PDataArray type="UInt8"  Name="types"        NumberOfComponents="1" format="ascii"/>
+</PCells>
+<PPointData>
+<PDataArray type="Float64" Name="F" NumberOfComponents="3" ComponentName0="0" ComponentName1="1" ComponentName2="2" format="ascii" />
+</PPointData>
+<PCellData>
+	<PDataArray type="Int32" Name="attribute" NumberOfComponents="1" format="ascii"/>
+</PCellData>
+<Piece Source="proc000000.vtu"/>
+</PUnstructuredGrid>
+</VTKFile>