From 13285892276411d5a6f0c367cc813bcf0b50d16e Mon Sep 17 00:00:00 2001
From: Fan Fei <fei2@llnl.gov>
Date: Fri, 23 Jan 2026 13:00:55 -0800
Subject: [PATCH] enable the subsolver to be sequentially coupled

---
 .../multiphysics/CoupledSolver.hpp            | 81 ++++++++++++++++++-
 1 file changed, 77 insertions(+), 4 deletions(-)

diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp
index c6da8277ae7..d1211db5c6f 100644
--- a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp
+++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp
@@ -272,6 +272,52 @@ class CoupledSolver : public PhysicsSolverBase
 
   }
 
+  real64
+  nonlinearImplicitStep( real64 const & time_n,
+                         real64 const & dt,
+                         integer const cycleNumber,
+                         DomainPartition & domain ) override final
+  {
+    GEOS_MARK_FUNCTION;
+
+    if( getNonlinearSolverParameters().couplingType() == NonlinearSolverParameters::CouplingType::FullyImplicit )
+    {
+      return PhysicsSolverBase::nonlinearImplicitStep( time_n, dt, cycleNumber, domain );
+    }
+    else if( getNonlinearSolverParameters().couplingType() == NonlinearSolverParameters::CouplingType::Sequential )
+    {
+      return sequentiallyNonlinearImplicitStep( time_n, dt, cycleNumber, domain );
+    }
+    else
+    {
+      GEOS_ERROR( getDataContext() << ": Invalid coupling type option." );
+      return 0;
+    }
+  }
+
+  virtual void
+  setupSystem( DomainPartition & domain,
+               DofManager & dofManager,
+               CRSMatrix< real64, globalIndex > & localMatrix,
+               ParallelVector & rhs,
+               ParallelVector & solution,
+               bool const setSparsity = true ) override
+  {
+    GEOS_MARK_FUNCTION;
+
+    if( getNonlinearSolverParameters().couplingType() == NonlinearSolverParameters::CouplingType::FullyImplicit )
+    {
+      PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity );
+    }
+    else if( getNonlinearSolverParameters().couplingType() == NonlinearSolverParameters::CouplingType::Sequential )
+    {
+      sequentiallySetupSystem( domain );
+    }
+    else
+    {
+      GEOS_ERROR( getDataContext() << ": Invalid coupling type option." );
+    }
+  }
 
   virtual void
   updateAndWriteConvergenceStep( real64 const & time_n, real64 const & dt,
@@ -471,6 +517,21 @@ class CoupledSolver : public PhysicsSolverBase
   {
     GEOS_MARK_FUNCTION;
 
+    sequentiallySetupSystem( domain );
+
+    implicitStepSetup( time_n, dt, domain );
+
+    real64 const stepDt = sequentiallyNonlinearImplicitStep( time_n, dt, cycleNumber, domain );
+
+    implicitStepComplete( time_n, stepDt, domain );
+
+    return stepDt;
+  }
+
+  virtual void sequentiallySetupSystem( DomainPartition & domain )
+  {
+    GEOS_MARK_FUNCTION;
+
     // Only build the sparsity pattern if the mesh has changed
     Timestamp const meshModificationTimestamp = getMeshModificationTimestamp( domain );
     forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
@@ -485,9 +546,23 @@ class CoupledSolver : public PhysicsSolverBase
         solver->setSystemSetupTimestamp( meshModificationTimestamp );
       }
     } );
+  }
 
-    implicitStepSetup( time_n, dt, domain );
-
+  /**
+   * @brief Nonlinear implicit step for sequentially coupled solver. It solves a nonlinear system of
+   * equations using a sequential approach.
+   *
+   * @param time_n the current time
+   * @param dt timestep size
+   * @param cycleNumber
+   * @param domain the domain partition
+   * @return real64 size of the accepted timestep
+   */
+  virtual real64 sequentiallyNonlinearImplicitStep( real64 const & time_n,
+                                                    real64 const & dt,
+                                                    integer const cycleNumber,
+                                                    DomainPartition & domain )
+  {
     NonlinearSolverParameters & solverParams = getNonlinearSolverParameters();
     integer const maxNumberDtCuts = solverParams.m_maxTimeStepCuts;
     real64 const dtCutFactor = solverParams.m_timeStepCutFactor;
@@ -607,8 +682,6 @@ class CoupledSolver : public PhysicsSolverBase
       }
     }
 
-    implicitStepComplete( time_n, stepDt, domain );
-
     return stepDt;
   }