Unblock upcoming nonlinear solves and accelerate MFEM tests

framework/include/mfem/problem/MFEMProblem.h

@@ -188,7 +188,7 @@ class MFEMProblem : public ExternalProblem
   /**
    * Return the MPI communicator associated with this FE problem's mesh.
    */
-  MPI_Comm getComm() { return mesh().getMFEMParMesh().GetComm(); }
+  MPI_Comm getComm() { return getProblemData().comm; }
 
   /**
    * Displace the mesh, if mesh displacement is enabled.

framework/include/mfem/solvers/MFEMHyprePatch.h

@@ -0,0 +1,48 @@
+//* 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
+
+// To be removed as soon as we have a fix upstream in either hypre or mfem
+
+#ifdef MOOSE_MFEM_ENABLED
+
+namespace mfem
+{
+namespace patched
+{
+/**
+ * Patch for mfem::HypreGMRES to reset preconditioning matrix at every nonlinear/time iteration
+ */
+class HypreGMRES : public mfem::HypreGMRES
+{
+public:
+  using mfem::HypreGMRES::HypreGMRES;
+  void SetOperator(const mfem::Operator & op)
+  {
+    mfem::HypreGMRES::SetOperator(op);
+    HYPRE_GMRESSetPrecondMatrix(HYPRE_Solver(*this), nullptr);
+  }
+};
+
+/**
+ * Patch for mfem::HyprePCG to reset preconditioning matrix at every nonlinear/time iteration
+ */
+class HyprePCG : public mfem::HyprePCG
+{
+public:
+  using mfem::HyprePCG::HyprePCG;
+  void SetOperator(const mfem::Operator & op)
+  {
+    mfem::HyprePCG::SetOperator(op);
+    HYPRE_PCGSetPrecondMatrix(HYPRE_Solver(*this), nullptr);
+  }
+};
+}
+}
+
+#endif

framework/include/mfem/solvers/MFEMSolverBase.h

@@ -15,6 +15,7 @@
 #include <mfem.hpp>
 #include "libmesh/restore_warnings.h"
 #include <memory>
+#include "MFEMHyprePatch.h"
 
 /**
  * Base class for wrapping mfem::Solver-derived classes.

framework/src/mfem/problem/MFEMProblem.C

@@ -34,6 +34,8 @@ MFEMProblem::MFEMProblem(const InputParameters & params) : ExternalProblem(param
 {
   // Initialise Hypre for all MFEM problems.
   mfem::Hypre::Init();
+  // Disable multithreading for all MFEM problems (including any libMesh or MFEM subapps).
+  libMesh::libMeshPrivateData::_n_threads = 1;
   setMesh();
 }
 
@@ -50,8 +52,8 @@ MFEMProblem::setMesh()
   auto pmesh = mesh().getMFEMParMeshPtr();
   getProblemData().pmesh = pmesh;
   getProblemData().comm = pmesh->GetComm();
-  MPI_Comm_size(pmesh->GetComm(), &(getProblemData().num_procs));
-  MPI_Comm_rank(pmesh->GetComm(), &(getProblemData().myid));
+  getProblemData().num_procs = pmesh->GetNRanks();
+  getProblemData().myid = pmesh->GetMyRank();
 }
 
 void
@@ -76,7 +78,7 @@ MFEMProblem::addMFEMSolver(const std::string & user_object_name,
 void
 MFEMProblem::addMFEMNonlinearSolver()
 {
-  auto nl_solver = std::make_shared<mfem::NewtonSolver>(getProblemData().comm);
+  auto nl_solver = std::make_shared<mfem::NewtonSolver>(getComm());
 
   // Defaults to one iteration, without further nonlinear iterations
   nl_solver->SetRelTol(0.0);

framework/src/mfem/solvers/MFEMCGSolver.C

@@ -38,8 +38,7 @@ MFEMCGSolver::MFEMCGSolver(const InputParameters & parameters) : MFEMSolverBase(
 void
 MFEMCGSolver::constructSolver(const InputParameters &)
 {
-  auto solver =
-      std::make_unique<mfem::CGSolver>(getMFEMProblem().mesh().getMFEMParMesh().GetComm());
+  auto solver = std::make_unique<mfem::CGSolver>(getMFEMProblem().getComm());
   solver->SetRelTol(getParam<mfem::real_t>("l_tol"));
   solver->SetAbsTol(getParam<mfem::real_t>("l_abs_tol"));
   solver->SetMaxIter(getParam<int>("l_max_its"));

framework/src/mfem/solvers/MFEMGMRESSolver.C

@@ -38,8 +38,7 @@ MFEMGMRESSolver::MFEMGMRESSolver(const InputParameters & parameters) : MFEMSolve
 void
 MFEMGMRESSolver::constructSolver(const InputParameters &)
 {
-  auto solver =
-      std::make_unique<mfem::GMRESSolver>(getMFEMProblem().mesh().getMFEMParMesh().GetComm());
+  auto solver = std::make_unique<mfem::GMRESSolver>(getMFEMProblem().getComm());
   solver->SetRelTol(getParam<mfem::real_t>("l_tol"));
   solver->SetAbsTol(getParam<mfem::real_t>("l_abs_tol"));
   solver->SetMaxIter(getParam<int>("l_max_its"));

framework/src/mfem/solvers/MFEMHypreFGMRES.C

@@ -37,8 +37,7 @@ MFEMHypreFGMRES::MFEMHypreFGMRES(const InputParameters & parameters) : MFEMSolve
 void
 MFEMHypreFGMRES::constructSolver(const InputParameters &)
 {
-  auto solver =
-      std::make_unique<mfem::HypreFGMRES>(getMFEMProblem().mesh().getMFEMParMesh().GetComm());
+  auto solver = std::make_unique<mfem::HypreFGMRES>(getMFEMProblem().getComm());
   solver->SetTol(getParam<mfem::real_t>("l_tol"));
   solver->SetMaxIter(getParam<int>("l_max_its"));
   solver->SetKDim(getParam<int>("kdim"));
@@ -62,8 +61,7 @@ MFEMHypreFGMRES::updateSolver(mfem::ParBilinearForm & a, mfem::Array<int> & tdof
   {
     checkSpectralEquivalence(a);
     mfem::ParLORDiscretization lor_disc(a, tdofs);
-    auto lor_solver = new mfem::LORSolver<mfem::HypreFGMRES>(
-        lor_disc, getMFEMProblem().mesh().getMFEMParMesh().GetComm());
+    auto lor_solver = new mfem::LORSolver<mfem::HypreFGMRES>(lor_disc, getMFEMProblem().getComm());
     lor_solver->GetSolver().SetTol(getParam<mfem::real_t>("l_tol"));
     lor_solver->GetSolver().SetMaxIter(getParam<int>("l_max_its"));
     lor_solver->GetSolver().SetKDim(getParam<int>("kdim"));

framework/src/mfem/solvers/MFEMHypreGMRES.C

@@ -39,8 +39,7 @@ MFEMHypreGMRES::MFEMHypreGMRES(const InputParameters & parameters) : MFEMSolverB
 void
 MFEMHypreGMRES::constructSolver(const InputParameters &)
 {
-  auto solver =
-      std::make_unique<mfem::HypreGMRES>(getMFEMProblem().mesh().getMFEMParMesh().GetComm());
+  auto solver = std::make_unique<mfem::patched::HypreGMRES>(getMFEMProblem().getComm());
   solver->SetTol(getParam<mfem::real_t>("l_tol"));
   solver->SetAbsTol(getParam<mfem::real_t>("l_abs_tol"));
   solver->SetMaxIter(getParam<int>("l_max_its"));
@@ -65,8 +64,7 @@ MFEMHypreGMRES::updateSolver(mfem::ParBilinearForm & a, mfem::Array<int> & tdofs
   {
     checkSpectralEquivalence(a);
     mfem::ParLORDiscretization lor_disc(a, tdofs);
-    auto lor_solver = new mfem::LORSolver<mfem::HypreGMRES>(
-        lor_disc, getMFEMProblem().mesh().getMFEMParMesh().GetComm());
+    auto lor_solver = new mfem::LORSolver<mfem::HypreGMRES>(lor_disc, getMFEMProblem().getComm());
     lor_solver->GetSolver().SetTol(getParam<mfem::real_t>("l_tol"));
     lor_solver->GetSolver().SetAbsTol(getParam<mfem::real_t>("l_abs_tol"));
     lor_solver->GetSolver().SetMaxIter(getParam<int>("l_max_its"));

framework/src/mfem/solvers/MFEMHyprePCG.C

@@ -38,8 +38,7 @@ MFEMHyprePCG::MFEMHyprePCG(const InputParameters & parameters) : MFEMSolverBase(
 void
 MFEMHyprePCG::constructSolver(const InputParameters &)
 {
-  auto solver =
-      std::make_unique<mfem::HyprePCG>(getMFEMProblem().mesh().getMFEMParMesh().GetComm());
+  auto solver = std::make_unique<mfem::patched::HyprePCG>(getMFEMProblem().getComm());
   solver->SetTol(getParam<mfem::real_t>("l_tol"));
   solver->SetAbsTol(getParam<mfem::real_t>("l_abs_tol"));
   solver->SetMaxIter(getParam<int>("l_max_its"));
@@ -63,8 +62,7 @@ MFEMHyprePCG::updateSolver(mfem::ParBilinearForm & a, mfem::Array<int> & tdofs)
   {
     checkSpectralEquivalence(a);
     mfem::ParLORDiscretization lor_disc(a, tdofs);
-    auto lor_solver = new mfem::LORSolver<mfem::HyprePCG>(
-        lor_disc, getMFEMProblem().mesh().getMFEMParMesh().GetComm());
+    auto lor_solver = new mfem::LORSolver<mfem::HyprePCG>(lor_disc, getMFEMProblem().getComm());
     lor_solver->GetSolver().SetTol(getParam<mfem::real_t>("l_tol"));
     lor_solver->GetSolver().SetAbsTol(getParam<mfem::real_t>("l_abs_tol"));
     lor_solver->GetSolver().SetMaxIter(getParam<int>("l_max_its"));

framework/src/mfem/solvers/MFEMSolverBase.C

@@ -48,11 +48,6 @@ MFEMSolverBase::setPreconditioner(T & solver)
         mooseError("hypre solver preconditioners must themselves be hypre solvers");
     else
       solver.SetPreconditioner(mfem_pre);
-
-    if constexpr (std::is_same_v<mfem::HypreGMRES, T>)
-      HYPRE_GMRESSetPrecondMatrix(HYPRE_Solver(solver), NULL);
-    else if constexpr (std::is_same_v<mfem::HyprePCG, T>)
-      HYPRE_PCGSetPrecondMatrix(HYPRE_Solver(solver), NULL);
   }
 }
 
@@ -62,6 +57,9 @@ template void MFEMSolverBase::setPreconditioner(mfem::HypreFGMRES &);
 template void MFEMSolverBase::setPreconditioner(mfem::HypreGMRES &);
 template void MFEMSolverBase::setPreconditioner(mfem::HyprePCG &);
 
+template void MFEMSolverBase::setPreconditioner(mfem::patched::HypreGMRES &);
+template void MFEMSolverBase::setPreconditioner(mfem::patched::HyprePCG &);
+
 void
 MFEMSolverBase::checkSpectralEquivalence(mfem::ParBilinearForm & blf) const
 {

framework/src/mfem/solvers/MFEMSuperLU.C

@@ -32,8 +32,7 @@ MFEMSuperLU::MFEMSuperLU(const InputParameters & parameters) : MFEMSolverBase(pa
 void
 MFEMSuperLU::constructSolver(const InputParameters &)
 {
-  _solver = std::make_unique<Moose::MFEM::SuperLUSolver>(
-      getMFEMProblem().mesh().getMFEMParMesh().GetComm());
+  _solver = std::make_unique<Moose::MFEM::SuperLUSolver>(getMFEMProblem().getComm());
 }
 
 void

framework/src/mfem/submeshes/MFEMCutTransitionSubMesh.C

@@ -66,10 +66,8 @@ MFEMCutTransitionSubMesh::buildSubMesh()
 void
 MFEMCutTransitionSubMesh::labelMesh(mfem::ParMesh & parent_mesh)
 {
-  int mpi_comm_rank = 0;
-  int mpi_comm_size = 1;
-  MPI_Comm_rank(getMFEMProblem().getComm(), &mpi_comm_rank);
-  MPI_Comm_size(getMFEMProblem().getComm(), &mpi_comm_size);
+  int mpi_comm_rank = getMFEMProblem().getProblemData().myid;
+  int mpi_comm_size = getMFEMProblem().getProblemData().num_procs;
 
   // First determine face normal based on the first boundary element found on the cut
   // to use when determining orientation relative to the cut

test/tests/mfem/kernels/tests

@@ -14,27 +14,6 @@
     max_threads = 1
     platform = 'linux'
   []
-  [MFEMCurlCurlLOR]
-    type = XMLDiff
-    input = curlcurl.i
-    xmldiff = 'OutputData/CurlCurlLOR/Run0/Run0.pvd
-                OutputData/CurlCurlLOR/Run0/Cycle000001/proc000000.vtu'
-    cli_args = 'FESpaces/HCurlFESpace/fec_order=THIRD '
-               'FESpaces/HCurlFESpace/closed_basis=GaussLobatto '
-               'FESpaces/HCurlFESpace/open_basis=IntegratedGLL '
-               'Solver/type=MFEMGMRESSolver Solver/l_tol=1e-14 '
-               'Preconditioner/ams/low_order_refined=true '
-               'Outputs/ParaViewDataCollection/file_base=OutputData/CurlCurlLOR'
-    requirement = 'The system shall have the ability to solve a LOR definite Maxwell problem with Nedelec elements of the first kind using MFEM.'
-    capabilities = 'mfem'
-    compute_devices = 'cpu cuda'
-    valgrind = heavy
-    max_parallel = 1 # schemadiff with multiple ranks
-    recover = false
-    max_threads = 1
-    platform = 'linux'
-    abs_zero = 1e-6
-  []
   [MFEMDiffusion]
     type = XMLDiff
     input = diffusion.i
@@ -74,22 +53,6 @@
     max_parallel = 1
     recover = false
   []
-  [MFEMDiffusionLOR]
-    type = XMLDiff
-    input = diffusion.i
-    xmldiff = 'OutputData/DiffusionLOR/Run0/Run0.pvd '
-              'OutputData/DiffusionLOR/Run0/Cycle000001/proc000000.vtu'
-    cli_args = 'FESpaces/H1FESpace/fec_order=SECOND '
-               'Solver/type=MFEMCGSolver Solver/preconditioner=jacobi '
-               'Preconditioner/jacobi/low_order_refined=true '
-               'Outputs/ParaViewDataCollection/file_base=OutputData/DiffusionLOR'
-    requirement = 'The system shall have the ability to solve a diffusion problem with Low-Order-Refined preconditioning set up from MOOSE and produce the same result as a native MFEM run.'
-    capabilities = 'mfem'
-    compute_devices = 'cpu cuda'
-    max_parallel = 1 # schemadiff with multiple ranks
-    recover = false
-    abs_zero = 1e-6
-  []
   [MFEMGradDiv]
     type = XMLDiff
     input = graddiv.i
@@ -102,24 +65,6 @@
     recover = false
     errors = 'foobarbaz' # allow cuda errors I guess?
   []
-  [MFEMGradDivLOR]
-    type = XMLDiff
-    input = graddiv.i
-    xmldiff = 'OutputData/GradDivLOR/Run0/Run0.pvd
-                OutputData/GradDivLOR/Run0/Cycle000001/proc000000.vtu'
-    cli_args = 'Mesh/file=../mesh/beam-hex.mesh FESpaces/HDivFESpace/fec_order=THIRD '
-               'FESpaces/HDivFESpace/closed_basis=GaussLobatto '
-               'FESpaces/HDivFESpace/open_basis=IntegratedGLL '
-               'Preconditioner/ADS/low_order_refined=true '
-               'Outputs/ParaViewDataCollection/file_base=OutputData/GradDivLOR'
-    requirement = 'The system shall have the ability to solve a LOR grad-div problem with Raviart-Thomas elements using MFEM.'
-    capabilities = 'mfem'
-    compute_devices = 'cpu cuda'
-    valgrind = heavy
-    max_parallel = 1 # schemadiff with multiple ranks
-    recover = false
-    abs_zero = 1e-6
-  []
   [MFEMHeatConduction]
     type = XMLDiff
     input = heattransfer.i
@@ -175,25 +120,6 @@
     recover = false
     restep = false
   []
-  [MFEMHeatConductionLOR]
-    type = XMLDiff
-    input = heattransfer.i
-    xmldiff = 'OutputData/HeatConductionLOR/Run0/Run0.pvd
-                OutputData/HeatConductionLOR/Run0/Cycle000001/proc000000.vtu
-                OutputData/HeatConductionLOR/Run0/Cycle000004/proc000000.vtu'
-    cli_args = 'BCs/active="bottom top_dirichlet" '
-               'Executioner/dt=0.25 Executioner/end_time=1.0 '
-               'FESpaces/H1FESpace/fec_order=SECOND '
-               'Solver/type=MFEMCGSolver Solver/preconditioner=jacobi '
-               'Preconditioner/jacobi/low_order_refined=true '
-               'Outputs/ParaViewDataCollection/file_base=OutputData/HeatConductionLOR'
-    requirement = 'The system shall have the ability to solve a transient heat conduction problem with Low-Order-Refined preconditioning set up from MOOSE and produce the same result as a native MFEM run.'
-    capabilities = 'mfem'
-    compute_devices = 'cpu cuda'
-    max_parallel = 1 # schemadiff with multiple ranks
-    recover = false
-    restep = false
-  []
   [MFEMHeatTransfer]
     type = XMLDiff
     input = heattransfer.i