MFEM Non-linear support

framework/doc/content/source/mfem/functions/MFEMGradientGridFunction.md

@@ -0,0 +1,16 @@
+# MFEMGradientGridFunction
+
+!if! function=hasCapability('mfem')
+
+## Summary
+
+This is class which can be use as vector cofficients to MFEM problems.
+
+## Overview
+
+`MFEMGradientGridFunction` is intended to allow the specification of `mfem::VectorCoefficient` object to add to the MFEM problem in a manner consistent with the standard MOOSE Materials system.
+
+!if-end!
+
+!else
+!include mfem/mfem_warning.md

framework/doc/content/source/mfem/kernels/MFEMDomainLFGradKernel.md

@@ -0,0 +1,35 @@
+# MFEMDomainLFGradKernel
+
+!if! function=hasCapability('mfem')
+
+!syntax description /Kernels/MFEMDomainLFGradKernel
+
+## Overview
+
+Adds the domain integrator for integrating the linear form
+
+!equation
+(\vec{f}, \nabla v)_\Omega \,\,\, \forall v \in V
+
+where $v \in H^1$ is the test variable and $\vec{f}$ is a
+vector forcing coefficient.
+
+This term arises from the weak form of the forcing term
+
+!equation
+-\nabla \cdot \vec{f}
+
+## Example Input File Syntax
+
+!listing mfem/kernels/nldiffusion.i block=/Kernels
+
+!syntax parameters /Kernels/MFEMDomainLFGradKernel
+
+!syntax inputs /Kernels/MFEMDomainLFGradKernel
+
+!syntax children /Kernels/MFEMDomainLFGradKernel
+
+!if-end!
+
+!else
+!include mfem/mfem_warning.md

framework/doc/content/source/mfem/kernels/MFEMMixedScalarWeakDivergenceKernel.md

@@ -0,0 +1,35 @@
+# MFEMMixedScalarWeakDivergenceKernel
+
+!if! function=hasCapability('mfem')
+
+!syntax description /Kernels/MFEMMixedScalarWeakDivergenceKernel
+
+## Overview
+
+Adds the domain integrator for integrating the bilinear form
+
+!equation
+(-\vec{\lambda} u, \nabla v)_\Omega \,\,\, \forall v \in V
+
+where $v \in H^1$ is the test variable and $\vec{\lambda}$ is a
+vector forcing coefficient.
+
+This term arises from the weak form of the forcing term
+
+!equation
+\nabla \cdot(\vec{\lambda} u)
+
+## Example Input File Syntax
+
+!listing mfem/kernels/nldiffusion.i block=/Kernels
+
+!syntax parameters /Kernels/MFEMMixedScalarWeakDivergenceKernel
+
+!syntax inputs /Kernels/MFEMMixedScalarWeakDivergenceKernel
+
+!syntax children /Kernels/MFEMMixedScalarWeakDivergenceKernel
+
+!if-end!
+
+!else
+!include mfem/mfem_warning.md

framework/include/mfem/equation_systems/EquationSystem.h

@@ -32,7 +32,6 @@ class EquationSystem : public mfem::Operator
 
 public:
   friend class EquationSystemProblemOperator;
-  friend class TimeDomainEquationSystemProblemOperator;
 
   EquationSystem() = default;
   ~EquationSystem() override;
@@ -63,6 +62,8 @@ class EquationSystem : public mfem::Operator
   /// Update all essentially constrained true DoF markers and values on boundaries
   virtual void ApplyEssentialBCs();
 
+  virtual void BuildNonLinearActions();
+
   /// Perform trivial eliminations of coupled variables lacking corresponding test variables
   virtual void EliminateCoupledVariables();
 
@@ -71,7 +72,8 @@ class EquationSystem : public mfem::Operator
   /// Build mixed bilinear forms (off-diagonal Jacobian contributions)
   virtual void BuildMixedBilinearForms();
   /// Build all forms comprising this EquationSystem
-  virtual void BuildEquationSystem();
+  virtual void BuildEquationSystem(Moose::MFEM::GridFunctions & gridfunctions,
+                                   mfem::Array<int> & btoffsets);
 
   /// Form Jacobian operator based on on- and off-diagonal bilinear form contributions, populate
   /// solution and RHS vectors of true DoFs, and apply constraints
@@ -97,6 +99,8 @@ class EquationSystem : public mfem::Operator
                                 mfem::BlockVector & trueX,
                                 mfem::BlockVector & trueRHS);
 
+  void UpdateJacobian() const;
+
   /// Build linear system, with essential boundary conditions accounted for
   virtual void BuildJacobian(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS);
 
@@ -155,7 +159,7 @@ class EquationSystem : public mfem::Operator
   // Components of weak form. // Named according to test variable
   Moose::MFEM::NamedFieldsMap<mfem::ParBilinearForm> _blfs;
   Moose::MFEM::NamedFieldsMap<mfem::ParLinearForm> _lfs;
-  Moose::MFEM::NamedFieldsMap<mfem::ParNonlinearForm> _nlfs;
+  Moose::MFEM::NamedFieldsMap<mfem::ParLinearForm> _nlfs;
   Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>
       _mblfs; // named according to trial variable
 
@@ -177,6 +181,12 @@ class EquationSystem : public mfem::Operator
       Moose::MFEM::NamedFieldsMap<
           Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
 
+  void ApplyDomainNLActionIntegrators(
+      const std::string & test_var_name,
+      std::shared_ptr<mfem::ParLinearForm> form,
+      Moose::MFEM::NamedFieldsMap<
+          Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
+
   template <class FormType>
   void ApplyBoundaryBLFIntegrators(
       const std::string & trial_var_name,
@@ -212,8 +222,16 @@ class EquationSystem : public mfem::Operator
   Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMEssentialBC>>> _essential_bc_map;
 
   mutable mfem::OperatorHandle _jacobian;
+  mutable mfem::Vector _trueRHS;
+  mutable mfem::BlockVector _trueBlockRHS, _trueBlockSol, _BlockResidual;
+
+  Moose::MFEM::GridFunctions * _gfuncs;
+  mfem::Array<int> * _block_true_offsets;
+  mfem::Array<int> empty_tdof;
 
   mfem::AssemblyLevel _assembly_level;
+
+  bool _non_linear = false;
 };
 
 template <class FormType>
@@ -264,6 +282,30 @@ EquationSystem::ApplyDomainLFIntegrators(
   }
 }
 
+inline void
+EquationSystem::ApplyDomainNLActionIntegrators(
+    const std::string & test_var_name,
+    std::shared_ptr<mfem::ParLinearForm> form,
+    Moose::MFEM::NamedFieldsMap<
+        Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map)
+{
+  if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(test_var_name))
+  {
+    auto kernels = kernels_map.GetRef(test_var_name).GetRef(test_var_name);
+    for (auto & kernel : kernels)
+    {
+      mfem::LinearFormIntegrator * integ = kernel->createNLActionIntegrator();
+      if (integ != nullptr)
+      {
+        _non_linear = true;
+        kernel->isSubdomainRestricted()
+            ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
+            : form->AddDomainIntegrator(std::move(integ));
+      }
+    }
+  }
+}
+
 template <class FormType>
 void
 EquationSystem::ApplyBoundaryBLFIntegrators(
@@ -322,14 +364,16 @@ EquationSystem::ApplyBoundaryLFIntegrators(
 class TimeDependentEquationSystem : public EquationSystem
 {
 public:
+  friend class TimeDomainEquationSystemProblemOperator;
   TimeDependentEquationSystem(const Moose::MFEM::TimeDerivativeMap & time_derivative_map);
 
   /// Initialise
   virtual void Init(Moose::MFEM::GridFunctions & gridfunctions,
                     mfem::AssemblyLevel assembly_level) override;
 
   virtual void SetTimeStep(mfem::real_t dt);
-  virtual void UpdateEquationSystem();
+  virtual void UpdateEquationSystem(Moose::MFEM::GridFunctions & gridfunctions,
+                                    mfem::Array<int> & btoffsets);
 
   virtual void AddKernel(std::shared_ptr<MFEMKernel> kernel) override;
   virtual void BuildBilinearForms() override;

framework/include/mfem/executioners/MFEMProblemSolve.h

@@ -35,6 +35,10 @@ class MFEMProblemSolve : public SolveObject
 protected:
   MFEMProblem & _mfem_problem;
   std::vector<std::shared_ptr<Moose::MFEM::ProblemOperatorBase>> & _problem_operators;
+  unsigned int _nl_max_its;
+  mfem::real_t _nl_abs_tol;
+  mfem::real_t _nl_rel_tol;
+  unsigned int _print_level;
 };
 
 #endif

framework/include/mfem/kernels/MFEMDomainLFGradKernel.h

@@ -0,0 +1,34 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#ifdef MOOSE_MFEM_ENABLED
+
+#pragma once
+#include "MFEMKernel.h"
+
+/**
+ * \f[
+ * (f, v)
+ * \f]
+ */
+class MFEMDomainLFGradKernel : public MFEMKernel
+{
+public:
+  static InputParameters validParams();
+
+  MFEMDomainLFGradKernel(const InputParameters & parameters);
+
+  virtual mfem::LinearFormIntegrator * createNLActionIntegrator() override;
+
+protected:
+  mfem::Coefficient & _coef;
+  mfem::ScalarVectorProductCoefficient * _product_coeff;
+};
+
+#endif

framework/include/mfem/kernels/MFEMKernel.h

@@ -29,6 +29,7 @@ class MFEMKernel : public MFEMGeneralUserObject, public MFEMBlockRestrictable
 
   /// Create a new MFEM integrator to apply to the weak form. Ownership managed by the caller.
   virtual mfem::LinearFormIntegrator * createLFIntegrator() { return nullptr; }
+  virtual mfem::LinearFormIntegrator * createNLActionIntegrator() { return nullptr; }
   virtual mfem::BilinearFormIntegrator * createBFIntegrator() { return nullptr; }
 
   /// Get name of the test variable labelling the weak form this kernel is added to

framework/include/mfem/kernels/MFEMMixedScalarWeakDivergenceKernel.h

@@ -0,0 +1,34 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#ifdef MOOSE_MFEM_ENABLED
+
+#pragma once
+#include "MFEMKernel.h"
+
+/**
+ * \f[
+ * (k \vec \nabla u, \vec \nabla v)
+ * \f]
+ */
+class MFEMMixedScalarWeakDivergenceKernel : public MFEMKernel
+{
+public:
+  static InputParameters validParams();
+
+  MFEMMixedScalarWeakDivergenceKernel(const InputParameters & parameters);
+
+  virtual mfem::BilinearFormIntegrator * createBFIntegrator() override;
+
+protected:
+  mfem::Coefficient & _coef;
+  mfem::ScalarVectorProductCoefficient * _product_coef;
+};
+
+#endif

framework/include/mfem/problem/MFEMProblem.h

@@ -22,7 +22,6 @@ class MFEMProblem : public ExternalProblem
   MFEMProblem(const InputParameters & params);
   virtual ~MFEMProblem() {}
 
-  virtual void initialSetup() override;
   virtual void externalSolve() override {}
   virtual void syncSolutions(Direction) override {}
 
@@ -160,10 +159,12 @@ class MFEMProblem : public ExternalProblem
                      InputParameters & parameters);
 
   /**
-   * Add the nonlinear solver to the system. TODO: allow user to specify solver options,
-   * similar to the linear solvers.
+   * Add the nonlinear solver to the system.
    */
-  void addMFEMNonlinearSolver();
+  void addMFEMNonlinearSolver(unsigned int nl_max_its,
+                              mfem::real_t nl_abs_tol,
+                              mfem::real_t nl_rel_tol,
+                              unsigned int print_level);
 
   /**
    * Method used to get an mfem FEC depending on the variable family specified in the input file.