MFEM Non-linear support for steady state

framework/doc/content/source/mfem/equation_systems/EquationSystem.md

@@ -32,6 +32,18 @@ form. [`mfem::BilinearFormIntegrators`](https://mfem.org/bilininteg/) add contri
 $A_{ij}(\varphi_i, \phi_j)$ and [`mfem::LinearFormIntegrators`](https://mfem.org/lininteg/) add
 contributions to $b_i(\varphi_i)$ when assembled.
 
+[MFEMKernels](source/mfem/kernels/MFEMKernel.md) can also contribute to domain integrators for 
+non-linear actions. This allows to form the residual $\mathcal{L}(u)$ for non-linear Newton's 
+methood as shown below
+
+!equation
+{\mathbf{J}\left(\vec{u}_n\right) \delta \vec{u}_{n+1}=-\vec{R}\left(\vec{u}_n\right)}
+
+!equation
+{\vec{u}_{n+1}=\vec{u}_n+\delta \vec{u}_{n+1}}
+
+where $\mathbf{J}$ is the Jacobian, and $\delta \vec{u}$ is the incremental solution.
+
 !if-end!
 
 !else

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

@@ -22,7 +22,9 @@ variable names is the same as the set of trial variable names for a square syste
 
 `MFEMKernel` is a purely virtual base class. Derived classes should override the `createBFIntegrator`
 and/or the `createLFIntegrator` methods to return a `BilinearFormIntegrator` and/or a
-`LinearFormIntegrator` (respectively) to add to the `EquationSystem`.
+`LinearFormIntegrator` (respectively) to add to the `EquationSystem`. Derived classes could also 
+override `createNLAIntegrator` to solve non-linear problems using 
+[EquationSystem](source/mfem/equation_systems/EquationSystem.md).
 
 !if-end!
 

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

@@ -0,0 +1,37 @@
+# MFEMNLDiffusionKernel
+
+!if! function=hasCapability('mfem')
+
+## Overview
+
+Adds the domain integrator for integrating the non-linear action
+
+!equation
+(k(u)\vec\nabla v, \vec\nabla v)_\Omega \,\,\, \forall v \in V
+
+Adds the domain integrator for integrating the bilinear form
+
+!equation
+(k(u)\vec\nabla v, \vec\nabla v)_\Omega + (k'(u) v, \vec\nabla u \vec\nabla v)_\Omega \,\,\, \forall v \in V
+
+where $u, v \in H^1$ and $k(u)$ is a scalar non-linear diffusivity coefficient.
+
+The above terms arises from the weak form of the non-linear operator
+
+!equation
+- \vec\nabla \cdot \left( k(u) \vec\nabla u \right)
+
+## Example Input File Syntax
+
+!listing mfem/kernels/nldiffusion.i block=/Kernels
+
+!syntax parameters /Kernels/MFEMNLDiffusionKernel
+
+!syntax inputs /Kernels/MFEMNLDiffusionKernel
+
+!syntax children /Kernels/MFEMNLDiffusionKernel
+
+!if-end!
+
+!else
+!include mfem/mfem_warning.md

framework/include/mfem/bcs/MFEMIntegratedBC.h

@@ -24,6 +24,9 @@ class MFEMIntegratedBC : public MFEMBoundaryCondition
   /// Create MFEM integrator to apply to the RHS of the weak form. Ownership managed by the caller.
   virtual mfem::LinearFormIntegrator * createLFIntegrator() = 0;
 
+  /// Create MFEM integrator to apply non-linear residual form. Ownership managed by the caller.
+  virtual mfem::LinearFormIntegrator * createNLAIntegrator() { return nullptr; };
+
   /// Create MFEM integrator to apply to the LHS of the weak form. Ownership managed by the caller.
   virtual mfem::BilinearFormIntegrator * createBFIntegrator() = 0;
 

framework/include/mfem/equation_systems/ComplexEquationSystem.h

@@ -28,6 +28,9 @@ class ComplexEquationSystem : public EquationSystem
                     ComplexGridFunctions & cmplx_gridfunctions,
                     mfem::AssemblyLevel assembly_level) override;
 
+  ///Nonlinear Mult (Used by Newton-solver not necessarily nonlinear)
+  virtual void Mult(const mfem::Vector & x, mfem::Vector & y) const override;
+
   /// Build all forms comprising this EquationSystem
   virtual void BuildEquationSystem() override;
 

framework/include/mfem/equation_systems/EquationSystem.h

@@ -5,7 +5,7 @@
 //* 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
+//* https://www.gnu.org/licenses/lgpl-2.1.html*/
 
 #ifdef MOOSE_MFEM_ENABLED
 
@@ -85,17 +85,21 @@ class EquationSystem : public mfem::Operator
                                 mfem::Array<int> & global_ess_markers);
   /// Update all essentially constrained true DoF markers and values on boundaries
   virtual void ApplyEssentialBCs();
-
   /// Perform trivial eliminations of coupled variables lacking corresponding test variables
   virtual void EliminateCoupledVariables();
   /// Build linear forms and eliminate constrained DoFs
   virtual void BuildLinearForms();
+  /// Build non-linear action forms
+  virtual void BuildNonLinearActionForms();
   /// Build bilinear forms (diagonal Jacobian contributions)
   virtual void BuildBilinearForms();
   /// Build mixed bilinear forms (off-diagonal Jacobian contributions)
   virtual void BuildMixedBilinearForms();
   /// Build all forms comprising this EquationSystem
   virtual void BuildEquationSystem();
+  // Reassemble the Jacobian matrix/Operator (Assuming Gridfunction dependant Operator)
+  virtual void ReassembleJacobian(mfem::BlockVector & x, mfem::BlockVector & rhs);
+  virtual void UpdateJacobian() const;
 
   /// Form linear system and jacobian operator based on on- and off-diagonal bilinear form
   /// contributions, populate solution and RHS vectors of true DoFs, and apply constraints.
@@ -130,6 +134,11 @@ class EquationSystem : public mfem::Operator
       std::shared_ptr<mfem::ParLinearForm> form,
       NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
 
+  void ApplyDomainNLAFIntegrators(
+      const std::string & test_var_name,
+      std::shared_ptr<mfem::ParLinearForm> form,
+      NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
+
   template <class FormType>
   void ApplyBoundaryBLFIntegrators(
       const std::string & trial_var_name,
@@ -145,6 +154,12 @@ class EquationSystem : public mfem::Operator
       NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
           integrated_bc_map);
 
+  void ApplyBoundaryNLAFIntegrators(
+      const std::string & test_var_name,
+      std::shared_ptr<mfem::ParLinearForm> form,
+      NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
+          integrated_bc_map);
+
   /// Names of all trial variables of kernels and boundary conditions
   /// added to this EquationSystem.
   std::vector<std::string> _coupled_var_names;
@@ -167,8 +182,8 @@ class EquationSystem : public mfem::Operator
   NamedFieldsMap<mfem::ParBilinearForm> _blfs;
   NamedFieldsMap<mfem::ParLinearForm> _lfs;
   NamedFieldsMap<mfem::ParNonlinearForm> _nlfs;
-  NamedFieldsMap<NamedFieldsMap<mfem::ParMixedBilinearForm>>
-      _mblfs; // named according to trial variable
+  NamedFieldsMap<mfem::ParLinearForm> _nlAs;
+  NamedFieldsMap<NamedFieldsMap<mfem::ParMixedBilinearForm>> _mblfs; // named according to trial var
 
   /// Gridfunctions holding essential constraints from Dirichlet BCs
   std::vector<std::unique_ptr<mfem::ParGridFunction>> _var_ess_constraints;
@@ -185,10 +200,18 @@ class EquationSystem : public mfem::Operator
   /// Named according to test variable.
   NamedFieldsMap<std::vector<std::shared_ptr<MFEMEssentialBC>>> _essential_bc_map;
 
+  // Operator handle for the jacobian matrix
   mutable mfem::OperatorHandle _jacobian;
-
   mfem::AssemblyLevel _assembly_level;
 
+  // Temporary vectors used for non-linear action
+  // assembly process
+  mutable mfem::BlockVector _trueBlockSol, _blockForces, _blockResidual;
+  Moose::MFEM::GridFunctions * _gfuncs;
+  mfem::Array<int> * _block_true_offsets = NULL;
+  mfem::Array<int> empty_tdof;
+  bool _non_linear = false;
+
 private:
   friend class EquationSystemProblemOperator;
   /// Disallowed inherited method
@@ -246,6 +269,29 @@ EquationSystem::ApplyDomainLFIntegrators(
   }
 }
 
+inline void
+EquationSystem::ApplyDomainNLAFIntegrators(
+    const std::string & test_var_name,
+    std::shared_ptr<mfem::ParLinearForm> form,
+    NamedFieldsMap<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->createNLAIntegrator();
+      if (integ)
+      {
+        _non_linear = true;
+        kernel->isSubdomainRestricted()
+            ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
+            : form->AddDomainIntegrator(std::move(integ));
+      }
+    }
+  }
+}
+
 template <class FormType>
 void
 EquationSystem::ApplyBoundaryBLFIntegrators(
@@ -301,6 +347,31 @@ EquationSystem::ApplyBoundaryLFIntegrators(
   }
 }
 
+inline void
+EquationSystem::ApplyBoundaryNLAFIntegrators(
+    const std::string & test_var_name,
+    std::shared_ptr<mfem::ParLinearForm> form,
+    NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
+        integrated_bc_map)
+{
+  if (integrated_bc_map.Has(test_var_name) &&
+      integrated_bc_map.Get(test_var_name)->Has(test_var_name))
+  {
+    auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(test_var_name);
+    for (auto & bc : bcs)
+    {
+      mfem::LinearFormIntegrator * integ = bc->createNLAIntegrator();
+      if (integ)
+      {
+        _non_linear = true;
+        bc->isBoundaryRestricted()
+            ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
+            : form->AddBoundaryIntegrator(std::move(integ));
+      }
+    }
+  }
+}
+
 } // namespace Moose::MFEM
 
 #endif

framework/include/mfem/executioners/MFEMProblemSolve.h

@@ -36,6 +36,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/MFEMKernel.h

@@ -30,6 +30,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::BilinearFormIntegrator * createBFIntegrator() { return nullptr; }
+  virtual mfem::LinearFormIntegrator * createNLAIntegrator() { return nullptr; }
 
   /// Get name of the test variable labelling the weak form this kernel is added to
   const VariableName & getTestVariableName() const { return _test_var_name; }

framework/include/mfem/kernels/MFEMNLDiffusionKernel.h

@@ -0,0 +1,38 @@
+//* 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(u) \vec \nabla u, \vec \nabla v)
+ * \f]
+ */
+class MFEMNLDiffusionKernel : public MFEMKernel
+{
+public:
+  static InputParameters validParams();
+
+  MFEMNLDiffusionKernel(const InputParameters & parameters);
+
+  virtual mfem::BilinearFormIntegrator * createBFIntegrator() override;
+  virtual mfem::LinearFormIntegrator * createNLAIntegrator() override;
+
+protected:
+  mfem::Coefficient & _coef;
+  mfem::ScalarVectorProductCoefficient * _product_coef_res;
+  mfem::ScalarVectorProductCoefficient * _product_coef_jac;
+  mfem::SumIntegrator * _sum;
+  mfem::ConstantCoefficient * _one;
+};
+
+#endif

framework/include/mfem/problem/MFEMProblem.h

@@ -5,7 +5,7 @@
 //* 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
+//* https://www.gnu.org/licenses/lgpl-2.1.html*/
 
 #ifdef MOOSE_MFEM_ENABLED
 
@@ -192,7 +192,10 @@ class MFEMProblem : public ExternalProblem
    * Add the nonlinear solver to the system. TODO: allow user to specify solver options,
    * similar to the linear solvers.
    */
-  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.

framework/src/mfem/equation_systems/ComplexEquationSystem.C

@@ -269,6 +269,15 @@ ComplexEquationSystem::FormSystemMatrix(mfem::OperatorHandle & op,
   op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
 }
 
+// Equation system Mult
+void
+ComplexEquationSystem::Mult(const mfem::Vector & x, mfem::Vector & residual) const
+{
+  _jacobian->Mult(x, residual);
+  x.HostRead();
+  residual.HostRead();
+}
+
 void
 ComplexEquationSystem::RecoverComplexFEMSolution(
     mfem::BlockVector & trueX,