Fix for linear and non-linear problems #31772

Author: karthichockalingam

framework/src/mfem/equation_systems/EquationSystem.C

@@ -272,7 +272,7 @@ EquationSystem::FormLegacySystem(mfem::OperatorHandle & op,
     mfem::Vector aux_x, aux_rhs;
     mfem::HypreParMatrix * aux_a = new mfem::HypreParMatrix;
     blf->FormLinearSystem(
-        _ess_tdof_lists.at(i), *(_var_ess_constraints.at(i)), *lf, *aux_a, aux_x, aux_rhs, true);
+        _ess_tdof_lists.at(i), *(_var_ess_constraints.at(i)), *lf, *aux_a, aux_x, aux_rhs);
     _h_blocks(i, i) = aux_a;
     trueX.GetBlock(i) = aux_x;
     trueRHS.GetBlock(i) = aux_rhs;
@@ -339,44 +339,40 @@ EquationSystem::Mult(const mfem::Vector & sol, mfem::Vector & residual) const
 
   _BlockResidual = 0.0;
 
-  for (unsigned int i = 0; i < _test_var_names.size(); i++)
-  {
-    auto & test_var_name = _test_var_names.at(i);
-
-    int offset = _BlockResidual.GetBlock(i).Size();
-    mfem::Vector b(offset);
+ if (_non_linear)
+ {
+    for (unsigned int i = 0; i < _test_var_names.size(); i++)
+    {
+      auto & test_var_name = _test_var_names.at(i);
 
-    auto lf = _lfs.GetShared(test_var_name);
-    lf->Assemble();
-    lf->ParallelAssemble(b);
+      int offset = _BlockResidual.GetBlock(i).Size();
+      mfem::Vector b(offset);
+      
+      auto lf = _lfs.GetShared(test_var_name);
+      lf->Assemble();
+      lf->ParallelAssemble(b);
 
-    auto nlf = _nlfs.GetShared(test_var_name);
-    nlf->Assemble();
-    nlf->ParallelAssemble(_BlockResidual.GetBlock(i));
+      auto nlf = _nlfs.GetShared(test_var_name);
+      nlf->Assemble();
+      nlf->ParallelAssemble(_BlockResidual.GetBlock(i));
 
-    _BlockResidual.GetBlock(i) -= b;
-    _BlockResidual.GetBlock(i) *= -1;
+      _BlockResidual.GetBlock(i) -= b;
+      _BlockResidual.GetBlock(i) *= -1;
 
-    if (_non_linear)
       _BlockResidual.GetBlock(i).SetSubVector(_ess_tdof_lists.at(i), 0.0);
-  }
+    }
 
-  if (!_non_linear)
-  {
-    const_cast<EquationSystem *>(this)->FormLinearSystem(_jacobian, _trueBlockSol, _BlockResidual);
-    const_cast<EquationSystem *>(this)->CopyVec(_BlockResidual, residual);
-  }
-  else
-  {
     const_cast<EquationSystem *>(this)->CopyVec(_BlockResidual, residual);
     const_cast<EquationSystem *>(this)->FormLinearSystem(_jacobian, _trueBlockSol, _BlockResidual);
-  }
+ }
 
   residual *= -1.0;
 
   if (!_non_linear)
   {
-    _jacobian->AddMult(sol, residual);
+    //_jacobian->AddMult(sol, residual);
+    residual = 0.0;
+    _jacobian->Mult(sol, residual);
   }
 
   sol.HostRead();

framework/src/mfem/problem_operators/EquationSystemProblemOperator.C

@@ -54,7 +54,11 @@ EquationSystemProblemOperator::Solve()
   zero_vec = 0.0;
   _problem_data.nonlinear_solver->SetSolver(_problem_data.jacobian_solver->getSolver());
   _problem_data.nonlinear_solver->SetOperator(*GetEquationSystem());
-  _problem_data.nonlinear_solver->Mult(zero_vec, _true_x);
+
+  if(!(GetEquationSystem()->_non_linear))
+    _problem_data.nonlinear_solver->Mult(_true_rhs, _true_x);
+  else
+    _problem_data.nonlinear_solver->Mult(zero_vec, _true_x);
 
   GetEquationSystem()->RecoverFEMSolution(_true_x, _problem_data.gridfunctions);
 }

framework/src/mfem/problem_operators/TimeDomainEquationSystemProblemOperator.C

@@ -73,7 +73,13 @@ TimeDomainEquationSystemProblemOperator::ImplicitSolve(const mfem::real_t dt,
 
   _problem_data.nonlinear_solver->SetSolver(_problem_data.jacobian_solver->getSolver());
   _problem_data.nonlinear_solver->SetOperator(*GetEquationSystem());
-  _problem_data.nonlinear_solver->Mult(zero_vec, dX_dt);
+
+  if(!(GetEquationSystem()->_non_linear))
+    _problem_data.nonlinear_solver->Mult(_true_rhs, dX_dt);
+  else
+    _problem_data.nonlinear_solver->Mult(zero_vec, dX_dt);
+
+  //_problem_data.nonlinear_solver->Mult(zero_vec, dX_dt);
   SetTrialVariablesFromTrueVectors();
 }
 

test/tests/mfem/kernels/tests

@@ -34,6 +34,7 @@
     max_threads = 1
     platform = 'linux'
     abs_zero = 1e-6
+    max_time = 1000  
   []
   [MFEMDiffusion]
     type = XMLDiff