fixed cmake

Merge branch 'master' into multiapp
Fixed const correctness
2026-01-06 18:59:44 -08:00 · 2026-01-06 18:53:05 -08:00 · 2026-01-06 18:49:38 -08:00 · 2025-11-19 14:28:38 -08:00 · 2025-11-19 14:24:23 -08:00 · 2025-11-19 14:24:04 -08:00
21 changed files with 6105 additions and 17 deletions
@@ -2137,14 +2137,14 @@ private:
   int trial_dofs1D, test_dofs1D, quad1D;

 public:
-   GradientIntegrator() :
-      Q{NULL}, trial_maps{NULL}, test_maps{NULL}, geom{NULL}
+   GradientIntegrator(const IntegrationRule *ir = NULL) :
+      BilinearFormIntegrator(ir), Q{NULL}, trial_maps{NULL}, test_maps{NULL}, geom{NULL}
   { }
-   GradientIntegrator(Coefficient *q_) :
-      Q{q_}, trial_maps{NULL}, test_maps{NULL}, geom{NULL}
+   GradientIntegrator(Coefficient *q_, const IntegrationRule *ir = NULL) :
+      BilinearFormIntegrator(ir), Q{q_}, trial_maps{NULL}, test_maps{NULL}, geom{NULL}
   { }
-   GradientIntegrator(Coefficient &q) :
-      Q{&q}, trial_maps{NULL}, test_maps{NULL}, geom{NULL}
+   GradientIntegrator(Coefficient &q, const IntegrationRule *ir = NULL) :
+      BilinearFormIntegrator(ir), Q{&q}, trial_maps{NULL}, test_maps{NULL}, geom{NULL}
   { }

   void AssembleElementMatrix2(const FiniteElement &trial_fe,
@@ -2612,6 +2612,11 @@ protected:
   Vector pa_data;

 public:
+
+   VectorMassIntegrator(const IntegrationRule *ir)
+      : BilinearFormIntegrator(ir), vdim(-1), Q_order(0), Q(NULL), VQ(NULL),
+        MQ(NULL) { }
+
   /// Construct an integrator with coefficient 1.0
   VectorMassIntegrator() = default;

@@ -3038,14 +3043,15 @@ private:
   int trial_dofs1D, test_dofs1D, quad1D;

 public:
-   VectorDivergenceIntegrator() :
-      Q(NULL), trial_maps(NULL), test_maps(NULL), geom(NULL)
+   VectorDivergenceIntegrator(const IntegrationRule *ir = NULL) :
+      BilinearFormIntegrator(ir), Q(NULL), trial_maps(NULL), test_maps(NULL),
+      geom(NULL)
   {  }
-   VectorDivergenceIntegrator(Coefficient *q_) :
-      Q(q_), trial_maps(NULL), test_maps(NULL), geom(NULL)
+   VectorDivergenceIntegrator(Coefficient *q_, const IntegrationRule *ir = NULL) :
+      BilinearFormIntegrator(ir), Q(q_), trial_maps(NULL), test_maps(NULL), geom(NULL)
   { }
-   VectorDivergenceIntegrator(Coefficient &q) :
-      Q(&q), trial_maps(NULL), test_maps(NULL), geom(NULL)
+   VectorDivergenceIntegrator(Coefficient &q, const IntegrationRule *ir = NULL) :
+      BilinearFormIntegrator(ir), Q(&q), trial_maps(NULL), test_maps(NULL), geom(NULL)
   { }

   void AssembleElementMatrix2(const FiniteElement &trial_fe,
@@ -3263,8 +3269,9 @@ private:
   void SetUpQuadratureSpaceAndCoefficients(const FiniteElementSpace &fes);

 public:
-   ElasticityIntegrator(Coefficient &l, Coefficient &m)
-   { lambda = &l; mu = &m; }
+   ElasticityIntegrator(Coefficient &l, Coefficient &m,
+                        const IntegrationRule *ir = NULL) :
+      BilinearFormIntegrator(ir) { lambda = &l; mu = &m; }
   /** With this constructor $\lambda = q_l m$ and $\mu = q_m m$
       if $dim q_l + 2 q_m = 0$ then $tr(\sigma) = 0$. */
   ElasticityIntegrator(Coefficient &m, real_t q_l, real_t q_m)
@@ -652,6 +652,16 @@ void forall(func_t f,
 class FDJacobian : public Operator
 {
 public:
+   FDJacobian(const Operator &op, real_t fixed_eps = 0.0) :
+      Operator(op.Height(), op.Width()),
+      op(op),
+      x(op.Width()),
+      fixed_eps(fixed_eps)
+   {
+      f.SetSize(Height());
+      xpev.SetSize(Width());
+   }
+
   FDJacobian(const Operator &op, const Vector &x, real_t fixed_eps = 0.0) :
      Operator(op.Height(), op.Width()),
      op(op),
@@ -671,6 +681,18 @@ public:
                    MPI_COMM_WORLD);
   }

+   void Update(const Vector &x_new)
+   {
+      x = x_new;
+      f.UseDevice(x.UseDevice());
+      xpev.UseDevice(x.UseDevice());
+
+      op.Mult(x, f);
+      const real_t xnorm_local = x.Norml2();
+      MPI_Allreduce(&xnorm_local, &xnorm, 1, MPITypeMap<real_t>::mpi_type, MPI_SUM,
+                    MPI_COMM_WORLD);
+   }
+
   void Mult(const Vector &v, Vector &y) const override
   {
      // See [1] for choice of eps.
@@ -758,7 +758,7 @@ void VectorConvectionNLFIntegrator::AssembleElementVector(
   EF.UseExternalData(elfun.GetData(), nd, dim);
   ELV.UseExternalData(elvect.GetData(), nd, dim);

-   Vector vec1(dim), vec2(dim);
+   Vector vec1(dim), vec2(dim), vecvq(dim);
   const IntegrationRule *ir = GetIntegrationRule(el, T);
   ELV = 0.0;
   for (int i = 0; i < ir->GetNPoints(); i++)
@@ -772,6 +772,11 @@ void VectorConvectionNLFIntegrator::AssembleElementVector(

      MultAtB(EF, dshape, gradEF);
      EF.MultTranspose(shape, vec1);
+      if (VQ)
+      {
+         VQ->Eval(vecvq, T, ip);
+         vec1 += vecvq;
+      }
      gradEF.Mult(vec1, vec2);
      vec2 *= w;
      AddMultVWt(shape, vec2, ELV);
@@ -797,7 +802,7 @@ void VectorConvectionNLFIntegrator::AssembleElementGrad(
   EF.UseExternalData(elfun.GetData(), nd, dim);

   real_t w;
-   Vector vec1(dim), vec2(dim), vec3(nd);
+   Vector vec1(dim), vec2(dim), vec3(nd), vecvq(dim);

   const IntegrationRule *ir = GetIntegrationRule(el, trans);

@@ -822,6 +827,12 @@ void VectorConvectionNLFIntegrator::AssembleElementGrad(
      MultAtB(EF, dshapex, gradEF);
      EF.MultTranspose(shape, vec1);

+      if (VQ)
+      {
+         VQ->Eval(vecvq, trans, ip);
+         vec1 += vecvq;
+      }
+
      trans.AdjugateJacobian().Mult(vec1, vec2);

      vec2 *= w;
@@ -381,6 +381,7 @@ class VectorConvectionNLFIntegrator : public NonlinearFormIntegrator
 {
 private:
   Coefficient *Q{};
+   VectorCoefficient *VQ{};
   DenseMatrix dshape, dshapex, EF, gradEF, ELV, elmat_comp;
   Vector shape;
   // PA extension
@@ -390,7 +391,13 @@ private:
   int dim, ne, nq;

 public:
-   VectorConvectionNLFIntegrator(Coefficient &q): Q(&q) { }
+   VectorConvectionNLFIntegrator(Coefficient &q,
+                                 const IntegrationRule *ir = NULL) :
+      NonlinearFormIntegrator(ir), Q(&q) { }
+
+   VectorConvectionNLFIntegrator(Coefficient &q, VectorCoefficient &vq,
+                                 const IntegrationRule *ir = NULL) :
+      NonlinearFormIntegrator(ir), Q(&q), VQ(&vq) { }

   VectorConvectionNLFIntegrator() = default;

@@ -193,6 +193,9 @@ public:
   /// Returns how many State vectors the ODE requires
   virtual int GetStateSize() { return 0; };

+   /// Returns the associated TimeDependentOperator
+   TimeDependentOperator* GetTimeDependentOperator() { return f; }
+
   ///@brief Returns @a true if the ODESolver supports the given
   /// #ImplicitVariableType, @a var, and returns @a false otherwise.
   ///@note Should be overriden in ODESolver that calls TimeDependentOperator::ImplicitSolve().
@@ -851,6 +851,216 @@ void SLI(const Operator &A, Solver &B, const Vector &b, Vector &x,
   sli.Mult(b, x);
 }

+real_t FPIRelaxation::Dot(const Vector &x, const Vector &y) const
+{
+   if (dot_oper) { return dot_oper->Eval(x,y); } // Use custom inner product (if provided)
+#ifndef MFEM_USE_MPI
+   return (x * y);
+#else
+   return InnerProduct(comm, x, y);
+#endif
+}
+
+real_t FPIRelaxation::SquaredDistance(const Vector &x, const Vector &y) const
+{
+#ifndef MFEM_USE_MPI
+   return x.DistanceSquaredTo(y);
+#else
+   return DistanceSquared(comm, x, y);
+#endif
+}
+
+void AitkenRelaxation::Init()
+{
+   FPIRelaxation::Init();
+   MemoryType mt = GetMemoryType(oper->GetMemoryClass());
+   rold.SetSize(oper->Width(), mt);
+   rold.UseDevice(true);
+   rold = 0.0;
+   rold_nsq = 0.0;
+}
+
+real_t AitkenRelaxation::Eval(const Vector &state, const Vector &residual,
+                              real_t res_norm, real_t rfactor)
+{
+   real_t rdot  = Dot(rold, residual);
+   real_t rnsq  = res_norm * res_norm;
+   real_t denom = rnsq - 2.0*rdot + rold_nsq; // ||rold - residual||^2
+   real_t ratio = (rdot - rold_nsq) / denom;
+
+   rold_nsq = rnsq;
+   rold     = residual;
+   if (ratio == 0.0) { return rfactor; } // Avoid factor = 0.0 at first call
+   return Clamp(-rfactor * ratio);
+}
+
+void SteepestDescentRelaxation::Init()
+{
+   FPIRelaxation::Init();
+   MemoryType mt = GetMemoryType(oper->GetMemoryClass());
+   z.SetSize(oper->Width(), mt);
+   z.UseDevice(true);
+}
+
+real_t SteepestDescentRelaxation::Eval(const Vector &state,
+                                       const Vector &residual,
+                                       real_t res_norm, real_t rfactor)
+{
+   MFEM_VERIFY(oper,"Operator not set; set using SetOperator(Operator&)")
+
+   Operator *J = &oper->GetGradient(state);
+   real_t num = res_norm * res_norm;
+   J->Mult(residual, z); // rold = F'(x) * rnew;
+   real_t denom = Dot(z, residual);
+   return Clamp(num/denom);
+}
+
+/// Fixed point iteration solver: x <- f(x)
+void FPISolver::UpdateVectors()
+{
+   MemoryType mt = GetMemoryType(oper->GetMemoryClass());
+
+   r.SetSize(width, mt);
+   r.UseDevice(true);
+
+   z.SetSize(width, mt);
+   z.UseDevice(true);
+}
+
+void FPISolver::SetOperator(const Operator &op)
+{
+   IterativeSolver::SetOperator(op);
+   UpdateVectors();
+   if (!relax_method)
+   {
+      relax_method = new FPIRelaxation(); // Default relaxation strategy
+      relax_owned = true;
+   }
+   relax_method->SetOperator(*oper);
+}
+
+void FPISolver::SetRelaxation(real_t rfactor, FPIRelaxation *relaxation,
+                              bool own)
+{
+   relax_factor = rfactor;
+   if (relaxation)
+   {
+      if (relax_method && relax_owned) { delete relax_method; }
+      relax_method = relaxation;
+      relax_owned  = own;
+   }
+#ifdef MFEM_USE_MPI
+   if (relax_method) { relax_method->SetComm(this->GetComm()); }
+#endif
+}
+
+/// Iterative solution of the (non)linear system using Fixed Point Iteration
+void FPISolver::Mult(const Vector &b, Vector &x) const
+{
+   int i;
+   real_t factor = relax_factor;
+   real_t r0, nom, nom0, nomold = 1, cf, fac_old = factor;
+
+   if (iterative_mode)
+   {
+      oper->Mult(x, z);  // z = F(x)
+      subtract(z, x, r); // r = z - x
+   }
+   else
+   {
+      x = 0.0;
+      oper->Mult(x, r);  // r = F(x)
+   }
+
+   nom0 = nom = sqrt(Dot(r, r));
+   initial_norm = nom0;
+
+   if (print_options.iterations | print_options.first_and_last)
+   {
+      mfem::out << "   Iteration : " << setw(3) << right << 0 << "  ||Br|| = "
+                << nom << (print_options.first_and_last ? " ..." : "") << '\n';
+   }
+
+   r0 = std::max(nom*rel_tol, abs_tol);
+   if (nom <= r0)
+   {
+      converged = true;
+      final_iter = 0;
+      final_norm = nom;
+      return;
+   }
+
+   // start iteration
+   converged = false;
+   final_iter = max_iter;
+   relax_method->Init();
+   for (i = 1; true; )
+   {
+      factor = relax_method->Eval(x,r,nom,fac_old);
+      x.Add(factor, r); // x = x + factor * r
+
+      oper->Mult(x, z);  // z = F(x)
+      subtract(z, x, r); // r = z - x
+      nom = sqrt(Dot(r, r));
+
+      cf = nom/nomold;
+      nomold = nom;
+      fac_old = factor;
+
+      bool done = false;
+      if (nom < r0)
+      {
+         converged = true;
+         final_iter = i;
+         done = true;
+      }
+
+      if (++i > max_iter)
+      {
+         done = true;
+      }
+
+      if (print_options.iterations || (done && print_options.first_and_last))
+      {
+         mfem::out << "   Iteration : " << setw(3) << right << (i-1)
+                   << "  ||r|| = " << setw(11) << left << nom
+                   << "\trlx. fac.: " << fac_old << '\n';
+      }
+
+      if (done) { break; }
+   }
+
+   if (print_options.summary || (print_options.warnings && !converged))
+   {
+      const auto rf = pow (nom/nom0, 1.0/final_iter);
+      mfem::out << "FPI: Number of iterations: " << final_iter << '\n'
+                << "Conv. rate: " << cf << '\n'
+                << "Average reduction factor: "<< rf << '\n';
+   }
+   if (print_options.warnings && !converged)
+   {
+      mfem::out << "FPI: No convergence!" << '\n';
+   }
+
+   final_norm = nom;
+}
+
+void FPI(const Operator &A, const Vector &b, Vector &x,
+         int print_iter, int max_num_iter,
+         real_t RTOLERANCE, real_t ATOLERANCE,
+         real_t relax_factor, FPIRelaxation *relax_method)
+{
+   MFEM_PERF_FUNCTION;
+
+   FPISolver fpi;
+   fpi.SetPrintLevel(print_iter);
+   fpi.SetMaxIter(max_num_iter);
+   fpi.SetRelTol(sqrt(RTOLERANCE));
+   fpi.SetAbsTol(sqrt(ATOLERANCE));
+   fpi.SetOperator(A);
+   fpi.SetRelaxation(relax_factor,relax_method,false);
+   fpi.Mult(b, x);
+}

 void CGSolver::UpdateVectors()
 {
@@ -622,6 +622,184 @@ void SLI(const Operator &A, Solver &B, const Vector &b, Vector &x,
         real_t RTOLERANCE = 1e-12, real_t ATOLERANCE = 1e-24);


+
+/**
+   @brief A class to handle fixed point iteration relaxation methods.
+   This class provides a base for implementing various relaxation strategies
+   for fixed point iteration solvers.
+ */
+class FPIRelaxation
+{
+#ifdef MFEM_USE_MPI
+private:
+   MPI_Comm comm = MPI_COMM_NULL;
+#endif
+
+protected:
+   const Operator *oper = nullptr;
+   real_t lbnd = std::numeric_limits<real_t>::lowest();
+   real_t ubnd = std::numeric_limits<real_t>::max();
+   real_t abs_lbnd = 0.0;
+   InnerProductOperator *dot_oper = nullptr;
+
+public:
+   FPIRelaxation() = default;
+
+#ifdef MFEM_USE_MPI
+   FPIRelaxation(MPI_Comm comm_) {comm = comm_;}
+   void SetComm(MPI_Comm comm_) {comm = comm_;}
+#endif
+
+   /// @brief Set the operator for the relaxation method.
+   virtual void SetOperator(const Operator &op) {oper = &op;}
+
+   /// @brief Initialize the relaxation method.
+   virtual void Init() {}
+
+   /// @brief Set the lower bound for the relaxation factor.
+   virtual void SetLowerBound(real_t lb) { lbnd = lb; }
+
+   /// @brief Set the upper bound for the relaxation factor.
+   virtual void SetUpperBound(real_t ub) { ubnd = ub; }
+
+   /// @brief Set both lower and upper bounds for the relaxation factor.
+   virtual void SetBounds(real_t lb, real_t ub)
+   { lbnd = lb; ubnd = ub;}
+
+   /// @brief Set the absolute bound for the relaxation factor.
+   virtual void SetAbsoluteLowerBound(real_t abs_lb)
+   { abs_lbnd = std::abs(abs_lb); }
+
+   /// @brief Set a user-defined inner product operator (not owned)
+   void SetInnerProduct(InnerProductOperator *ipo) { dot_oper = ipo; }
+
+   /// @brief Clamp the relaxation factor to the specified range.
+   virtual real_t Clamp(real_t factor) const
+   {
+      real_t afac = (abs_lbnd == 0.0) ? factor :
+                    std::copysign(std::max(std::abs(factor), abs_lbnd), factor);
+      return std::max(std::min(afac, ubnd), lbnd);
+   }
+
+   /**
+    @brief Compute the relaxation factor for Fixed Point Iteration.
+
+      @param state Current state vector
+      @param residual Current residual vector
+      @param res_norm Norm of the current residual vector
+      @param rfactor Current relaxation factor
+      @return real_t The computed relaxation factor
+   */
+   virtual real_t Eval(const Vector &state, const Vector &residual,
+                       real_t res_norm, real_t rfactor)
+   {
+      return rfactor; // Default implementation returns the fixed factor
+   };
+
+   real_t Dot(const Vector &x, const Vector &y) const;
+
+   real_t SquaredDistance(const Vector &x, const Vector &y) const;
+
+   virtual ~FPIRelaxation() {}
+};
+
+/**
+   @brief A class to implement Aitken relaxation for Fixed Point Iteration
+   to accelerate the convergence of fixed point iterations.
+ */
+class AitkenRelaxation : public FPIRelaxation
+{
+protected:
+   Vector rold; // Old residual vector
+   real_t rold_nsq = 0.0;
+
+public:
+   AitkenRelaxation() = default;
+
+#ifdef MFEM_USE_MPI
+   AitkenRelaxation(MPI_Comm comm_) : FPIRelaxation(comm_) {}
+#endif
+
+   void Init() override;
+
+   /**
+    @brief Compute the Aitken relaxation factor for Fixed Point Iteration.
+
+      @param state Current state vector
+      @param residual Current residual vector
+      @param res_norm Norm of the current residual vector
+      @param rfactor Current relaxation factor
+      @return The computed relaxation factor
+   */
+   real_t Eval(const Vector &state, const Vector &residual, real_t res_norm,
+               real_t rfactor) override;
+};
+
+/**
+   @brief A class to implement steepest descent relaxation for Fixed Point Iteration
+   to accelerate the convergence of fixed point iterations.
+ */
+class SteepestDescentRelaxation : public FPIRelaxation
+{
+protected:
+   Vector z;
+
+public:
+   SteepestDescentRelaxation() = default;
+
+#ifdef MFEM_USE_MPI
+   SteepestDescentRelaxation(MPI_Comm comm_) : FPIRelaxation(comm_) {}
+#endif
+
+   void Init() override;
+
+   /**
+      @brief Compute the steepest descent relaxation factor for Fixed Point Iteration.
+
+      @param state Current state vector
+      @param residual Current residual vector
+      @param res_norm Norm of the current residual vector
+      @param rfactor Current relaxation factor
+      @return The computed relaxation factor
+    */
+   real_t Eval(const Vector &state, const Vector &residual, real_t res_norm,
+               real_t rfactor) override;
+};
+
+/// Fixed point iteration solver: x <- f(x)
+class FPISolver : public IterativeSolver
+{
+protected:
+   FPIRelaxation *relax_method = nullptr; ///< Relaxation strategy for FPI
+   bool relax_owned = false;
+   mutable real_t relax_factor = 1.0;
+   mutable Vector r, z;
+
+   void UpdateVectors();
+
+public:
+   FPISolver() { }
+
+#ifdef MFEM_USE_MPI
+   FPISolver(MPI_Comm comm_) : IterativeSolver(comm_) { }
+#endif
+
+   virtual void SetOperator(const Operator &op) override;
+
+   void SetRelaxation(real_t rfactor, FPIRelaxation *relaxation = nullptr,
+                      bool own = false);
+
+   /// Iterative solution of the (non)linear system using Fixed Point Iteration
+   /// b is not used
+   void Mult(const Vector &b, Vector &x) const override;
+};
+
+/// Fixed point iteration. (tolerances are squared)
+void FPI(const Operator &A, const Vector &b, Vector &x,
+         int print_iter = 0, int max_num_iter = 1000,
+         real_t RTOLERANCE = 1e-12, real_t ATOLERANCE = 1e-24,
+         real_t relax_factor = 1.0, FPIRelaxation *relax_method = nullptr);
+
 /// Conjugate gradient method
 class CGSolver : public IterativeSolver
 {
@@ -777,8 +777,31 @@ inline real_t InnerProduct(MPI_Comm comm, const Vector &x, const Vector &y)
   MPI_Allreduce(&loc_prod, &glb_prod, 1, MFEM_MPI_REAL_T, MPI_SUM, comm);
   return glb_prod;
 }
+
+/// Returns the square of the Euclidean distance between two vectors in parallel
+/** In parallel this computes the square of the Euclidean distance of the local
+    vectors, producing identical results on each MPI rank.*/
+inline real_t DistanceSquared(MPI_Comm comm, const real_t *x, const real_t *y,
+                              const int n)
+{
+   real_t d = DistanceSquared(x, y, n);
+   real_t glb_d;
+   MPI_Allreduce(&d, &glb_d, 1, MFEM_MPI_REAL_T, MPI_SUM, comm);
+   return glb_d;
+}
+
+inline real_t DistanceSquared(MPI_Comm comm, const Vector &x, const Vector &y)
+{
+   MFEM_ASSERT(x.Size() == y.Size(), "Incompatible vector sizes.");
+   real_t d = x.DistanceSquaredTo(y);
+   real_t glb_d;
+   MPI_Allreduce(&d, &glb_d, 1, MFEM_MPI_REAL_T, MPI_SUM, comm);
+   return glb_d;
+}
 #endif

+
+
 } // namespace mfem

 #endif
@@ -362,6 +362,22 @@ void NavierSolver::UpdateTimestepHistory(real_t dt)
   un_gf.SetFromTrueDofs(un);
 }

+void NavierSolver::Step(Vector &up, real_t &t, real_t &dt)
+{
+   Array<int> offsets({0, vfes->GetTrueVSize(), pfes->GetTrueVSize()});
+   offsets.PartialSum();
+   BlockVector upb(up.GetData(), offsets);
+
+   un_gf.SetFromTrueDofs(upb.GetBlock(0));
+   un_next_gf.SetFromTrueDofs(upb.GetBlock(0));
+   pn_gf.SetFromTrueDofs(upb.GetBlock(1));
+
+   Step(t, dt, (int) t/dt,false);
+
+   un_gf.GetTrueDofs(upb.GetBlock(0));
+   pn_gf.GetTrueDofs(upb.GetBlock(1));
+}
+
 void NavierSolver::Step(real_t &time, real_t dt, int current_step,
                        bool provisional)
 {
@@ -182,6 +182,8 @@ public:
    */
   void Step(real_t &time, real_t dt, int cur_step, bool provisional = false);

+   void Step(Vector &up, real_t &t, real_t &dt);
+
   /// Return a pointer to the provisional velocity ParGridFunction.
   ParGridFunction *GetProvisionalVelocity() { return &un_next_gf; }

@@ -0,0 +1,58 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+set(MESH_FILES
+  backward-facing-step.msh
+  channel-cylinder.msh
+)
+
+# Add a target to copy the mesh files from the source directory; used by sample
+# runs.
+set(SRC_MESH_FILES)
+foreach(MESH_FILE ${MESH_FILES})
+  list(APPEND SRC_MESH_FILES ${CMAKE_CURRENT_SOURCE_DIR}/${MESH_FILE})
+endforeach()
+add_custom_command(OUTPUT data_is_copied
+  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_FILES} .
+  COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
+  COMMENT "Copying multiapp miniapps data files ...")
+add_custom_target(copy_miniapps_multiapp_data DEPENDS data_is_copied)
+
+
+list(APPEND MULTIAPP_COMMON_SOURCES
+	multiapp.cpp)
+
+list(APPEND MULTIAPP_COMMON_HEADERS
+	multiapp.hpp)
+
+set(MULTIAPP_COMMON_FILES
+	EXTRA_SOURCES ${MULTIAPP_COMMON_SOURCES}
+	EXTRA_HEADERS ${MULTIAPP_COMMON_HEADERS})
+
+
+# Parallel apps.
+if (MFEM_USE_MPI)
+  add_mfem_miniapp(cht-BFS
+    MAIN cht-BFS.cpp
+    ${MFEM_MINIAPPS_COMMON_HEADERS} ${MULTIAPP_COMMON_FILES} ../navier/navier_solver.hpp ../navier/navier_solver.cpp
+    LIBRARIES mfem-common)
+  add_dependencies(cht-BFS copy_miniapps_multiapp_data)
+
+  add_mfem_miniapp(fsi
+    MAIN fsi.cpp
+    ${MFEM_MINIAPPS_COMMON_HEADERS} ${MULTIAPP_COMMON_FILES}
+    LIBRARIES mfem-common)
+  add_dependencies(fsi copy_miniapps_multiapp_data)
+
+# Add parallel tests.
+  # if (MFEM_ENABLE_TESTING)
+  # endif()
+endif()
@@ -0,0 +1,323 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+
+#ifndef MFEM_ELASTICITY_HPP
+#define MFEM_ELASTICITY_HPP
+
+#include "mfem.hpp"
+#include "../multiapp.hpp"
+
+
+using namespace mfem;
+
+/**
+ * @brief Elasticity time dependent operator
+ *
+ *              ρ•dv/dt = (-∇•σ + f)
+ *              du/dt = v
+ *
+ * where σ is the stress tensor, f is the body force,
+ * and ρ is the density.
+ */
+class Elasticity : public Application
+{
+public:
+
+   // Mesh and finite element space
+   ParMesh &mesh;
+   ParFiniteElementSpace &fes;
+   Array<int> offsets;
+
+   /// Essential and natural dof array.
+   Array<int> ess_attr, nat_attr;
+   Array<int> ess_tdofs, nat_tdofs;
+
+   /// Material properties
+   ConstantCoefficient density, mu, lambda;
+
+   /// Grid functions for the displacement, velocity, and traction
+   mutable ParGridFunction x_gf; ///< displacement
+   mutable ParGridFunction u_gf; ///< velocity (dx/dt)
+   mutable ParGridFunction stress_gf; ///< traction
+   mutable ParGridFunction bc_send_gf; ///< Grid functions for transfering BCs
+   
+   /// Mass and Stiffness forms
+   mutable ParBilinearForm Mform, Kform, Kform_e;
+
+   /// RHS form
+   mutable ParLinearForm bform;
+   mutable Vector b;
+
+   /// Mass and Stiffness operators
+   mutable HypreParMatrix Mmat, Kmat, Kmat_e, Mrhomat, Mmat_e;
+   mutable HypreParMatrix *T = nullptr;
+     
+   /// Force 
+   VectorGridFunctionCoefficient fcoeff;
+   ScalarVectorProductCoefficient scaled_fcoeff;
+
+   /// Mass matrix and implicit solver
+   GMRESSolver implicit_solver;   
+   CGSolver M_solver;
+
+   /// Preconditioner
+   HypreSmoother M_prec;
+   mutable HypreBoomerAMG *amg = nullptr;
+   mutable Solver *pc = nullptr;
+
+    /// Auxiliary variables
+   mutable Vector z;
+   real_t current_dt = -1.0;
+   bool updated = false;
+
+public:
+   Elasticity(ParFiniteElementSpace &fes_,
+              Array<int> ess_attr_,
+              Array<int> nat_attr_,
+              real_t density_ = 1.0,
+              real_t mu_ = 1.0,
+              real_t lambda_ = 1.0) :
+              Application(2*fes_.GetTrueVSize()),
+              mesh(*fes_.GetParMesh()),
+              fes(fes_),
+              ess_attr(ess_attr_),
+              nat_attr(nat_attr_),
+              density(density_),
+              mu(mu_),
+              lambda(lambda_),
+              x_gf(&fes), u_gf(&fes), stress_gf(&fes),
+              bc_send_gf(&fes),
+              Mform(&fes), Kform(&fes),
+              Kform_e(&fes), bform(&fes),
+              fcoeff(&stress_gf), scaled_fcoeff(-1.0, fcoeff),
+              implicit_solver(mesh.GetComm()),
+              M_solver(mesh.GetComm())
+   {
+        fes.GetEssentialTrueDofs(ess_attr, ess_tdofs);
+        fes.GetEssentialTrueDofs(nat_attr, nat_tdofs);
+
+        x_gf = 0.0;
+        u_gf = 0.0;
+        stress_gf = 0.0;
+        bc_send_gf = 0.0;
+
+        // Setup field collection for output and transfer
+        field_collection.SetName("Elasticity");
+        field_collection.AddField("Displacement", &x_gf);
+        field_collection.AddField("Velocity", &u_gf);
+        field_collection.AddField("Traction", &stress_gf);
+        field_collection.AddSourceField("Displacement_BC", &bc_send_gf);
+        field_collection.AddSourceField("Velocity_BC", &bc_send_gf);
+
+        offsets = Array<int>({0, fes.GetTrueVSize(), fes.GetTrueVSize()});
+        offsets.PartialSum();
+
+        Mform.AddDomainIntegrator(new VectorMassIntegrator(density));
+        Kform.AddDomainIntegrator(new ElasticityIntegrator(lambda, mu));
+        Kform_e.AddDomainIntegrator(new ElasticityIntegrator(lambda, mu));
+
+        if(nat_attr.Size() > 0)
+        {
+            bform.AddBoundaryIntegrator(new VectorBoundaryLFIntegrator(scaled_fcoeff), nat_attr);
+        }
+
+        b.SetSize(fes.GetTrueVSize());
+        z.SetSize(fes.GetTrueVSize());
+        Assemble();
+        BuildSolvers();
+   }
+
+   /// Assemble linear and bilinear forms; called if mesh is updated
+   void Assemble() override 
+   {
+        AssembleLinearForms();
+        AssembleBilinearForms();
+   }
+
+   /// Assemble linear forms for traction
+   void AssembleLinearForms()
+   {
+        bform.Assemble();
+        b.SetSize(fes.GetTrueVSize());
+        bform.ParallelAssemble(b);
+   }
+
+   /// Assemble bilinear forms for mass and stiffness
+   void AssembleBilinearForms()
+   {
+        Mform.Assemble();
+        Kform.Assemble();
+        Kform_e.Assemble();
+
+        Array<int> empty;
+        Mform.FormSystemMatrix(ess_tdofs, Mrhomat);
+        Kform_e.FormSystemMatrix(ess_tdofs, Kmat_e);
+        Kform.FormSystemMatrix(empty, Kmat);
+   }
+
+   /// Update finite element space and re-assemble forms
+   /// if the mesh has changed
+   void Update() override
+   {
+        fes.Update();
+
+        u_gf.Update();
+        x_gf.Update();
+        stress_gf.Update();
+        bc_send_gf.Update();
+    
+        Mform.Update();
+        Kform.Update();
+        Kform_e.Update();
+        bform.Update();
+
+        Assemble();
+        updated = true;
+   }
+
+   /// Build implicit solver and AMG preconditioner
+   void BuildSolvers()
+   {
+        M_solver.iterative_mode = false;
+        M_solver.SetRelTol(1e-8);
+        M_solver.SetAbsTol(1e-8);
+        M_solver.SetMaxIter(1000);
+        M_solver.SetPrintLevel(0);
+        M_solver.SetOperator(Mrhomat);
+        M_prec.SetType(HypreSmoother::Jacobi);
+        M_prec.SetOperator(Mrhomat);
+        M_solver.SetPreconditioner(M_prec);
+
+        if(amg) delete amg;
+        amg = new HypreBoomerAMG;
+        HYPRE_BoomerAMGSetSmoothType(*amg, 5);
+        amg->SetOperator(Kmat_e);
+        amg->SetSystemsOptions(2, true);
+        amg->SetElasticityOptions(&fes);
+        amg->SetPrintLevel(0);
+        pc = amg;
+
+        implicit_solver.iterative_mode = false;
+        implicit_solver.SetRelTol(1e-4);
+        implicit_solver.SetAbsTol(0.0);
+        implicit_solver.SetMaxIter(500);
+        implicit_solver.SetKDim(200);
+        implicit_solver.SetPrintLevel(0);
+        implicit_solver.SetPreconditioner(*pc);
+   }
+
+   /// Apply operator
+   void Mult(const Vector &u, Vector &k) const override
+   {
+        BlockVector ub(u.GetData(), offsets);
+        BlockVector kb(k.GetData(), offsets);
+
+        Vector &vel = ub.GetBlock(0);
+        Vector &pos = ub.GetBlock(1);
+
+        Vector &ku = kb.GetBlock(0);
+        Vector &kx = kb.GetBlock(1);
+        
+        Kmat.Mult(pos,z);
+        z.Neg();
+        z.Add(1.0, b);
+
+        kx = vel; // dx/dt = u
+        M_solver.Mult(z, ku);
+
+        ku.SetSubVector(ess_tdofs, 0.0);
+        kx.SetSubVector(ess_tdofs, 0.0);
+   }
+
+   /// Solve implicit system in Schur complement form
+   void ImplicitSolve(const real_t dt, const Vector &u, Vector &k)
+   {
+        BlockVector ub(u.GetData(), offsets);
+        BlockVector kb(k.GetData(), offsets);
+
+        Vector &vel = ub.GetBlock(0);
+        Vector &pos = ub.GetBlock(1);
+
+        Vector &ku = kb.GetBlock(0);
+        Vector &kx = kb.GetBlock(1);
+
+        AssembleLinearForms();
+        Kmat.Mult(pos,z);
+        z.Neg();
+        z.Add(1.0, b);
+        Kmat.AddMult(vel,z,-dt); // z = z  - dt*K*vel
+
+        if((current_dt != dt) || updated)
+        {
+            if (T) { delete T; }
+            current_dt = dt;
+            T = Add(1.0, Mrhomat, dt*dt, Kmat_e);
+            implicit_solver.SetOperator(*T);
+            amg->SetOperator(*T);
+            updated = false;
+        }
+
+        implicit_solver.Mult(z, ku);
+        add(1.0, vel, dt, ku, kx);
+
+        ku.SetSubVector(ess_tdofs, 0.0);
+        kx.SetSubVector(ess_tdofs, 0.0);
+   }
+
+   /// Computes residual of the elasticity equations
+   void ImplicitMult(const Vector &x, const Vector &k, Vector &v ) const override
+   {
+        BlockVector xb(x.GetData(), offsets);
+        BlockVector kb(k.GetData(), offsets);
+        BlockVector vb(v.GetData(), offsets);
+
+        // Apply the implicit operator
+        Mrhomat.Mult(kb.GetBlock(0), vb.GetBlock(0)); // v = A(x) + M*k
+        vb.GetBlock(0).Add(-1.0, b);
+        Kmat.AddMult(xb.GetBlock(1), vb.GetBlock(0), 1.0); // v = K*x
+
+        vb.GetBlock(1) = kb.GetBlock(1); // v = u
+        vb.GetBlock(1) -= xb.GetBlock(0); // v = k
+
+        vb.GetBlock(0).SetSubVector(ess_tdofs, 0.0);
+        vb.GetBlock(1).SetSubVector(ess_tdofs, 0.0);
+   }
+
+    void PreProcess(Vector &x) override {}
+
+    void PostProcess(Vector &x) override {}
+
+
+   void Transfer(const Vector &x) override
+   {
+        BlockVector xb(x.GetData(), offsets);
+        field_collection.Transfer("Velocity_BC", xb.GetBlock(0));
+        field_collection.Transfer("Displacement_BC", xb.GetBlock(1));
+   }
+
+   void Transfer(const Vector &u, const Vector &k, real_t dt = 0.0) override
+   {
+        BlockVector kb(k.GetData(), offsets);        
+        field_collection.Transfer("Velocity_BC", kb.GetBlock(0));
+        field_collection.Transfer("Displacement_BC", kb.GetBlock(1));
+   }
+
+   ~Elasticity() override
+   {
+        if (amg) { delete amg; }
+        if (T) { delete T; }
+   }
+};
+
+
+
+#endif
@@ -0,0 +1,235 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+
+#ifndef MFEM_MESH_MORPHER_HPP
+#define MFEM_MESH_MORPHER_HPP
+
+#include "mfem.hpp"
+#include "../multiapp.hpp"
+
+
+using namespace mfem;
+
+
+/**
+ * @brief Mesh morphing is modeled as a time dependent vector diffusion equation
+ *
+ *              dx/dt = κΔx
+ *
+ * with vector diffusion coefficient, κ.
+ */
+class MeshDiffusion : public Application
+{
+public: 
+
+   // Mesh and finite element space
+   ParMesh &mesh;
+   ParFiniteElementSpace &fes;
+
+   /// Essential true dof array.
+   Array<int> ess_attr, ess_tdofs;
+
+   /// Diffusion coefficient
+   ConstantCoefficient kappa;
+
+   /// Grid function for the mesh displacement variable
+   mutable ParGridFunction x_gf; 
+   mutable ParGridFunction x_gf_bc;
+   mutable ParGridFunction u_gf_bc;
+   mutable ParGridFunction bc_send_gf; ///< Grid functions for transfering BCs
+
+   /// Mass and Stiffness forms
+   mutable ParBilinearForm Mform, Kform,  Kform_e;
+
+   /// Mass and Stiffness operators
+   OperatorHandle M, K;
+   mutable HypreParMatrix Mmat, Kmat,  Mmat_e, Kmat_e;
+
+   /// Mass matrix and implicit solver
+   mutable CGSolver M_solver;
+   mutable GMRESSolver implicit_solver;
+   HypreParMatrix *T = nullptr; // T = M + dt K
+
+   /// Mass matrix preconditioner
+   HypreSmoother M_prec;
+   HypreSmoother T_prec;
+
+   real_t current_dt = -1.0; 
+
+   /// Auxiliary vectors
+   mutable Vector z, zv;
+   bool updated = false;
+
+public:
+
+     MeshDiffusion(ParFiniteElementSpace &fes_,
+                   Array<int> ess_attr_,
+                   real_t kappa_ = 1.0e0) : 
+                   Application(fes_.GetTrueVSize()),
+                   mesh(*fes_.GetParMesh()),
+                   fes(fes_),
+                   ess_attr(ess_attr_),
+                   kappa(kappa_),
+                   x_gf(&fes),
+                   x_gf_bc(&fes), u_gf_bc(&fes),
+                   bc_send_gf(&fes),
+                   Mform(&fes), Kform(&fes),Kform_e(&fes),
+                   M_solver(mesh.GetComm()),
+                   implicit_solver(mesh.GetComm()),
+                   z(fes_.GetTrueVSize())
+     {
+          fes.GetEssentialTrueDofs(ess_attr, ess_tdofs);
+
+          x_gf = 0.0;
+          x_gf_bc = 0.0;
+          u_gf_bc = 0.0;
+          bc_send_gf = 0.0;
+
+          // Setup field collection for output and transfer
+          field_collection.SetName("Mesh-Diffusion");
+          field_collection.AddSourceField("Displacement", &x_gf);
+          field_collection.AddSourceField("dxdt", &bc_send_gf);
+          field_collection.AddField("Displacement_BC", &x_gf_bc);
+          field_collection.AddField("Velocity_BC", &u_gf_bc);
+
+          Mform.AddDomainIntegrator(new VectorMassIntegrator);
+          Kform.AddDomainIntegrator(new VectorDiffusionIntegrator(kappa));
+          Kform_e.AddDomainIntegrator(new VectorDiffusionIntegrator(kappa));
+
+          Assemble();
+          BuildSolvers();
+     }
+
+     void Assemble() override 
+     {
+          Mform.Assemble();
+          Kform.Assemble();
+          Kform_e.Assemble();
+
+          Array<int> empty;
+          Mform.FormSystemMatrix(ess_tdofs, Mmat_e);
+          Kform_e.FormSystemMatrix(ess_tdofs, Kmat_e);
+          Kform.FormSystemMatrix(empty, Kmat);
+     }
+
+     void Update() override
+     {
+          fes.Update();
+          x_gf.Update();
+          x_gf_bc.Update();
+          u_gf_bc.Update();
+          bc_send_gf.Update();
+
+          Mform.Update();
+          Kform.Update();
+          Kform_e.Update();
+
+          Assemble();
+          updated = true;
+     }
+
+     void BuildSolvers()
+     {
+          M_solver.iterative_mode = false;
+          M_solver.SetRelTol(1e-8);
+          M_solver.SetAbsTol(0.0);
+          M_solver.SetMaxIter(100);
+          M_solver.SetPrintLevel(0);
+          M_solver.SetOperator(Mmat_e);
+          M_prec.SetType(HypreSmoother::Jacobi);
+          M_solver.SetPreconditioner(M_prec);
+
+          implicit_solver.iterative_mode = false;
+          implicit_solver.SetRelTol(1e-8);
+          implicit_solver.SetAbsTol(1e-5);
+          implicit_solver.SetMaxIter(500);
+          implicit_solver.SetPrintLevel(0);
+          T_prec.SetType(HypreSmoother::Jacobi);
+          implicit_solver.SetPreconditioner(T_prec);
+     }
+
+     void Mult(const Vector &x, Vector &k) const override
+     {
+          Kmat.Mult(x, z);
+          z.Neg();
+          M_solver.Mult(z, k);
+          k.SetSubVector(ess_tdofs, 0.0);
+     }
+
+     void ImplicitSolve(const real_t dt, const Vector &x, Vector &k)
+     {
+          if((current_dt != dt) || updated)
+          {
+               if (T) delete T;
+               T = Add(1.0, Mmat_e, dt, Kmat_e);
+               implicit_solver.SetOperator(*T);
+               current_dt = dt;
+               updated = false;
+          }
+
+          Kmat.Mult(x, z);
+          z.Neg();
+          implicit_solver.Mult(z, k);
+
+          if(IsCoupled())
+          {
+               x_gf_bc.GetTrueDofs(z); // contains the BC in terms of kx (dx/dt)
+               for (int i = 0; i < ess_tdofs.Size(); i++)
+               {
+                    int idx = ess_tdofs[i];
+                    k(idx) = z(idx);
+               }
+          }
+          else
+          {
+               k.SetSubVector(ess_tdofs, 0.0);
+          }
+     }
+
+     void ImplicitMult(const Vector &x, const Vector &k, Vector &v ) const override
+     {}
+
+
+     void PreProcess(Vector &x) override
+     {}
+
+     void PostProcess(Vector &x) override {}
+
+     void Transfer(const Vector &x) override
+     {
+          field_collection.Transfer("Displacement", x);
+
+          zv.SetSize(fes.GetTrueVSize());
+          ParGridFunction *dxdt = field_collection.GetField("Velocity_BC");
+          dxdt->GetTrueDofs(zv);
+          bc_send_gf.GetTrueDofs(z);
+          for (int i = 0; i < ess_tdofs.Size(); i++)
+          {
+               int idx = ess_tdofs[i];
+               z(idx) = zv(idx);
+          }
+          field_collection.Transfer("dxdt", z);
+     }
+
+     void Transfer(const Vector &u, const Vector &k, real_t dt = 0.0) override
+     {
+          field_collection.Transfer("dxdt", k);
+     }
+
+     ~MeshDiffusion() override
+     {
+          if (T) delete T;
+     }
+};
+
+
+#endif
@@ -0,0 +1,503 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+
+#ifndef MFEM_NAVIER_STOKES_HPP
+#define MFEM_NAVIER_STOKES_HPP
+
+#include "mfem.hpp"
+#include "../multiapp.hpp"
+
+
+using namespace mfem;
+
+/// @brief Compute deviatoric stress on a boundary
+class DeviatoricStressCoefficient : public VectorCoefficient
+{
+protected:
+    ParGridFunction *p_gf  = nullptr;  ///< Grid function for pressure
+    ParGridFunction *u_gf  = nullptr;  ///< Grid function for velocity
+    Coefficient *viscosity = nullptr;  ///< Coefficient for the kinematic viscosity
+    Coefficient *density   = nullptr;  ///< Coefficient for the density
+    bool scaled_pressure   = false;    ///< True if pressure is density-scaled
+
+    DenseMatrix dudx, tau; ///< Velocity gradient tensor
+    Vector normal;
+
+public:
+    // For velocity and pressure case
+    DeviatoricStressCoefficient(ParGridFunction *p_gf_,
+                                ParGridFunction *u_gf_,
+                                Coefficient *viscosity_,
+                                Coefficient *density_,
+                                bool scaled_pressure_ = false) :
+                                VectorCoefficient(p_gf_->FESpace()->GetMesh()->Dimension()),
+                                p_gf(p_gf_), u_gf(u_gf_),
+                                viscosity(viscosity_), density(density_), scaled_pressure(scaled_pressure_)
+    {
+        int ndim = p_gf->FESpace()->GetMesh()->Dimension();
+        normal.SetSize(ndim);
+        if(u_gf)
+        {
+            dudx.SetSize(ndim, vdim);
+            tau.SetSize(ndim, vdim);
+        }
+    }
+
+    // For pressure-only case
+    DeviatoricStressCoefficient(ParGridFunction *p_gf_) :
+                                DeviatoricStressCoefficient(p_gf_,NULL,NULL,NULL) {}
+
+    /// Evaluate the deviatoric stress tensor at the given point.
+    void Eval(Vector &v, ElementTransformation &T,
+              const IntegrationPoint &ip) override
+    {
+        v.SetSize(vdim);
+        v=0.0;
+
+        const DenseMatrix &jacobian = T.Jacobian();
+
+        MFEM_ASSERT( (jacobian.Height() - 1 == jacobian.Width()), 
+                     "Incorrect Jacobian dimension. Coefficient only supported for boundary elements.");
+
+        CalcOrtho(jacobian, normal);
+        const double scale = normal.Norml2();
+        normal /= scale;
+
+        real_t rho = density ? density->Eval(T, ip) : 1.0;
+
+        if(p_gf)
+        {   // Add pressure term to the deviatoric stress tensor
+            real_t pscale = scaled_pressure ? rho : 1.0;
+            for (int i = 0; i < vdim; ++i)
+            {
+                v(i) = -pscale*p_gf->GetValue(T, ip) * normal(i);
+            }
+        }
+
+        if(u_gf)
+        {
+            real_t nu  = viscosity ? viscosity->Eval(T, ip) : 1.0;
+            u_gf->GetVectorGradient(T, dudx);
+
+            for (int i = 0; i < vdim; i++)
+            {
+                for (int j = 0; j < vdim; j++)
+                {
+                    tau(i,j) = rho*nu*(dudx(i,j) + dudx(j,i));
+                }
+            }
+            tau.AddMult(normal, v);
+        }
+    }
+};
+
+/**
+ * @brief Navier-stokes time dependent operator
+ *  du/dt = - div(u u) + nu * laplacian(u) - grad(p) + f
+ *       div(u) = 0
+ * 
+ *   M du/dt = [ -N(u) + L, -G][u]
+ *         0 = [ D        ,  0][p]
+ */
+class NavierStokes : public Application
+{
+protected:
+
+  class NewtonResidual : public TimeDependentOperator
+  {
+  protected:
+    NavierStokes *app;
+    const Vector *u = nullptr; ///< Pointer to the input vector, used in Mult() methods
+    mutable future::FDJacobian grad;
+    real_t dt= 0.0;     ///< Time step size, used for time-dependent applications
+    mutable Vector upk; ///< Temporary vector for u + dt*k
+
+  public:
+    NewtonResidual(NavierStokes *app_) : TimeDependentOperator(app_->Width()), 
+                                         app(app_), grad(*this,1e-6), upk(Width()) {}
+    virtual void SetTimeStep(real_t dt_){ dt = dt_;}
+    virtual void SetState(const Vector *u_){u = u_;}
+    virtual void Mult(const Vector &k, Vector &y) const override
+    {
+        add(1.0,*u, dt, k, upk); // upk = u + dt*k
+        app->ImplicitMult(upk, k, y); // y = f(upk,k,t)
+    }
+    Operator& GetGradient(const Vector &k) const override {
+        // grad.Update(k);
+        // return const_cast<future::FDJacobian&>(grad);
+        add(1.0,*u, dt, k, upk); // upk = u + dt*k 
+        return app->GetGradient(upk);
+    }
+};
+
+public:
+
+   // Mesh and finite element spaces
+   ParMesh &mesh;
+   ParFiniteElementSpace &u_fes, &p_fes;
+   Array<int> offsets; ///< Offsets for the velocity and pressure spaces
+
+   /// Essential true dof array. Relevant for eliminating boundary conditions
+   Array<int> u_ess_attr, p_ess_attr;
+   Array<int> u_ess_tdofs, p_ess_tdofs;
+
+   ConstantCoefficient viscosity, density, inv_density, compressibility; // Fluid properties
+
+   mutable ParGridFunction u_gf, p_gf; ///< Grid function for the velocity and pressure
+   mutable ParGridFunction u_gf_bc, p_gf_bc; ///< Grid functions for enforcing/transfering BCs
+
+   mutable ParBilinearForm Mpform, Muform, Kuform; ///< Mass and Stiffness forms
+   mutable ParMixedBilinearForm Gform, Dform;   ///< Divergence and Gradient forms
+   mutable ParNonlinearForm Nform; ///< Nonlinear form for the Navier-Stokes equations
+   mutable ParNonlinearForm Nform_e; ///< Gradient for nonlinear form for the Navier-Stokes equations
+
+   /// Mass matrix and implicit solver
+   CGSolver M_solver;
+   FGMRESSolver linear_solver;
+   NewtonSolver newton_solver;
+
+   mutable OperatorHandle Mumat, Mpmat;
+   mutable OperatorHandle Kumat, Dmat, Gmat;
+   mutable OperatorHandle Kumat_e, Dmat_e, Gmat_e, Mpmat_e;
+
+   BlockOperator *implicit_op   = nullptr; ///< Block operator for the Navier-Stokes system
+   BlockOperator *implicit_grad = nullptr; ///< Block linear operator for the Navier-Stokes gradient
+   NewtonResidual *newton_residual = nullptr; ///< Newton residual operator
+
+   /// Mass matrix preconditioner
+   HypreSmoother M_prec;
+   mutable BlockLowerTriangularPreconditioner *pc = nullptr;
+   mutable Solver *Nu_pc = nullptr;
+   mutable Vector z, pz, uz;
+
+   ConstantCoefficient zero, one; ///< Zero coefficient for boundary conditions
+   ConstantCoefficient inv_dtc;
+
+   IntegrationRules intrules;
+   IntegrationRule ir, ir_nl;
+
+   DeviatoricStressCoefficient stress_coeff;
+   ParGridFunction stress_gf;
+   VectorCoefficient *ale_velocity; ///< ALE velocity for moving mesh problems
+
+   real_t current_dt = -1.0; ///< Current time in the simulation
+   bool updated = false;
+
+public:
+   NavierStokes(ParFiniteElementSpace &u_fes_,
+                ParFiniteElementSpace &p_fes_,
+                Array<int> u_ess_attr,
+                Array<int> p_ess_attr,
+                real_t density_ = 1.0,
+                real_t viscosity_ = 1.0,
+                real_t compressibility_ = 1.0,
+                bool scaled_pressure = true,
+                VectorCoefficient *ale_velocity_ = nullptr) :
+                Application(u_fes_.GetTrueVSize()+p_fes_.GetTrueVSize()),
+                mesh(*p_fes_.GetParMesh()),
+                u_fes(u_fes_), p_fes(p_fes_),
+                u_ess_attr(u_ess_attr), p_ess_attr(p_ess_attr),
+                viscosity(viscosity_), density(density_),
+                inv_density(scaled_pressure ? 1.0 : 1.0/density_),
+                compressibility(compressibility_),
+                u_gf(&u_fes), p_gf(&p_fes),
+                u_gf_bc(&u_fes), p_gf_bc(&p_fes),
+                Mpform(&p_fes), Muform(&u_fes),
+                Kuform(&u_fes), Gform(&p_fes, &u_fes),
+                Dform(&u_fes, &p_fes),
+                Nform(&u_fes), Nform_e(&u_fes),
+                linear_solver(mesh.GetComm()),
+                newton_solver(mesh.GetComm()),
+                zero(0.0), one(1.0), inv_dtc(1.0),
+                stress_coeff(&p_gf, &u_gf, &viscosity, &density, scaled_pressure),
+                stress_gf(&u_fes), ale_velocity(ale_velocity_)
+   {
+        u_fes.GetEssentialTrueDofs(u_ess_attr, u_ess_tdofs);
+        p_fes.GetEssentialTrueDofs(p_ess_attr, p_ess_tdofs);
+
+        u_gf = 0.0;
+        p_gf = 0.0;
+        stress_gf = 0.0;
+        u_gf_bc   = 0.0;
+        p_gf_bc   = 0.0;
+
+        // Setup field collection for output and transfer
+        field_collection.SetName("Navier-Stokes");
+        field_collection.AddField("Velocity", &u_gf);
+        field_collection.AddField("Pressure", &p_gf);
+        field_collection.AddSourceField("Stress", &stress_gf);
+        field_collection.AddField("Velocity_BC", &u_gf_bc);
+        field_collection.AddField("Pressure_BC", &p_gf_bc);
+
+        offsets = Array<int>({0, u_fes.GetTrueVSize(), p_fes.GetTrueVSize()});
+        offsets.PartialSum();
+
+        auto geom_type = u_fes.GetFE(0)->GetGeomType();
+        ir    = intrules.Get(geom_type, (int)(2*(u_fes.GetOrder(0)+1) - 3));
+        ir_nl = intrules.Get(geom_type, (int)(ceil(1.5 * 2*(u_fes.GetOrder(0)+1) - 3)));
+
+        Mpform.AddDomainIntegrator(new MassIntegrator(compressibility,&ir));
+        Muform.AddDomainIntegrator(new VectorMassIntegrator(&ir));
+        // Kuform.AddDomainIntegrator(new VectorDiffusionIntegrator(viscosity, &ir));
+        Gform.AddDomainIntegrator(new GradientIntegrator(inv_density, &ir));
+        Dform.AddDomainIntegrator(new VectorDivergenceIntegrator(&ir));
+        Nform.AddDomainIntegrator(new VectorConvectionNLFIntegrator(one, *ale_velocity, &ir_nl));
+        Nform.AddDomainIntegrator(new VectorDiffusionIntegrator(viscosity, &ir));
+
+        Nform_e.AddDomainIntegrator(new VectorConvectionNLFIntegrator(one, *ale_velocity, &ir_nl));
+        Nform_e.AddDomainIntegrator(new VectorDiffusionIntegrator(viscosity, &ir));
+        Nform_e.AddDomainIntegrator(new VectorMassIntegrator(inv_dtc,&ir));
+
+        // Muform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+        Mpform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+        Gform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+        Dform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+        Nform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+
+        Assemble();
+        BuildSolvers();
+   }
+
+   void Assemble() override
+   {
+        AssembleLinearForms();
+        AssembleNonlinearForm();
+        AssembleBilinearForms();
+   }
+
+   void AssembleLinearForms() {}
+
+   void AssembleBilinearForms()
+   {
+        Mpform.Assemble();
+        Muform.Assemble();
+        // Kuform.Assemble();
+        Dform.Assemble();
+        Gform.Assemble();
+
+        // Muform.FormSystemMatrix(u_ess_tdofs, Mumat);
+        Mpform.FormSystemMatrix(p_ess_tdofs, Mpmat_e);
+        Gform.FormRectangularSystemMatrix(p_ess_tdofs, u_ess_tdofs, Gmat_e);
+        Dform.FormRectangularSystemMatrix(u_ess_tdofs, p_ess_tdofs, Dmat_e);
+
+        Array<int> empty;
+        Mpform.FormSystemMatrix(empty, Mpmat);
+        Muform.FormSystemMatrix(empty, Mumat);
+        Gform.FormRectangularSystemMatrix(empty, empty, Gmat);
+        Dform.FormRectangularSystemMatrix(empty, empty, Dmat);
+
+        if(!implicit_op) implicit_op = new BlockOperator(offsets);
+        implicit_op->SetBlock(0, 0, &Nform);
+        implicit_op->SetBlock(0, 1, Gmat.Ptr());
+        implicit_op->SetBlock(1, 0, Dmat.Ptr());
+
+        if(!implicit_grad) implicit_grad = new BlockOperator(offsets);
+        implicit_grad->SetBlock(0, 1, Gmat_e.Ptr());
+        implicit_grad->SetBlock(1, 0, Dmat_e.Ptr());
+        implicit_grad->SetBlock(1, 1, Mpmat_e.Ptr());
+   }
+
+   void AssembleNonlinearForm()
+   {
+        Nform.Setup();
+        Nform_e.SetEssentialTrueDofs(u_ess_tdofs);
+        Nform_e.Setup();
+   }
+
+   void Update() override
+   {
+        u_fes.Update();
+        p_fes.Update();
+
+        u_gf.Update();
+        p_gf.Update();
+        stress_gf.Update();
+        p_gf_bc.Update();
+        u_gf_bc.Update();
+
+        Mpform.Update();
+        Muform.Update();
+        // Kuform.Update();
+        Dform.Update();
+        Gform.Update();
+        Nform.Update();
+        Nform_e.Update();
+
+        Assemble();
+        updated = true;
+   }
+
+   void BuildSolvers()
+   {
+        M_solver.iterative_mode = false;
+        M_solver.SetRelTol(1e-8);
+        M_solver.SetMaxIter(100);
+        M_prec.SetType(HypreSmoother::Jacobi);
+        M_solver.SetPreconditioner(M_prec);
+        M_solver.SetOperator(*Mumat.Ptr());
+
+        linear_solver.iterative_mode = false;
+        linear_solver.SetRelTol(1e-4);
+        linear_solver.SetAbsTol(1e-4);
+        linear_solver.SetMaxIter(500);
+        linear_solver.SetKDim(300);
+        linear_solver.SetPrintLevel(0);
+        // pc = new BlockLowerTriangularPreconditioner(offsets);
+        // pc->SetBlock(1, 0, Dmat_e.Ptr());
+        // linear_solver.SetPreconditioner(*pc);
+
+        if(newton_residual) delete newton_residual;
+        newton_residual = new NewtonResidual(this);
+        newton_solver.iterative_mode = true;
+        newton_solver.SetRelTol(0.0);
+        newton_solver.SetAbsTol(1e-4);
+        newton_solver.SetMaxIter(30);
+        newton_solver.SetPrintLevel(0);
+        newton_solver.SetOperator(*newton_residual);
+        newton_solver.SetPreconditioner(linear_solver);
+   }
+
+   void UpdatePreconditioner() const
+   {
+      if(Nu_pc) delete Nu_pc;
+
+      auto amg = new HypreBoomerAMG;
+      HYPRE_BoomerAMGSetSmoothType(*amg, 5);
+      amg->SetOperator(*static_cast<HypreParMatrix*>(&implicit_grad->GetBlock(0, 0)));
+      amg->SetSystemsOptions(2, true);
+      amg->SetPrintLevel(0);
+      Nu_pc = amg;
+      pc->SetBlock(0, 0, Nu_pc);
+   }
+
+   Operator& GetGradient(const Vector &x) const override
+   {
+        BlockVector xb(x.GetData(), offsets);
+        implicit_grad->SetBlock(0, 0, &Nform_e.GetGradient(xb.GetBlock(0)));
+        implicit_grad->SetBlockCoef(0,0,current_dt);
+        implicit_grad->SetBlockCoef(0,1,current_dt);
+        implicit_grad->SetBlockCoef(1,0,current_dt);
+
+        if(pc) UpdatePreconditioner();
+
+        return *implicit_grad;
+   }
+
+   void Mult(const Vector &x, Vector &k) const override
+   {
+        BlockVector xb(x.GetData(), offsets);
+        BlockVector kb(k.GetData(), offsets);
+
+        implicit_op->Mult(xb, kb);
+
+        M_solver.Mult(kb.GetBlock(0), z);
+        kb.GetBlock(0) = z;
+
+        kb.GetBlock(0).SetSubVector(u_ess_tdofs, 0.0);
+        kb.GetBlock(1).SetSubVector(p_ess_tdofs, 0.0);
+   }
+
+   void ImplicitSolve(const real_t dt, const Vector &u, Vector &k) override
+   {
+        BlockVector ub(u.GetData(), offsets);
+        BlockVector kb(k.GetData(), offsets);
+
+        if((current_dt != dt) || updated)
+        {
+            inv_dtc.constant = 1.0/dt;
+            current_dt = dt;
+            AssembleNonlinearForm();
+            updated = false;
+        }
+
+        if(IsCoupled())
+        {   // Enforce velocity BCs through initial guess
+            Vector &ku = kb.GetBlock(0);
+            u_gf_bc.GetTrueDofs(z);
+            for (int i = 0; i < u_ess_tdofs.Size(); i++)
+            {
+                int idx = u_ess_tdofs[i];
+                ku(idx) = z(idx);
+            }
+        }
+        else
+        {
+            kb.GetBlock(0) = 0.0;
+            kb.GetBlock(1) = 0.0;
+        }
+
+        Vector zero_vec;
+        newton_residual->SetTimeStep(dt);
+        newton_residual->SetState(&u);
+        newton_solver.Mult(zero_vec, k); // Solve the nonlinear system
+   }
+
+   // Computes residual of the Navier-Stokes equations
+   void ImplicitMult(const Vector &x, const Vector &k, Vector &v ) const override
+   {
+        BlockVector xb(x.GetData(), offsets);
+        BlockVector kb(k.GetData(), offsets);
+        BlockVector vb(v.GetData(), offsets);
+
+        // Apply the implicit operator
+        implicit_op->Mult(xb, vb); // v = A(x)
+        Mumat->AddMult(kb.GetBlock(0), vb.GetBlock(0), 1.0); // v = A(x) + M*k
+        Mpmat->AddMult(kb.GetBlock(1), vb.GetBlock(1), 1.0); // v = A(x) + M*k
+
+        vb.GetBlock(0).SetSubVector(u_ess_tdofs, 0.0);
+        vb.GetBlock(1).SetSubVector(p_ess_tdofs, 0.0);
+   }
+
+    void PreProcess(Vector &x) override {}
+
+    void PostProcess(Vector &x) override {}
+
+    void Transfer(const Vector &x) override
+    {
+        BlockVector xb(x.GetData(), offsets);
+        u_gf.SetFromTrueDofs(xb.GetBlock(0));
+        p_gf.SetFromTrueDofs(xb.GetBlock(1));
+        stress_gf.ProjectBdrCoefficient(stress_coeff,u_ess_attr);
+        field_collection.Transfer("Stress");
+    }
+
+    void Transfer(const Vector &u, const Vector &k, real_t dt = 0.0) override
+    {
+        BlockVector ub(u.GetData(), offsets);
+        BlockVector kb(k.GetData(), offsets);
+
+        uz.SetSize(u_fes.GetTrueVSize());
+        pz.SetSize(p_fes.GetTrueVSize());
+
+        // compute stage-updated velocity and pressure
+        add(1.0,ub.GetBlock(0), dt, kb.GetBlock(0), uz);
+        add(1.0,ub.GetBlock(1), dt, kb.GetBlock(1), pz);
+
+        u_gf.SetFromTrueDofs(uz);
+        p_gf.SetFromTrueDofs(pz);
+        stress_gf.ProjectBdrCoefficient(stress_coeff,u_ess_attr);
+        field_collection.Transfer("Stress");
+    }
+
+   ~NavierStokes()
+   {
+        if(implicit_op) delete implicit_op;
+        if(implicit_grad) delete implicit_grad;
+        if(newton_residual) delete newton_residual;
+        if(Nu_pc) delete Nu_pc;
+        if(pc) delete pc;
+    }
+};
+
+
+
+#endif
@@ -0,0 +1,57 @@
+SetFactory("OpenCASCADE");
+Mesh.ElementOrder = 2;
+Mesh.SecondOrderLinear = 0;
+
+lc = 1.0;
+
+// channel height
+H = 1.0;
+// step height
+h = H / 2.0;
+// slab thickness
+b = 4.0 * h;
+// channel length
+L = 30.0 * H;
+
+Point(1) = {0, 0, 0, lc};
+Point(2) = {L, 0, 0, lc};
+Point(3) = {L, b, 0, lc};
+Point(4) = {L, b+h, 0, lc};
+Point(5) = {L, b+H, 0, lc};
+Point(6) = {0, b+H, 0, lc};
+Point(7) = {0, b+h, 0, lc};
+Point(8) = {0, b, 0, lc};
+
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 5};
+Line(5) = {5, 6};
+Line(6) = {6, 7};
+Line(7) = {7, 8};
+Line(8) = {8, 1};
+Line(9) = {8, 3};
+Line(10) = {7, 4};
+
+Curve Loop(1) = {10, 4, 5, 6};
+Plane Surface(1) = {1};
+Curve Loop(2) = {9, 3, -10, 7};
+Plane Surface(2) = {2};
+Curve Loop(3) = {9, -2, -1, -8};
+Plane Surface(3) = {3};
+
+Physical Surface("fluid", 1) = {1, 2};
+Physical Surface("solid", 2) = {3};
+Physical Curve("inlet", 1) = {6};
+Physical Curve("outlet", 2) = {3, 4};
+Physical Curve("interface", 3) = {9};
+Physical Curve("wall fluid", 4) = {7, 5};
+Physical Curve("wall side", 5) = {8, 2};
+Physical Curve("wall bottom", 6) = {1};
+
+Transfinite Surface {1:3};
+Recombine Surface {1:3};
+
+Transfinite Curve {1, -5, 10, 9} = 8 Using Progression 1;
+Transfinite Curve {6, 4, 7, 3} = 1 Using Progression 1;
+Transfinite Curve {8, 2} = 2 Using Progression 1;
@@ -0,0 +1,173 @@
+$MeshFormat
+2.2 0 8
+$EndMeshFormat
+$PhysicalNames
+8
+1 1 "inlet"
+1 2 "outlet"
+1 3 "interface"
+1 4 "wall fluid"
+1 5 "wall side"
+1 6 "wall bottom"
+2 1 "fluid"
+2 2 "solid"
+$EndPhysicalNames
+$Nodes
+105
+1 0 0 0
+2 30 0 0
+3 30 2 0
+4 30 2.5 0
+5 30 3 0
+6 0 3 0
+7 0 2.5 0
+8 0 2 0
+9 4.285714285714286 0 0
+10 8.571428571428571 0 0
+11 12.85714285714286 0 0
+12 17.14285714285714 0 0
+13 21.42857142857143 0 0
+14 25.71428571428571 0 0
+15 2.142857142857143 0 0
+16 6.428571428571429 0 0
+17 10.71428571428571 0 0
+18 15 0 0
+19 19.28571428571428 0 0
+20 23.57142857142857 0 0
+21 27.85714285714285 0 0
+22 30 0.9999999999999999 0
+23 30 2.25 0
+24 30 2.75 0
+25 25.71428571428572 3 0
+26 21.42857142857143 3 0
+27 17.14285714285715 3 0
+28 12.85714285714286 3 0
+29 8.571428571428573 3 0
+30 4.285714285714288 3 0
+31 27.85714285714286 3 0
+32 23.57142857142857 3 0
+33 19.28571428571428 3 0
+34 15 3 0
+35 10.71428571428572 3 0
+36 6.428571428571431 3 0
+37 2.142857142857146 3 0
+38 0 2.75 0
+39 0 2.25 0
+40 0 1 0
+41 4.285714285714286 2 0
+42 8.571428571428571 2 0
+43 12.85714285714286 2 0
+44 17.14285714285714 2 0
+45 21.42857142857143 2 0
+46 25.71428571428571 2 0
+47 2.142857142857143 2 0
+48 6.428571428571429 2 0
+49 10.71428571428571 2 0
+50 15 2 0
+51 19.28571428571428 2 0
+52 23.57142857142857 2 0
+53 27.85714285714285 2 0
+54 4.285714285714286 2.5 0
+55 8.571428571428571 2.5 0
+56 12.85714285714286 2.5 0
+57 17.14285714285714 2.5 0
+58 21.42857142857143 2.5 0
+59 25.71428571428571 2.5 0
+60 2.142857142857143 2.5 0
+61 6.428571428571429 2.5 0
+62 10.71428571428571 2.5 0
+63 15 2.5 0
+64 19.28571428571428 2.5 0
+65 23.57142857142857 2.5 0
+66 27.85714285714285 2.5 0
+67 4.285714285714286 2.75 0
+68 2.142857142857146 2.75 0
+69 8.571428571428573 2.75 0
+70 6.428571428571431 2.75 0
+71 12.85714285714286 2.75 0
+72 10.71428571428572 2.75 0
+73 17.14285714285714 2.75 0
+74 15 2.75 0
+75 21.42857142857143 2.75 0
+76 19.28571428571428 2.75 0
+77 25.71428571428572 2.75 0
+78 23.57142857142857 2.75 0
+79 27.85714285714286 2.75 0
+80 4.285714285714285 2.25 0
+81 2.142857142857142 2.25 0
+82 8.571428571428571 2.25 0
+83 6.428571428571431 2.25 0
+84 12.85714285714286 2.25 0
+85 10.71428571428571 2.25 0
+86 17.14285714285714 2.25 0
+87 15 2.25 0
+88 21.42857142857143 2.25 0
+89 19.28571428571428 2.25 0
+90 25.71428571428571 2.25 0
+91 23.57142857142857 2.25 0
+92 27.85714285714285 2.25 0
+93 4.285714285714285 1 0
+94 2.142857142857142 1 0
+95 8.571428571428571 1 0
+96 6.428571428571431 1 0
+97 12.85714285714286 1 0
+98 10.71428571428571 1 0
+99 17.14285714285714 1 0
+100 15 1 0
+101 21.42857142857143 1 0
+102 19.28571428571428 1 0
+103 25.71428571428571 1 0
+104 23.57142857142857 1 0
+105 27.85714285714285 1 0
+$EndNodes
+$Elements
+48
+1 8 2 6 1 1 9 15
+2 8 2 6 1 9 10 16
+3 8 2 6 1 10 11 17
+4 8 2 6 1 11 12 18
+5 8 2 6 1 12 13 19
+6 8 2 6 1 13 14 20
+7 8 2 6 1 14 2 21
+8 8 2 5 2 2 3 22
+9 8 2 2 3 3 4 23
+10 8 2 2 4 4 5 24
+11 8 2 4 5 5 25 31
+12 8 2 4 5 25 26 32
+13 8 2 4 5 26 27 33
+14 8 2 4 5 27 28 34
+15 8 2 4 5 28 29 35
+16 8 2 4 5 29 30 36
+17 8 2 4 5 30 6 37
+18 8 2 1 6 6 7 38
+19 8 2 4 7 7 8 39
+20 8 2 5 8 8 1 40
+21 8 2 3 9 8 41 47
+22 8 2 3 9 41 42 48
+23 8 2 3 9 42 43 49
+24 8 2 3 9 43 44 50
+25 8 2 3 9 44 45 51
+26 8 2 3 9 45 46 52
+27 8 2 3 9 46 3 53
+28 10 2 1 1 7 54 30 6 60 67 37 38 68
+29 10 2 1 1 54 55 29 30 61 69 36 67 70
+30 10 2 1 1 55 56 28 29 62 71 35 69 72
+31 10 2 1 1 56 57 27 28 63 73 34 71 74
+32 10 2 1 1 57 58 26 27 64 75 33 73 76
+33 10 2 1 1 58 59 25 26 65 77 32 75 78
+34 10 2 1 1 59 4 5 25 66 24 31 77 79
+35 10 2 1 2 8 41 54 7 47 80 60 39 81
+36 10 2 1 2 41 42 55 54 48 82 61 80 83
+37 10 2 1 2 42 43 56 55 49 84 62 82 85
+38 10 2 1 2 43 44 57 56 50 86 63 84 87
+39 10 2 1 2 44 45 58 57 51 88 64 86 89
+40 10 2 1 2 45 46 59 58 52 90 65 88 91
+41 10 2 1 2 46 3 4 59 53 23 66 90 92
+42 10 2 2 3 8 41 9 1 47 93 15 40 94
+43 10 2 2 3 41 42 10 9 48 95 16 93 96
+44 10 2 2 3 42 43 11 10 49 97 17 95 98
+45 10 2 2 3 43 44 12 11 50 99 18 97 100
+46 10 2 2 3 44 45 13 12 51 101 19 99 102
+47 10 2 2 3 45 46 14 13 52 103 20 101 104
+48 10 2 2 3 46 3 2 14 53 22 21 103 105
+$EndElements
@@ -0,0 +1,395 @@
+$MeshFormat
+2.2 0 8
+$EndMeshFormat
+$PhysicalNames
+9
+1 1 "inlet"
+1 2 "outlet"
+1 3 "top wall"
+1 4 "bottom wall"
+1 5 "beam wall"
+1 6 "cylinder"
+1 7 "beam cylinder curve"
+2 1 "fluid"
+2 2 "beam"
+$EndPhysicalNames
+$Nodes
+268
+1 0 0 0
+2 2.2 0 0
+3 2.2 0.41 0
+4 0 0.41 0
+5 0.2353553390593274 0.2353553390593274 0
+6 0.2353553390593274 0.1646446609406726 0
+7 0.1646446609406726 0.2353553390593274 0
+8 0.1646446609406726 0.1646446609406726 0
+9 0.2 0.25 0
+10 0.15 0.2 0
+11 0.2 0.15 0
+12 0.2489897948556636 0.19 0
+13 0.2489897948556636 0.21 0
+14 0.3310660171779822 0.3060660171779822 0
+15 0.3310660171779822 0.09393398282201787 0
+16 0.06893398282201785 0.3060660171779822 0
+17 0.06893398282201785 0.09393398282201787 0
+18 0.2 0.325 0
+19 0.04749999999999999 0.2 0
+20 0.2 0.07499999999999998 0
+21 0.6 0.21 0
+22 0.6 0.19 0
+23 0.3525 0.21 0
+24 0.3525 0.19 0
+25 0.2 0.41 0
+26 0.2 0 0
+27 0.06893398282201785 0.41 0
+28 0.06893398282201785 0 0
+29 0.3310660171779822 0 0
+30 0.3310660171779822 0.41 0
+31 0.6 0.41 0
+32 0.6 0 0
+33 0.6 0.3075 0
+34 0.6 0.1025 0
+35 2.2 0.3075 0
+36 2.2 0.1025 0
+37 2.2 0.21 0
+38 2.2 0.19 0
+39 0 0.2 0
+40 0 0.3060660171779822 0
+41 0 0.09393398282201787 0
+42 0.1808658283817455 0.2461939766255643 0
+43 0.1538060233744357 0.2191341716182545 0
+44 0.1538060233744357 0.1808658283817456 0
+45 0.1808658283817455 0.1538060233744357 0
+46 0.2191341716182545 0.1538060233744357 0
+47 0.2440369016676673 0.1763198122576657 0
+48 0.2440369016676674 0.2236801877423343 0
+49 0.2191341716182545 0.2461939766255643 0
+50 0.2 0.2875 0
+51 0.2 0.3675 0
+52 0.1344669914110089 0.41 0
+53 0.06893398282201785 0.3580330085889911 0
+54 0.1167893218813452 0.2707106781186548 0
+55 0.1344669914110089 0.3155330085889911 0
+56 0.2655330085889911 0.3155330085889911 0
+57 0.3310660171779822 0.3580330085889911 0
+58 0.2655330085889911 0.41 0
+59 0.2832106781186548 0.2707106781186548 0
+60 0.5103553390593271 0.41 0
+61 0.4207106781186545 0.41 0
+62 0.5551776695296635 0.41 0
+63 0.4655330085889908 0.41 0
+64 0.3758883476483184 0.41 0
+65 0.03446699141100892 0.41 0
+66 0 0.3580330085889911 0
+67 0.03446699141100892 0.3060660171779822 0
+68 0.05821699141100892 0.2530330085889911 0
+69 0.09875 0.2 0
+70 0.02374999999999999 0.2 0
+71 0 0.2530330085889911 0
+72 0 0.1469669914110089 0
+73 0.05821699141100892 0.1469669914110089 0
+74 0.03446699141100892 0.09393398282201787 0
+75 0 0.04696699141100893 0
+76 0.06893398282201785 0.04696699141100893 0
+77 0.03446699141100892 0 0
+78 0.1344669914110089 0 0
+79 0.2 0.03749999999999999 0
+80 0.1344669914110089 0.08446699141100893 0
+81 0.1167893218813452 0.1292893218813453 0
+82 0.2 0.1125 0
+83 0.2832106781186548 0.1292893218813453 0
+84 0.2655330085889911 0.08446699141100893 0
+85 0.2655330085889911 0 0
+86 0.3310660171779822 0.04696699141100893 0
+87 0.3417830085889911 0.1419669914110089 0
+88 0.3007448974278318 0.19 0
+89 0.3007448974278318 0.21 0
+90 0.3417830085889911 0.2580330085889911 0
+91 0.4207106781186546 0.3065440114519881 0
+92 0.5103553390593272 0.307022005725994 0
+93 0.3758883476483184 0.3063050143149851 0
+94 0.4655330085889909 0.306783008588991 0
+95 0.5551776695296636 0.307261002862997 0
+96 0.6 0.35875 0
+97 0.6 0.25875 0
+98 0.5174999999999998 0.21 0
+99 0.4349999999999998 0.21 0
+100 0.5587499999999999 0.21 0
+101 0.4762499999999998 0.21 0
+102 0.3937499999999999 0.21 0
+103 0.4350000000000002 0.19 0
+104 0.5175000000000003 0.19 0
+105 0.3937500000000001 0.19 0
+106 0.4762500000000002 0.19 0
+107 0.5587500000000001 0.19 0
+108 0.6 0.2 0
+109 0.6 0.14625 0
+110 0.5103553390593273 0.09964466094067262 0
+111 0.4207106781186548 0.09678932188134524 0
+112 0.5551776695296636 0.1010723304703363 0
+113 0.465533008588991 0.09821699141100892 0
+114 0.3758883476483185 0.09536165235168155 0
+115 0.420710678118655 0 0
+116 0.5103553390593276 0 0
+117 0.3758883476483186 0 0
+118 0.4655330085889914 0 0
+119 0.5551776695296639 0 0
+120 0.6 0.05125 0
+121 0.7769304775334696 0.41 0
+122 1.00694009683966 0.41 0
+123 1.305952603331686 0.41 0
+124 1.69466886165541 0.41 0
+125 0.6884652387667348 0.41 0
+126 0.8919352871865647 0.41 0
+127 1.156446350085673 0.41 0
+128 1.500310732493548 0.41 0
+129 1.947334430827705 0.41 0
+130 2.2 0.35875 0
+131 2.2 0.05125 0
+132 0.7769304775334696 0 0
+133 1.00694009683966 0 0
+134 1.305952603331686 0 0
+135 1.69466886165541 0 0
+136 0.6884652387667348 0 0
+137 0.8919352871865647 0 0
+138 1.156446350085673 0 0
+139 1.500310732493548 0 0
+140 1.947334430827705 0 0
+141 0.7769304775334696 0.1025 0
+142 1.00694009683966 0.1025 0
+143 1.305952603331686 0.1025 0
+144 1.69466886165541 0.1025 0
+145 0.6884652387667348 0.1025 0
+146 0.8919352871865647 0.1025 0
+147 1.156446350085673 0.1025 0
+148 1.500310732493548 0.1025 0
+149 1.947334430827705 0.1025 0
+150 0.7769304775334696 0.3075 0
+151 1.00694009683966 0.3075 0
+152 1.305952603331686 0.3075 0
+153 1.69466886165541 0.3075 0
+154 0.6884652387667348 0.3075 0
+155 0.8919352871865647 0.3075 0
+156 1.156446350085673 0.3075 0
+157 1.500310732493548 0.3075 0
+158 1.947334430827705 0.3075 0
+159 0.7769304775334696 0.21 0
+160 1.00694009683966 0.21 0
+161 1.305952603331686 0.21 0
+162 1.69466886165541 0.21 0
+163 0.6884652387667348 0.21 0
+164 0.8919352871865647 0.21 0
+165 1.156446350085673 0.21 0
+166 1.500310732493548 0.21 0
+167 1.947334430827705 0.21 0
+168 0.7769304775334696 0.19 0
+169 1.00694009683966 0.19 0
+170 1.305952603331686 0.19 0
+171 1.69466886165541 0.19 0
+172 0.6884652387667348 0.19 0
+173 0.8919352871865647 0.19 0
+174 1.156446350085673 0.19 0
+175 1.500310732493548 0.19 0
+176 1.947334430827705 0.19 0
+177 0.3525 0.2 0
+178 2.2 0.2 0
+179 2.2 0.25875 0
+180 2.2 0.14625 0
+181 0.25 0.2 0
+182 0.30125 0.2 0
+183 0.2929099551283292 0.2408565981656627 0
+184 0.4278553390593272 0.2582720057259941 0
+185 0.3848191738241591 0.2581525071574926 0
+186 0.5139276695296635 0.2585110028629971 0
+187 0.4708915042944953 0.2583915042944955 0
+188 0.5569638347648318 0.2586305014314985 0
+189 0.5175000000000001 0.2 0
+190 0.55875 0.2 0
+191 0.435 0.2 0
+192 0.47625 0.2 0
+193 0.39375 0.2 0
+194 0.2929099551283292 0.1591434018343373 0
+195 0.2423335901036228 0.2808634926072777 0
+196 0.2655330085889911 0.3627665042944955 0
+197 0.5103553390593272 0.358511002862997 0
+198 0.5551776695296635 0.3586305014314985 0
+199 0.4207106781186545 0.358272005725994 0
+200 0.4655330085889908 0.3583915042944955 0
+201 0.3758883476483184 0.3581525071574925 0
+202 0.1344669914110089 0.3627665042944955 0
+203 0.1576664098963772 0.2808634926072777 0
+204 0.03446699141100892 0.3580330085889911 0
+205 0.02910849570550446 0.2530330085889911 0
+206 0.1060115073927223 0.2360835901036228 0
+207 0.1060115073927223 0.1639164098963772 0
+208 0.03446699141100892 0.04696699141100893 0
+209 0.1576664098963772 0.1191365073927223 0
+210 0.1344669914110089 0.04223349570550446 0
+211 0.2423335901036228 0.1191365073927223 0
+212 0.2655330085889911 0.04223349570550446 0
+213 0.4278553390593275 0.1433946609406726 0
+214 0.3848191738241593 0.1426808261758408 0
+215 0.5139276695296637 0.1448223304703363 0
+216 0.4708915042944956 0.1441084957055045 0
+217 0.5569638347648319 0.1455361652351682 0
+218 0.5103553390593275 0.04982233047033631 0
+219 0.5551776695296637 0.05053616523516816 0
+220 0.4207106781186549 0.04839466094067262 0
+221 0.4655330085889912 0.04910849570550446 0
+222 0.3758883476483185 0.04768082617584078 0
+223 0.7769304775334696 0.2 0
+224 0.6884652387667349 0.2 0
+225 1.00694009683966 0.2 0
+226 0.8919352871865647 0.2 0
+227 1.305952603331686 0.2 0
+228 1.156446350085673 0.2 0
+229 1.69466886165541 0.2 0
+230 1.500310732493548 0.2 0
+231 1.947334430827705 0.2 0
+232 1.69466886165541 0.25875 0
+233 1.947334430827705 0.25875 0
+234 1.305952603331686 0.25875 0
+235 1.500310732493548 0.25875 0
+236 1.00694009683966 0.25875 0
+237 1.156446350085673 0.25875 0
+238 0.7769304775334696 0.25875 0
+239 0.8919352871865647 0.25875 0
+240 0.6884652387667347 0.25875 0
+241 0.7769304775334696 0.35875 0
+242 0.6884652387667349 0.35875 0
+243 1.00694009683966 0.35875 0
+244 0.8919352871865647 0.35875 0
+245 1.305952603331686 0.35875 0
+246 1.156446350085673 0.35875 0
+247 1.69466886165541 0.35875 0
+248 1.500310732493548 0.35875 0
+249 1.947334430827705 0.35875 0
+250 1.69466886165541 0.14625 0
+251 1.947334430827705 0.14625 0
+252 1.305952603331686 0.14625 0
+253 1.500310732493548 0.14625 0
+254 1.00694009683966 0.14625 0
+255 1.156446350085673 0.14625 0
+256 0.7769304775334696 0.14625 0
+257 0.8919352871865647 0.14625 0
+258 0.6884652387667349 0.14625 0
+259 0.7769304775334696 0.05125 0
+260 0.6884652387667349 0.05125 0
+261 1.00694009683966 0.05125 0
+262 0.8919352871865647 0.05125 0
+263 1.305952603331686 0.05125 0
+264 1.156446350085673 0.05125 0
+265 1.69466886165541 0.05125 0
+266 1.500310732493548 0.05125 0
+267 1.947334430827705 0.05125 0
+268 0.02910849570550446 0.1469669914110089 0
+$EndNodes
+$Elements
+106
+1 8 2 6 1 9 7 42
+2 8 2 6 2 7 10 43
+3 8 2 6 3 10 8 44
+4 8 2 6 4 8 11 45
+5 8 2 6 5 11 6 46
+6 8 2 6 6 6 12 47
+7 8 2 6 9 13 5 48
+8 8 2 6 10 5 9 49
+9 8 2 3 13 25 27 52
+10 8 2 3 19 30 25 58
+11 8 2 3 21 31 60 62
+12 8 2 3 21 60 61 63
+13 8 2 3 21 61 30 64
+14 8 2 3 22 27 4 65
+15 8 2 1 23 4 40 66
+16 8 2 1 28 39 40 71
+17 8 2 1 29 41 39 72
+18 8 2 1 32 41 1 75
+19 8 2 4 34 28 1 77
+20 8 2 4 35 28 26 78
+21 8 2 4 42 26 29 85
+22 8 2 5 45 24 12 88
+23 8 2 5 46 13 23 89
+24 8 2 5 51 21 98 100
+25 8 2 5 51 98 99 101
+26 8 2 5 51 99 23 102
+27 8 2 5 52 24 103 105
+28 8 2 5 52 103 104 106
+29 8 2 5 52 104 22 107
+30 8 2 5 53 22 21 108
+31 8 2 4 56 29 115 117
+32 8 2 4 56 115 116 118
+33 8 2 4 56 116 32 119
+34 8 2 3 58 31 121 125
+35 8 2 3 58 121 122 126
+36 8 2 3 58 122 123 127
+37 8 2 3 58 123 124 128
+38 8 2 3 58 124 3 129
+39 8 2 2 59 3 35 130
+40 8 2 2 61 36 2 131
+41 8 2 4 62 32 132 136
+42 8 2 4 62 132 133 137
+43 8 2 4 62 133 134 138
+44 8 2 4 62 134 135 139
+45 8 2 4 62 135 2 140
+46 8 2 2 68 37 38 178
+47 8 2 2 69 35 37 179
+48 8 2 2 70 38 36 180
+49 8 2 7 71 13 12 181
+50 10 2 2 1 13 23 24 12 89 177 88 181 182
+51 10 2 1 2 14 5 13 23 59 48 89 90 183
+52 10 2 1 3 23 14 91 99 90 93 184 102 185
+53 10 2 1 3 99 91 92 98 184 94 186 101 187
+54 10 2 1 3 98 92 33 21 186 95 97 100 188
+55 10 2 2 4 21 98 104 22 100 189 107 108 190
+56 10 2 2 4 98 99 103 104 101 191 106 189 192
+57 10 2 2 4 99 23 24 103 102 177 105 191 193
+58 10 2 1 5 15 24 12 6 87 88 47 83 194
+59 10 2 1 6 9 18 14 5 50 56 59 49 195
+60 10 2 1 7 18 25 30 14 51 58 57 56 196
+61 10 2 1 8 31 60 92 33 62 197 95 96 198
+62 10 2 1 8 60 61 91 92 63 199 94 197 200
+63 10 2 1 8 61 30 14 91 64 57 93 199 201
+64 10 2 1 9 18 25 27 16 51 52 53 55 202
+65 10 2 1 10 16 7 9 18 54 42 50 55 203
+66 10 2 1 11 27 16 40 4 53 67 66 65 204
+67 10 2 1 12 16 19 39 40 68 70 71 67 205
+68 10 2 1 13 19 10 7 16 69 43 54 68 206
+69 10 2 1 14 10 8 17 19 44 81 73 69 207
+70 10 2 1 16 41 1 28 17 75 77 76 74 208
+71 10 2 1 17 17 8 11 20 81 45 82 80 209
+72 10 2 1 18 20 17 28 26 80 76 78 79 210
+73 10 2 1 19 15 20 11 6 84 82 46 83 211
+74 10 2 1 20 15 20 26 29 84 79 85 86 212
+75 10 2 1 21 15 24 103 111 87 105 213 114 214
+76 10 2 1 21 111 103 104 110 213 106 215 113 216
+77 10 2 1 21 110 104 22 34 215 107 109 112 217
+78 10 2 1 22 34 32 116 110 120 119 218 112 219
+79 10 2 1 22 110 116 115 111 218 118 220 113 221
+80 10 2 1 22 111 115 29 15 220 117 86 114 222
+81 10 2 1 23 21 159 168 22 163 223 172 108 224
+82 10 2 1 23 159 160 169 168 164 225 173 223 226
+83 10 2 1 23 160 161 170 169 165 227 174 225 228
+84 10 2 1 23 161 162 171 170 166 229 175 227 230
+85 10 2 1 23 162 37 38 171 167 178 176 229 231
+86 10 2 1 24 35 37 162 153 179 167 232 158 233
+87 10 2 1 24 153 162 161 152 232 166 234 157 235
+88 10 2 1 24 152 161 160 151 234 165 236 156 237
+89 10 2 1 24 151 160 159 150 236 164 238 155 239
+90 10 2 1 24 150 159 21 33 238 163 97 154 240
+91 10 2 1 25 33 31 121 150 96 125 241 154 242
+92 10 2 1 25 150 121 122 151 241 126 243 155 244
+93 10 2 1 25 151 122 123 152 243 127 245 156 246
+94 10 2 1 25 152 123 124 153 245 128 247 157 248
+95 10 2 1 25 153 124 3 35 247 129 130 158 249
+96 10 2 1 26 38 36 144 171 180 149 250 176 251
+97 10 2 1 26 171 144 143 170 250 148 252 175 253
+98 10 2 1 26 170 143 142 169 252 147 254 174 255
+99 10 2 1 26 169 142 141 168 254 146 256 173 257
+100 10 2 1 26 168 141 34 22 256 145 109 172 258
+101 10 2 1 27 32 132 141 34 136 259 145 120 260
+102 10 2 1 27 132 133 142 141 137 261 146 259 262
+103 10 2 1 27 133 134 143 142 138 263 147 261 264
+104 10 2 1 27 134 135 144 143 139 265 148 263 266
+105 10 2 1 27 135 2 36 144 140 131 149 265 267
+106 10 2 1 28 17 19 39 41 73 70 72 74 268
+$EndElements
@@ -0,0 +1,960 @@
+/**
+ * Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+ * at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+ * LICENSE and NOTICE for details. LLNL-CODE-806117.
+ * 
+ * This file is part of the MFEM library. For more information and source code
+ * availability visit https://mfem.org.
+ * 
+ * MFEM is free software; you can redistribute it and/or modify it under the
+ * terms of the BSD-3 license. We welcome feedback and contributions, see file
+ * CONTRIBUTING.md for details.
+ * 
+ *            --------------------------------------------------
+ *                      Conjugate heat transfer miniapp
+ *            --------------------------------------------------
+ * 
+ * This is a miniapp demonstrates conjugate heat transfer by coupling different
+ * physics in different domains:
+ *    1) Incompressible Navier-Stokes equations in a fluid domain
+ *    2) Heat equation in fluid 
+ *    3) Heat equation solid domains
+ * 
+ * The test case is a benchmark Backward Facing Step (BFS) with a heated base.
+ * The following boundary conditions are applied:
+ *   1) Fluid inlet (attribute 1): parabolic velocity profile, T = 0.0
+ *   2) Fluid outlet (attribute 2): zero-pressure, -kappa grad(T)•n = 0.0
+ *   3) Fluid walls (attribute 3 & 4): no-slip, -kappa grad(T)•n = 0.0
+ *   4) Solid base (attribute 6): T = 1.0
+ *   5) Solid walls (attribute 5): -kappa grad(T)•n = 0.0
+ * 
+ *                                  4
+ *     -----------------------------------------------------------------
+ *   1 |                          fluid                                | 2
+ *   4 |                                                               | 2
+ *     ---------------------------- 3 ----------------------------------
+ *   5 |                                                               | 5
+ *   5 |                          solid                                | 5
+ *     -----------------------------------------------------------------
+ *                                  6
+ * 
+ * This example demonstrates nested coupling with the fluid flow and heat
+ * transfer solvers coupling where the heat transfer solver is itself a 
+ * coupled solver with the fluid and solid heat transfer solvers.
+ * 
+ * The fluid flow and heat transfer solvers are flow maps (i.e ODE Solvers) and 
+ * coupled with one-way, serial or parallel, coupling. The fluid and solid heat 
+ * transfer solvers can be partitioned-coupled or monolithically-coupled. 
+ * 
+ * For partitioned-coupling, the following boundary conditions are applied 
+ * at the fluid-solid interface:
+ *  1) Dirichlet: T = T_f  on attr 3 in solid domain
+ *  2) Neumann:   Q = -kappa grad(T_s)•n on attr 3 in fluid domain
+ * 
+ * For monolithic-coupling, the following conditions are imposed at the
+ * fluid-solid interface:
+ * 1) Continuity of temperature: T_f = T_s
+ * 2) Continuity of heat flux: -k_f grad(T_f)•n = k_s grad(T_s)•n
+ * 
+ * Sample run:
+ *    mpirun -np 6 cht-BFS -vs 500 -dt 1e-3 -tf 1000 -rs 2 -o 3 -ode 21 -scheme -1 -cht
+ */
+
+#include "mfem.hpp"
+#include "multiapp.hpp"
+#include "../navier/navier_solver.hpp"
+
+
+using namespace mfem;
+using namespace navier;
+using namespace std;
+
+struct BFSContext
+{
+   int ser_ref = 1;         // Serial mesh refinement
+   int order = 3;           // Finite element order
+   int ode_solver = 21;     // ODE solver
+   real_t dt = 1e-2;        // Time step size
+   real_t t_final = 3.0;    // Final time
+   int vis_steps = 100;     // Visualization steps
+   int couple_scheme = -1;  // Coupling scheme
+                            // -1: Monolithic, 0: Alternating Schwarz,
+                            // >0: Additive Schwarz with number of iterations
+   bool ht_only = false;     // Conjugate heat transfer on/off
+   bool visualization = true;// Visualization on/off
+
+   real_t Re = 800.0;         // Reynolds number
+   real_t Pr = 0.71;          // Prandtl number
+   real_t density = 1.0;      // Density
+   real_t kappa_ratio = 1e2;  // Conductivity ratio solid/fluid
+   bool checkres = false;     // Check results
+
+#if defined(MFEM_USE_DOUBLE)
+   real_t tol_T = 1e-4;
+   real_t tol_Q = 1e-4;
+#elif defined(MFEM_USE_SINGLE)
+   real_t tol_T = 1e-3;
+   real_t tol_Q = 1e-3;
+#else
+#error "Only single and double precision are supported!"
+   real_t tol_T = 0;
+   real_t tol_Q = 0;
+#endif   
+} ctx;
+
+void SetSolverParameters(IterativeSolver *solver, real_t rtol, real_t atol , int max_it, 
+                         int print_level, bool iterative_mode);
+
+/// Temperature profile 
+double temp_profile(const Vector& x)
+{
+    return x(1) == 0.0 ? 1.0 : 0.0;
+}
+
+/// Parabolic velocity profile for channel inlet
+void velocity_profile(const Vector &x, Vector &u)
+{
+   double xi = x(0), yi = x(1) - 2.5;
+   u = 0.0;
+   u(0) = 24.0*yi*(0.5-yi);
+}
+
+/**
+ * @brief Coefficient to compute normal heat flux Q = -kappa grad(T) • n
+ */
+class HeatFluxCoefficient : public Coefficient
+{
+protected:
+   ParGridFunction *T_gf = nullptr;
+   Coefficient *conductivity = nullptr;
+   Vector grad_T, normal;
+
+public:
+   HeatFluxCoefficient(ParGridFunction *T_gf_,
+                       Coefficient *conductivity_) : 
+                       Coefficient(), T_gf(T_gf_), 
+                       conductivity(conductivity_)
+                       {
+                        int dim = T_gf->FESpace()->GetMesh()->Dimension();
+                        grad_T.SetSize(dim);
+                        normal.SetSize(dim);
+                       }
+
+    real_t Eval(ElementTransformation &T,
+              const IntegrationPoint &ip) override
+    {
+        const DenseMatrix &jacobian = T.Jacobian();
+        
+        MFEM_ASSERT( (jacobian.Height() - 1 == jacobian.Width()), 
+                     "Incorrect Jacobian dimension. Coefficient only "
+                     "supported for boundary elements.");
+
+        CalcOrtho(jacobian, normal);
+        const double scale = normal.Norml2();
+        normal /= scale;
+
+        real_t kappa = conductivity ? conductivity->Eval(T, ip) : 1.0;
+        T_gf->GetGradient(T,grad_T);
+        real_t flux = -kappa*(grad_T * normal);
+        
+        return flux;
+    }                       
+};
+
+/**
+ * @brief Convection-diffusion time dependent operator
+ *
+ *              dT/dt = κΔT - α∇T•u
+ *
+ * Can also be used to create a diffusion or convection 
+ * only operator by setting α or κ to zero.
+ */
+class ConvectionDiffusion : public Application
+{
+public:
+
+   // Mesh and finite element space
+   ParMesh &mesh;
+   ParFiniteElementSpace &fes;
+
+   /// Essential and natural dof array.
+   Array<int> ess_attr, nat_attr;
+   Array<int> ess_tdofs, nat_tdofs;
+
+   /// Material properties
+   ConstantCoefficient diffusivity, kappa, alpha;
+
+   /// Grid functions for the temperature and heat flux
+   mutable ParGridFunction T_gf, Q_gf;
+
+   /// Used to store boundary condition data
+   mutable ParGridFunction T_gf_bc, Q_gf_bc; 
+   mutable GridFunctionCoefficient Q_fgc;
+   mutable HeatFluxCoefficient Q_coeff;
+
+   /// Mass form and Stiffness form. Might include 
+   /// diffusion, convection or both.
+   mutable ParBilinearForm Mform, Kform, Mform_e, Kform_e;
+
+   /// RHS form
+   mutable ParLinearForm bform;
+   mutable Vector b;
+
+   /// Mass and Stiffness operators
+   mutable HypreParMatrix Mmat, Kmat, Kmat_e, Mmat_e;
+
+   /// Mass matrix solver
+   CGSolver M_solver;
+   GMRESSolver implicit_solver;
+   HypreParMatrix *T = nullptr; // T = M + dt K
+
+   /// Preconditioners
+   HypreSmoother M_prec, T_prec;
+
+   /// Velocity coefficient
+   VectorCoefficient *velocity_coeff = nullptr;
+
+   /// Auxiliary variables
+   real_t current_dt = -1.0;
+   bool updated = false;
+   mutable Vector z, q;
+
+public:
+
+   ConvectionDiffusion(ParFiniteElementSpace &fes_,
+                        Array<int> ess_attr_,
+                        Array<int> nat_attr_,
+                        real_t diffusivity_ = 1.0,
+                        real_t kappa_ = 1.0,
+                        VectorCoefficient *velocity_coeff_ = nullptr,
+                        real_t alpha_ = 1.0)
+                        : Application(fes_.GetTrueVSize()),
+                        mesh(*fes_.GetParMesh()),
+                        fes(fes_),
+                        ess_attr(ess_attr_),
+                        nat_attr(nat_attr_),
+                        diffusivity(diffusivity_), 
+                        kappa(kappa_), alpha(-alpha_),
+                        T_gf(&fes), Q_gf(&fes),
+                        T_gf_bc(&fes), Q_gf_bc(&fes),
+                        Q_fgc(&Q_gf_bc), Q_coeff(&T_gf,&kappa),
+                        Mform(&fes), Kform(&fes),
+                        Mform_e(&fes), Kform_e(&fes),
+                        bform(&fes),
+                        M_solver(mesh.GetComm()),
+                        implicit_solver(mesh.GetComm()),
+                        velocity_coeff(velocity_coeff_)
+   {
+
+      fes.GetEssentialTrueDofs(ess_attr, ess_tdofs);
+      fes.GetEssentialTrueDofs(nat_attr, nat_tdofs);
+      
+      T_gf = 0.0;
+      Q_gf = 0.0;
+      T_gf_bc = 0.0;
+      Q_gf_bc = 0.0;
+
+      field_collection.AddSourceField("Temperature",&T_gf);
+      field_collection.AddSourceField("Flux",&Q_gf);
+      field_collection.AddField("Temperature_BC",&T_gf_bc);
+      field_collection.AddField("Flux_BC",&Q_gf_bc);
+
+      Mform.AddDomainIntegrator(new MassIntegrator);
+      Mform_e.AddDomainIntegrator(new MassIntegrator);
+
+      Kform.AddDomainIntegrator(new DiffusionIntegrator(diffusivity));
+      Kform_e.AddDomainIntegrator(new DiffusionIntegrator(diffusivity));
+
+      if(velocity_coeff)
+      {
+         Kform.AddDomainIntegrator(new ConvectionIntegrator(*velocity_coeff, alpha.constant));
+         Kform_e.AddDomainIntegrator(new ConvectionIntegrator(*velocity_coeff, alpha.constant));
+      }
+
+      if(nat_attr.Max() > 0)
+      {
+         bform.AddBoundaryIntegrator(new BoundaryLFIntegrator(Q_fgc),nat_attr);
+      }
+      
+      z.SetSize(fes.GetTrueVSize());
+      q.SetSize(fes.GetTrueVSize());
+      Assemble();   
+      BuildSolvers();
+   }
+
+   /// Assemble linear and bilinear forms; called if mesh is updated
+   void Assemble()
+   {
+      AssembleLinearForms();
+      AssembleBilinearForms();
+   }
+
+   void AssembleBilinearForms()
+   {
+      Mform.Assemble();
+      Kform.Assemble();
+      Mform_e.Assemble();
+      Kform_e.Assemble();
+
+      Array<int> empty;      
+      Mform.FormSystemMatrix(ess_tdofs, Mmat);
+      Kform.FormSystemMatrix(ess_tdofs, Kmat);
+      Mform_e.FormSystemMatrix(empty, Mmat_e);
+      Kform_e.FormSystemMatrix(empty, Kmat_e);
+   }
+
+   void AssembleLinearForms()
+   {
+      b.SetSize(fes.GetTrueVSize());
+      b = 0.0;
+      bform.Assemble();
+      bform.ParallelAssemble(b);
+   }
+
+   /// Update finite element space and re-assemble forms
+   /// if the mesh has changed
+   void Update() override 
+   {
+        fes.Update();
+
+        T_gf.Update();
+        T_gf_bc.Update();
+        Q_gf.Update();
+        Q_gf_bc.Update();
+    
+        Mform.Update();
+        Kform.Update();
+        Mform_e.Update();
+        Kform_e.Update();
+        bform.Update();
+
+        Assemble();
+        updated = true;
+   }
+
+   void BuildSolvers()
+   {
+      M_solver.iterative_mode = false;
+      M_solver.SetRelTol(1e-8);
+      M_solver.SetAbsTol(0.0);
+      M_solver.SetMaxIter(1000);
+      M_solver.SetPrintLevel(0);
+      M_prec.SetType(HypreSmoother::Jacobi);
+      M_solver.SetPreconditioner(M_prec);
+      M_solver.SetOperator(Mmat);
+
+      implicit_solver.iterative_mode = false;
+      implicit_solver.SetRelTol(1e-8);
+      implicit_solver.SetAbsTol(0.0);
+      implicit_solver.SetMaxIter(500);
+      implicit_solver.SetPrintLevel(0);
+      T_prec.SetType(HypreSmoother::Jacobi);
+      implicit_solver.SetPreconditioner(T_prec);
+   }
+
+   /// For explict time integration: k = M^{-1} ( -K u + b )
+   /// Used for partitioned coupling of explicit methods
+   void Mult(const Vector &u, Vector &k) const override
+   {      
+      Kmat.Mult(u, z);
+      z.Neg();
+      z.Add(1.0, b);
+      M_solver.Mult(z, k);
+      k.SetSubVector(ess_tdofs, 0.0);
+   }
+
+   /// For implicit time integration: k solves (M + dt K) k = -K u + b
+   /// Used for partitioned coupling of implicit methods
+   void ImplicitSolve(const real_t dt, const Vector &u, Vector &k)
+   {   
+      AssembleLinearForms();
+      if((current_dt != dt) || updated)
+      {
+         if (T) delete T;
+         T = Add(1.0, Mmat, dt, Kmat);
+         implicit_solver.SetOperator(*T);
+         updated = false;
+         current_dt = dt;
+      }
+
+      Kmat_e.Mult(u, z);
+      z.Neg();
+      z.Add(1.0, b);
+      implicit_solver.Mult(z, k);
+
+      if(IsCoupled() && nat_attr.Max() == 0)
+      {  // Apply interface conditions on temperature 
+         // Condition imposed by prescribing k = dT/dt
+         // T_gf_bc contains dT/dt from other domain
+         T_gf_bc.GetTrueDofs(z);
+         for (int i = 0; i < ess_tdofs.Size(); i++)
+         {
+            int idx = ess_tdofs[i];
+            k(idx) = z(idx);
+         }
+      }
+      else
+      {  // if uncoupled, apply standard essential BCs
+         k.SetSubVector(ess_tdofs, 0.0);
+      }
+   }
+
+   /// Computes the residual v = M k + K u - b + (Tf - Ts)_int + (Qf + Qs)_int
+   /// where ()_int represents the temperature and flux interface conditions
+   /// Used for monolithic coupling
+   void ImplicitMult(const Vector &u, const Vector &k, Vector &v ) const override
+   {
+      // v = M*k + K*u
+      Mmat_e.Mult(k, v);
+      Kmat_e.AddMult(u, v);
+
+      bform.Assemble();
+      bform.ParallelAssemble(b);
+      v.Add(-1.0,b); // v -= b
+
+      v.SetSubVector(ess_tdofs, 0.0); // Residual at uncoupled dofs is zero
+      
+      if(IsCoupled()) // Residual from interface conditions
+      {  // nat_dofs represent the coupled dofs
+         T_gf_bc.GetTrueDofs(z); // T from other domain
+         for (int i = 0; i < nat_tdofs.Size(); i++)
+         {
+               int idx = nat_tdofs[i];
+               v(idx) += (u(idx) - z(idx)); // T_f = T_s
+         }
+
+         Q_gf.GetTrueDofs(q); // Flux from T in this domain
+         Q_gf_bc.GetTrueDofs(z); // Flux from the other domain
+         for (int i = 0; i < nat_tdofs.Size(); i++)
+         {
+               int idx = nat_tdofs[i];
+               v(idx) += (q(idx) + z(idx)); // Q_f = -Q_s
+         }
+      }
+   }
+
+   /// Transfer temperature and flux boundary data
+   /// Called if application is monolithically coupled
+   void Transfer(const Vector &x) override
+   {
+      T_gf.SetFromTrueDofs(x);
+      Q_gf.ProjectBdrCoefficient(Q_coeff,ess_attr);
+      Q_gf.ProjectBdrCoefficient(Q_coeff,nat_attr);
+
+      Application::Transfer(); // Transfer all source fields
+   }   
+
+   /// Transfer temperature and flux boundary data
+   /// Called if application is partitioned coupled
+   void Transfer(const Vector &x, const Vector &k, real_t dt = 0.0) override
+   {
+      if(nat_attr.Max() == 0)
+      {  // This domain sends flux
+         add(1.0,x,dt,k,z); 
+         T_gf.SetFromTrueDofs(z);
+         Q_gf.ProjectBdrCoefficient(Q_coeff,ess_attr);
+         field_collection.Transfer("Flux"); // Only transfer flux
+      }
+      else
+      { // This domain is sending temperature
+         T_gf.SetFromTrueDofs(k);
+         field_collection.Transfer("Temperature"); // Only transfer temperature
+      }      
+   }   
+
+   ~ConvectionDiffusion() override
+   {
+      if(T) delete T;
+   }
+};
+
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init();
+   Hypre::Init();
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&ctx.order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&ctx.t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&ctx.dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&ctx.visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&ctx.vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.AddOption(&ctx.ode_solver, "-ode", "--ode-solver-type",
+                  "ODESolver id.");
+   args.AddOption(&ctx.ser_ref, "-rs", "--serial-refine",
+                  "Number of times to refine the mesh in serial.");                
+   args.AddOption(&ctx.ht_only, "-ht", "--heat-transfer-only",
+                  "-cht", "--conjugate-heat-transfer",
+                  "Conjugate heat transfer or heat transfer only.");
+   args.AddOption(&ctx.couple_scheme, "-scheme", "--coupling-scheme", 
+                 "Coupling scheme: -1 = Monolithic; 0 = Alt. Schw.; >0 = Add. Schw.");
+   args.AddOption(&ctx.checkres, "-cr", "--checkresult", "-no-cr", "--no-checkresult",
+                 "Enable or disable checking of the result. Returns -1 on failure.");
+   args.ParseCheck();
+
+
+   int order = ctx.order;
+   int ode_solver = ctx.ode_solver;
+
+   Mesh *serial_mesh = new Mesh("backward-facing-step.msh");
+   int dim = serial_mesh->Dimension();
+
+   for (int i = 0; i < ctx.ser_ref; ++i) { serial_mesh->UniformRefinement(); }
+   serial_mesh->SetCurvature(order, false, dim, Ordering::byNODES);
+
+   ParMesh parent_mesh = ParMesh(MPI_COMM_WORLD, *serial_mesh);
+   delete serial_mesh;
+   parent_mesh.UniformRefinement();
+   
+   
+   // Create the sub-domains and accompanying Finite Element spaces 
+   Array<int> domain_attributes(1);
+   domain_attributes[0] = 1;
+   auto fluid_mesh = ParSubMesh::CreateFromDomain(parent_mesh, domain_attributes);
+   fluid_mesh.SetAttributes();
+   fluid_mesh.EnsureNodes();
+   fluid_mesh.Finalize();
+
+   domain_attributes[0] = 2;
+   auto solid_mesh = ParSubMesh::CreateFromDomain(parent_mesh,domain_attributes);
+   solid_mesh.SetAttributes();
+   solid_mesh.EnsureNodes();
+   solid_mesh.Finalize();
+
+   // Set essential and natural boundary conditions attributes
+   Array<int> u_ess_attr, p_ess_attr, noslip_ess_attr;
+   Array<int> Ts_ess_attr, Ts_nat_attr;
+   Array<int> Tf_ess_attr, Tf_nat_attr;
+
+   if (solid_mesh.bdr_attributes.Size() > 0)
+   {
+      Ts_ess_attr.SetSize(solid_mesh.bdr_attributes.Max());
+      Ts_ess_attr = 0; 
+      Ts_ess_attr[2] = 1; // fluid-solid interface (T_fluid -> T_solid)
+      Ts_ess_attr[5] = 1; // bottom wall
+
+      Ts_nat_attr.SetSize(solid_mesh.bdr_attributes.Max());      
+      Ts_nat_attr = 0;
+      if(ctx.couple_scheme < 0)
+      { // If fully coupled, set nat. bc on interface for equality condition
+         Ts_nat_attr[2] = 1;
+      }
+   }
+   
+   if (fluid_mesh.bdr_attributes.Size() > 0)
+   {
+      u_ess_attr.SetSize(fluid_mesh.bdr_attributes.Max());
+      u_ess_attr = 0;
+      u_ess_attr[0] = 1; // inlet
+
+      p_ess_attr.SetSize(fluid_mesh.bdr_attributes.Max());
+      p_ess_attr = 0;
+      p_ess_attr[1] = 1; // outlet
+
+      noslip_ess_attr.SetSize(fluid_mesh.bdr_attributes.Max());
+      noslip_ess_attr = 1;
+      noslip_ess_attr[0] = 0; // inlet
+      noslip_ess_attr[1] = 0; // outlet
+
+      Tf_ess_attr.SetSize(fluid_mesh.bdr_attributes.Max());
+      Tf_ess_attr = 0;
+      Tf_ess_attr[0] = 1; // inlet
+
+      Tf_nat_attr.SetSize(fluid_mesh.bdr_attributes.Max());
+      Tf_nat_attr = 0;
+      Tf_nat_attr[2] = 1; // fluid-solid interface (Qs -> Qf)
+   }
+
+
+   // Finite element spaces for solid and fluid domains   
+   H1_FECollection ufec(order, dim);   // Velocity field (fluid domain)
+   H1_FECollection pfec(order-1, dim); // Pressure field (fluid domain)
+   H1_FECollection Tfec(order, dim);
+
+   ParFiniteElementSpace u_fes(&fluid_mesh, &ufec, dim, Ordering::byNODES);
+   ParFiniteElementSpace p_fes(&fluid_mesh, &pfec);   
+   ParFiniteElementSpace Tf_fes(&fluid_mesh, &Tfec);
+   ParFiniteElementSpace Tuf_fes(&fluid_mesh, &Tfec, dim, Ordering::byNODES);
+   ParFiniteElementSpace Ts_fes(&solid_mesh, &Tfec);   
+
+   // Set material properties
+   // From Backward Facing Step (BFS) Benchmark
+   real_t Re = ctx.Re;
+   real_t fluid_density = ctx.density;
+   real_t viscosity = fluid_density/Re;
+   real_t Pr = ctx.Pr ;
+   real_t fluid_alpha = 1.0;
+   real_t fluid_diffusivity = 1/(Re*Pr);
+   real_t fluid_kappa = fluid_diffusivity;
+   real_t solid_alpha = 0.0e0;
+   real_t kappa_ratio = ctx.kappa_ratio;
+   real_t solid_diffusivity = kappa_ratio*fluid_diffusivity;
+   real_t solid_kappa = fluid_kappa*kappa_ratio;
+
+   if(Mpi::Root())
+   {
+      std::cout << "Fluid Density: " << fluid_density << std::endl;
+      std::cout << "Viscosity: " << viscosity << std::endl;
+      std::cout << "Prandtl Number: " << Pr << std::endl;
+      std::cout << "Reynolds Number: " << Re << std::endl;
+      std::cout << "Solid Conductivity: " << solid_kappa << std::endl;
+      std::cout << "Fluid Diffusivity: " << fluid_diffusivity << std::endl;
+   }   
+
+   
+   // Navier miniapp
+   NavierSolver nse_miniapp(&fluid_mesh, order, viscosity);  
+
+   int max_bdf_order = 3;
+   nse_miniapp.EnablePA(true); 
+   nse_miniapp.SetMaxBDFOrder(max_bdf_order);
+
+   ParGridFunction &uf_gf = *nse_miniapp.GetCurrentVelocity();
+   ParGridFunction &p_gf = *nse_miniapp.GetCurrentPressure(); 
+   
+   Vector vzero(dim); vzero = 0.0;
+   ConstantCoefficient one_coeff(1.0);
+   Coefficient *zero_coeff = new ConstantCoefficient(0.0);   
+   VectorCoefficient *zerovec = new VectorConstantCoefficient(vzero);
+   VectorCoefficient *u_coeff = new VectorFunctionCoefficient(dim, velocity_profile);
+
+   // Set initial conditions in fluid    
+   p_gf.ProjectCoefficient(*zero_coeff);
+   uf_gf.ProjectCoefficient(*zerovec);
+   uf_gf.ProjectBdrCoefficient(*u_coeff,u_ess_attr);
+
+   nse_miniapp.AddVelDirichletBC(u_coeff, u_ess_attr);
+   nse_miniapp.AddVelDirichletBC(zerovec, noslip_ess_attr);
+   nse_miniapp.AddPresDirichletBC(zero_coeff, p_ess_attr);
+   nse_miniapp.Setup(ctx.dt);
+
+   /// Create navier block vector
+   Array<int> nse_offsets({0,uf_gf.ParFESpace()->GetTrueVSize(), p_gf.ParFESpace()->GetTrueVSize()});
+   nse_offsets.PartialSum();   
+   BlockVector up(nse_offsets);
+
+   uf_gf.GetTrueDofs(up.GetBlock(0));
+   p_gf.GetTrueDofs(up.GetBlock(1));
+
+
+   // Fluid Heat Transfer
+   ParGridFunction Tuf_gf(&Tuf_fes);  Tuf_gf = 0.0;
+   VectorGridFunctionCoefficient fluid_velocity(&Tuf_gf);
+   ConvectionDiffusion fluid_ht(Tf_fes, Tf_ess_attr, Tf_nat_attr, fluid_diffusivity, 
+                                fluid_kappa, &fluid_velocity, fluid_alpha);
+   
+   std::unique_ptr<ODESolver> Tf_odesolver = ODESolver::Select(ode_solver);
+   Tf_odesolver->Init(fluid_ht);
+
+   ParGridFunction &Tf_gf = *fluid_ht.Fields().GetField("Temperature");
+   ParGridFunction &Qf_gf = *fluid_ht.Fields().GetField("Flux");
+   ParGridFunction &Tf_gf_bc = *fluid_ht.Fields().GetField("Temperature_BC");
+   ParGridFunction &Qf_gf_bc = *fluid_ht.Fields().GetField("Flux_BC");
+
+   // Set initial conditions in fluid
+   Tf_gf.ProjectCoefficient(*zero_coeff);
+   Tf_gf.ProjectBdrCoefficient(*zero_coeff,Tf_ess_attr);
+
+   // Solid Heat Transfer
+   ConvectionDiffusion solid_ht(Ts_fes, Ts_ess_attr, Ts_nat_attr, solid_diffusivity, 
+                                solid_kappa, nullptr, solid_alpha);
+
+   std::unique_ptr<ODESolver> Ts_odesolver = ODESolver::Select(ode_solver);
+   Ts_odesolver->Init(solid_ht);
+
+   ParGridFunction &Ts_gf = *solid_ht.Fields().GetField("Temperature");
+   ParGridFunction &Qs_gf = *solid_ht.Fields().GetField("Flux");
+   ParGridFunction &Ts_gf_bc = *solid_ht.Fields().GetField("Temperature_BC");
+   ParGridFunction &Qs_gf_bc = *solid_ht.Fields().GetField("Flux_BC");
+
+   // Set initial conditions in solid
+   FunctionCoefficient temp_coeff(temp_profile);
+   Ts_gf.ProjectCoefficient(temp_coeff);
+
+
+   // Set up the coupled heat transfer multiapp
+   CoupledOperator ht_operator(2); // two coupled applications: fluid and solid heat transfer
+
+   Application* fl_ht_app = ht_operator.AddOperator(&fluid_ht);
+   Application* sl_ht_app = ht_operator.AddOperator(&solid_ht);
+
+   fl_ht_app->SetCoupled(true);
+   sl_ht_app->SetCoupled(true);
+
+   // ODE solver for the coupled operator   
+   std::unique_ptr<ODESolver> ht_odesolver = ODESolver::Select(ode_solver);
+   ht_odesolver->Init(ht_operator);
+
+   Array<int> ht_offsets({0,Tf_fes.GetTrueVSize(), Ts_fes.GetTrueVSize()});
+   ht_offsets.PartialSum();   
+   BlockVector Tfsv(ht_offsets);
+
+   Tf_gf.GetTrueDofs(Tfsv.GetBlock( fl_ht_app->GetOperatorIndex() ));
+   Ts_gf.GetTrueDofs(Tfsv.GetBlock( sl_ht_app->GetOperatorIndex() ));
+   
+   // Set up conjugate heat transfer app
+   CoupledOperator cht_app(2); // two coupled applications: navier, fluid-solid heat transfer
+
+   Application* nse_app = cht_app.AddOperator(&nse_miniapp,up.Size());  // (type-erased) Navier miniapp
+   Application* ht_app  = cht_app.AddOperator(ht_odesolver.get()); // Coupled fluid-solid heat transfer ODESolver;
+
+
+   // Set up field transfers
+   NativeTransfer Tf_Ts_map(Tf_fes, Ts_fes); // default map if none provided
+   // GSLibTransfer Tf_Ts_map(Tf_fes, Ts_fes);
+
+   /// Create link between fields in different apps
+   LinkedFields uf_to_Tuf_lf(&uf_gf, &Tuf_gf); // Navier velocity to fluid-heat convection velocity
+   LinkedFields Tf_to_Ts_lf(&Tf_gf, &Ts_gf_bc, &Tf_Ts_map);  // Fluid temperature to solid temperature
+   LinkedFields Qs_to_Qf_lf(&Qs_gf, &Qf_gf_bc);  // Solid heat flux to fluid heat flux
+
+   /// Different methods for adding the linked field to their source apps
+   nse_app->AddLinkedFields("Velocity",&uf_to_Tuf_lf);
+   fl_ht_app->AddLinkedFields("Temperature", &Tf_to_Ts_lf);
+   fl_ht_app->Fields().AddTargetField("Flux", &Qs_gf_bc, &Tf_Ts_map);
+   sl_ht_app->Fields().AddLinkedFields("Flux", &Qs_to_Qf_lf);
+   sl_ht_app->Fields().AddTargetField("Temperature", &Tf_gf_bc);
+
+   // Solvers
+   // Select solver for partitioned coupling
+   FPISolver fp_solver(MPI_COMM_WORLD); // For partitioned solves
+   AitkenRelaxation fp_relax;
+   SetSolverParameters(&fp_solver, 0.0, 5e-4, 500, 1, false);
+   fp_relax.SetBounds(0.0,1.0e-1);
+   fp_solver.SetRelaxation(5e-1, nullptr); // Use default relaxation method
+   // fp_solver.SetRelaxation(1e-1, &fp_relax);
+
+   // Select solver for monolithic/full coupling
+   GMRESSolver gmres_solver(MPI_COMM_WORLD); 
+   SetSolverParameters(&gmres_solver, 1e-7, 1e-7, 500, 0, false);
+   gmres_solver.SetKDim(300);
+   
+   NewtonSolver newton_solver(MPI_COMM_WORLD);  
+   SetSolverParameters(&newton_solver, 0.0, 1e-7, 100, 0, false);
+   newton_solver.SetSolver(gmres_solver);
+
+   /// Set coupling scheme and corresponding solvers   
+   /// The Navier miniapp and coupled heat transfer ODESolver (both flow maps) are coupled 
+   /// can be coupled in parallel (Additive Schwarz) or serial (Alternating Schwarz) but not
+   /// monolithically. The fluid and solid heat transfer Applications in the coupled heat 
+   /// transfer ODESolver can also be coupled monolithically.
+   if(ctx.couple_scheme == -1)
+   {
+      cht_app.SetCouplingScheme(CoupledOperator::Scheme::ALTERNATING_SCHWARZ);
+      ht_operator.SetCouplingScheme(CoupledOperator::Scheme::MONOLITHIC);
+      ht_operator.SetSolver(&newton_solver);
+   }
+   else if(ctx.couple_scheme == 0)
+   {
+      cht_app.SetCouplingScheme(CoupledOperator::Scheme::ALTERNATING_SCHWARZ);
+      ht_operator.SetCouplingScheme(CoupledOperator::Scheme::ALTERNATING_SCHWARZ);
+      ht_operator.SetSolver(&fp_solver);
+   }
+   else
+   {
+      cht_app.SetCouplingScheme(CoupledOperator::Scheme::ADDITIVE_SCHWARZ);
+      ht_operator.SetCouplingScheme(CoupledOperator::Scheme::ADDITIVE_SCHWARZ);
+      ht_operator.SetSolver(&fp_solver);
+   }
+
+   ht_operator.Assemble(false);
+   cht_app.Assemble(false);
+
+
+   auto nse_preprocess  = [&nse_offsets, &uf_gf, &p_gf](Vector &x) mutable { 
+         BlockVector up(x.GetData(), nse_offsets);
+         uf_gf.GetTrueDofs(up.GetBlock(0));
+         p_gf.GetTrueDofs(up.GetBlock(1));
+   };
+   auto nse_postprocess = [&nse_offsets, &uf_gf, &p_gf](Vector &x) mutable { 
+         BlockVector up(x.GetData(), nse_offsets);
+         uf_gf.SetFromTrueDofs(up.GetBlock(0));
+         p_gf.SetFromTrueDofs(up.GetBlock(1));
+   };
+
+   auto ht_preprocess = [&ht_offsets, &Tf_gf, &Ts_gf](Vector &x) mutable { 
+         BlockVector Tb(x.GetData(), ht_offsets);
+         Tf_gf.GetTrueDofs(Tb.GetBlock(0));
+         Ts_gf.GetTrueDofs(Tb.GetBlock(1));
+   };
+   auto ht_postprocess = [&ht_offsets, &Tf_gf, &Ts_gf](Vector &x) mutable { 
+         BlockVector Tb(x.GetData(), ht_offsets);
+         Tf_gf.SetFromTrueDofs(Tb.GetBlock(0));
+         Ts_gf.SetFromTrueDofs(Tb.GetBlock(1));
+   };
+
+   /// Set pre/post processing lambdas to corresponding apps
+   ht_app->SetPreProcessFunction(ht_preprocess);
+   ht_app->SetPostProcessFunction(ht_postprocess);
+
+   /// Not strictly necessary since Navier owns and updates GridFunctions 
+   /// internally but included here for completeness
+   nse_app->SetPreProcessFunction(nse_preprocess);
+   nse_app->SetPostProcessFunction(nse_postprocess);
+
+
+   // Set up the initial conditions in block vector for the
+   // coupled application in the correct order
+   int fl_id = nse_app->GetOperatorIndex();
+   int ht_id = ht_app->GetOperatorIndex();
+
+   BlockVector xb(cht_app.GetBlockOffsets());
+   xb.GetBlock(fl_id) = up;  // Fluid velocity and pressure    
+   xb.GetBlock(ht_id) = Tfsv;
+
+   
+   // Set up visualization 
+   ParaViewDataCollection *fluid_pv = nullptr;
+   ParaViewDataCollection *solid_pv = nullptr;
+
+   if(ctx.visualization)
+   {
+      fluid_pv = new ParaViewDataCollection("cht-BFS-fluid", &fluid_mesh);
+      solid_pv = new ParaViewDataCollection("cht-BFS-solid", &solid_mesh);
+
+      fluid_pv->SetLevelsOfDetail(order);
+      fluid_pv->SetDataFormat(VTKFormat::BINARY);
+      fluid_pv->SetHighOrderOutput(true);
+      fluid_pv->RegisterField("pressure",&p_gf);
+      fluid_pv->RegisterField("velocity",&uf_gf);
+      fluid_pv->RegisterField("convection",&Tuf_gf);
+      fluid_pv->RegisterField("Temperature",&Tf_gf);
+      fluid_pv->RegisterField("Flux",&Qf_gf);
+      
+      solid_pv->SetLevelsOfDetail(order);
+      solid_pv->SetDataFormat(VTKFormat::BINARY);
+      solid_pv->SetHighOrderOutput(true);
+      solid_pv->RegisterField("Temperature",&Ts_gf);
+      solid_pv->RegisterField("Flux",&Qs_gf);      
+   }
+      
+   auto save_callback = [&](int cycle, double t)
+   {
+      if(fluid_pv)
+      {
+         fluid_pv->SetCycle(cycle);
+         fluid_pv->SetTime(t);
+         fluid_pv->Save();
+      }
+
+      if(solid_pv)
+      {
+         solid_pv->SetCycle(cycle);
+         solid_pv->SetTime(t);
+         solid_pv->Save();
+      }
+   };
+   
+
+   if (Mpi::Root()) { 
+      out << "Starting time integration..." << std::endl; 
+   }
+
+   StopWatch timer;
+   timer.Start();
+
+   real_t t = 0.0;
+   bool last_step = false;
+   int tindex = 1;
+   save_callback(0, t);
+
+   last_step = false;
+   for (; !last_step; tindex++)
+   {
+      if (t + ctx.dt >= ctx.t_final - ctx.dt/2){ last_step = true; }
+
+      if(ctx.ht_only)
+      {
+         ht_odesolver->Step(Tfsv,t,ctx.dt);
+         Tf_gf.SetFromTrueDofs(Tfsv.GetBlock(0));
+         Ts_gf.SetFromTrueDofs(Tfsv.GetBlock(1));
+      }
+      else
+      {
+         cht_app.Step(xb, t, ctx.dt);
+      }
+
+      if (last_step || (tindex % ctx.vis_steps) == 0){
+         if (Mpi::Root()) { out << "step " << tindex << ", t = " << t << std::endl;}
+         save_callback(tindex, t);
+      }
+   }
+
+   timer.Stop();
+   if (Mpi::Root()){
+      out << "Total time: " << timer.RealTime() << " seconds." << std::endl;
+   }
+
+   /// Compute interface error
+   if(ctx.checkres)
+   {
+      Array<int> fl_int_attr(fluid_mesh.bdr_attributes.Max());
+      fl_int_attr[2] = 1; // fluid-solid interface
+      
+      /// Create submesh and FE space for the interface
+      ParSubMesh int_mesh = ParSubMesh::CreateFromBoundary(fluid_mesh, fl_int_attr);
+      ParFiniteElementSpace int_fes(&int_mesh, &Tfec);
+
+      /// Receiving grid functions on the interface
+      ParGridFunction fl_int_gf(&int_fes);
+      ParGridFunction sl_int_gf(&int_fes);
+
+      /// Maps and linked fields to transfer domain grid functions to the interface
+      NativeTransfer fl_to_int_submesh(Tf_fes, int_fes);
+      LinkedFields TQf_to_Tint_lf(&Tf_gf, &fl_int_gf, &fl_to_int_submesh);
+      LinkedFields TQs_to_Tint_lf(&Tf_gf_bc, &sl_int_gf, &fl_to_int_submesh);
+
+      cht_app.Transfer(xb); // Transfer all fields *_gf to their target fields *_gf_bc
+      
+      fl_int_gf = 0.0; sl_int_gf = 0.0;
+      TQf_to_Tint_lf.Transfer(); // Transfer Tf to interface
+      TQs_to_Tint_lf.Transfer(); // Transfer Ts (in Tf_gf_bc) to interface
+      real_t err_T = sqrt(DistanceSquared(int_mesh.GetComm(), fl_int_gf, sl_int_gf));
+
+      // Update sources in existing linked fields (can also create new ones)
+      TQf_to_Tint_lf.SetSource(&Qf_gf);
+      TQs_to_Tint_lf.SetSource(&Qf_gf_bc);
+
+      fl_int_gf = 0.0; sl_int_gf = 0.0;
+      TQf_to_Tint_lf.Transfer(); // Transfer Qf to interface
+      TQs_to_Tint_lf.Transfer(); // Transfer Qs (in Qf_gf_bc) to interface
+      real_t err_Q = sqrt(DistanceSquared(int_mesh.GetComm(), fl_int_gf, sl_int_gf));
+
+      if (sqrt(err_T) > ctx.tol_T || sqrt(err_Q) > ctx.tol_Q)
+      {
+         if (Mpi::Root())
+         {
+            mfem::out << "Result has a larger error than expected."
+                      << "T Error = " << sqrt(err_T)
+                      << ", Q Error = " << sqrt(err_Q)
+                      << std::endl;
+         }
+         return -1;
+      }
+   }
+
+   if(fluid_pv) delete fluid_pv;
+   if(solid_pv) delete solid_pv;
+   delete zero_coeff;
+   delete zerovec;
+   delete u_coeff;
+
+   return 0;
+}
+
+void SetSolverParameters(IterativeSolver *solver, real_t rtol, real_t atol , int max_it, 
+                         int print_level, bool iterative_mode)
+{
+    solver->SetRelTol(rtol);
+    solver->SetAbsTol(atol);
+    solver->SetMaxIter(max_it);
+    solver->SetPrintLevel(print_level);
+    solver->iterative_mode = iterative_mode;
+}
@@ -0,0 +1,616 @@
+/**
+ * Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+ * at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+ * LICENSE and NOTICE for details. LLNL-CODE-806117.
+ * 
+ * This file is part of the MFEM library. For more information and source code
+ * availability visit https://mfem.org.
+ * 
+ * MFEM is free software; you can redistribute it and/or modify it under the
+ * terms of the BSD-3 license. We welcome feedback and contributions, see file
+ * CONTRIBUTING.md for details.
+ * 
+ *            --------------------------------------------------
+ *                    Fluid-Structure Interaction miniapp
+ *            --------------------------------------------------
+ * 
+ * This miniapp simulates fluid-structure interaction (FSI) problems described in
+ * the paper:
+ *   
+ * with the incompressible Navier-Stokes equations, in Arbitrary Lagrangian-Eulerian (ALE)
+ * formulation, coupled with linear elasticity equations in the solid domain. 
+ * The coupling is done with a partitioned approach using the alternating or additive
+ * Schwarz method. The fluid mesh motion is handled with a mesh displacement diffusion
+ * approach. The geometry is a channel with a cylinder and a flexible beam attached to
+ * the cylinder downstream. 
+ * 
+ * The following boundary conditions are applied:
+ *  1) Channel inlet (attribute 1): parabolic velocity profile
+ *  2) Channel outlet (attribute 2): zero-pressure
+ *  3) Channel walls: no-slip
+ * 
+ * with the following interface conditions at the fluid-structure interface:
+ * 1) Continuity of velocity: u_f = u_s
+ * 2) Continuity of traction: -*pI + mu(grad(u_f)+grad(u_f)^T))•n = sigma_s•n
+ * 
+ * The velocity continuity is imposed using the stage-slope in the implicit multistage method
+ * ku_f = du_f/dt = du_s/dt = ku_s, where ku_f and ku_s are the fluid and solid stage-slopes
+ * 
+ * The mesh morphing is modeled as a displacement diffusion equation, 
+ * dx/dt = κΔx, with the ALE mesh velocity, w = dx/dt.
+ * 
+ * Sample run:
+ *    mpirun -np 6 ./fsi -vs 5 -dt 5e-3 -tf 10 -o 3 -rs 2 -ode 21 -U 1.0 -cs 1 -idir fsi-turek
+ */
+
+#include "mfem.hpp"
+#include "multiapp.hpp"
+#include "apps/navier_stokes.hpp"
+#include "apps/elasticity.hpp"
+#include "apps/mesh_morpher.hpp"
+
+#include <filesystem>
+
+using namespace mfem;
+using namespace std;
+
+// mpirun -np 6 ./fsi -vs 5 -dt 1e-2 -tf 5 -o 2 -rs 2 -ode 21
+// mpirun -np 6 ./fsi -vs 5 -dt 5e-3 -tf 10 -o 3 -rs 2 -ode 21 -U 1.0 -cs 0 -init
+
+void SetSolverParameters(IterativeSolver *solver, real_t rtol, real_t atol , int max_it, 
+                         int print_level, bool iterative_mode);
+
+bool ReadGridFunctionFromFile(const string &dirname, const string &gf_name,
+                              ParMesh &mesh, ParGridFunction &gf);
+
+void CollectArrays(std::vector<Array<int>*> &dof_arrays, Array<int> &tdof_array);
+
+int main(int argc, char *argv[])
+{
+    Mpi::Init();
+    Hypre::Init();
+    int num_procs = Mpi::WorldSize();
+    int myid = Mpi::WorldRank();
+
+    int order = 2;
+    int ser_ref = 0;    
+    int ode_solver_type = 21; // (21) BackwardEulerSolver
+                              // (22) SDIRK23Solver
+                              // (23) SDIRK33Solver
+                              // (34) SDIRK34Solver
+    real_t Uavg = 2.0;
+    real_t t_dev = 2.0;
+    real_t t_final = 1.0;
+    real_t dt = 1.0e-3;
+    real_t relax_factor = 1.0;
+    int couple_scheme = 1;
+    bool init = false;
+    bool lsave = false;
+    bool visualization = true;
+    int vis_steps = 10;
+    std::string init_dir= "";
+
+
+    OptionsParser args(argc, argv);
+    args.AddOption(&order, "-o", "--order",
+                    "Finite element order (polynomial degree).");
+    args.AddOption(&t_final, "-tf", "--t-final",
+                    "Final time; start time is 0.");
+    args.AddOption(&dt, "-dt", "--time-step",
+                    "Time step.");
+    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                    "--no-visualization",
+                    "Enable or disable GLVis visualization.");
+    args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                    "Visualize every n-th timestep.");
+    args.AddOption(&ode_solver_type, "-ode", "--ode-solver-type",
+                    "ODESolver id.");
+    args.AddOption(&ser_ref, "-rs", "--serial-refine",
+                    "Number of times to refine the mesh in serial.");
+    args.AddOption(&Uavg, "-U", "--velocity", "Mean velocity.");
+    args.AddOption(&t_dev, "-t_dev", "--transistion", "Developed flow transition time.");
+    args.AddOption(&couple_scheme, "-cs", "--coupling-scheme", 
+                    "Coupling scheme: -1 = Monolithic; 0 = Add. Schw.; >0 = Alt. Schw.");
+    args.AddOption(&relax_factor, "-rf", "--relaxation-factor",
+                    "Initial relaxation factor for the FPI solver.");    
+    args.AddOption(&init_dir, "-idir", "--init-directory", 
+                    "Directory containing intialization files. If doesn't exist, used to write init files.");
+    args.AddOption(&lsave, "-save", "--save-init", "-no-save",
+                    "--no-save", "Enable or disable saving initialization files.");
+    args.ParseCheck();
+
+
+    init = !init_dir.empty(); // if init_dir is provided, then init = true
+    
+    std::string mesh_file = "channel-cylinder.msh";
+    Mesh *serial_mesh = new Mesh(mesh_file);
+    int dim = serial_mesh->Dimension();
+
+    for (int i = 0; i < ser_ref; ++i) { serial_mesh->UniformRefinement(); }
+    serial_mesh->SetCurvature(order, false, dim, Ordering::byNODES);
+    serial_mesh->EnsureNCMesh();
+    ParMesh parent_mesh = ParMesh(MPI_COMM_WORLD, *serial_mesh);
+    delete serial_mesh;
+
+    // Create mesh (sub)domains
+    Array<int> domain_attributes(1);
+
+    // Create submesh for solid
+    domain_attributes[0] = 2;
+    auto solid_mesh = ParSubMesh::CreateFromDomain(parent_mesh, domain_attributes);
+    solid_mesh.SetAttributes();
+    solid_mesh.EnsureNodes();
+
+    // Create submesh for fluid
+    domain_attributes[0] = 1;
+    auto fluid_mesh = ParSubMesh::CreateFromDomain(parent_mesh, domain_attributes);
+    fluid_mesh.SetAttributes();
+    fluid_mesh.EnsureNodes();
+
+
+    // Set essential and natural boundary conditions attributes
+    Array<int> u_ess_attr, p_ess_attr, noslip_attr;
+    Array<int> xs_ess_attr, xs_nat_attr;
+    Array<int> xf_ess_attr;
+    Array<int> empty;
+
+    if (solid_mesh.bdr_attributes.Size() > 0)
+    {
+        xs_ess_attr.SetSize(solid_mesh.bdr_attributes.Max());
+        xs_nat_attr.SetSize(solid_mesh.bdr_attributes.Max());
+        xs_ess_attr = 0; xs_nat_attr = 0;
+        xs_nat_attr[4] = 1; // beam wall
+        xs_ess_attr[6] = 1; // beam cylinder curve
+    }
+
+    if (fluid_mesh.bdr_attributes.Size() > 0)
+    {
+        u_ess_attr.SetSize(fluid_mesh.bdr_attributes.Max());
+        p_ess_attr.SetSize(fluid_mesh.bdr_attributes.Max());
+        noslip_attr.SetSize(fluid_mesh.bdr_attributes.Max());
+        xf_ess_attr.SetSize(fluid_mesh.bdr_attributes.Max());
+
+        xf_ess_attr = 1;
+
+        u_ess_attr = 1;
+        u_ess_attr[1] = 0; // outlet
+
+        p_ess_attr = 0;
+        p_ess_attr[1] = 1; // outlet
+
+        noslip_attr = 1;
+        noslip_attr[0] = 0; // inlet
+        noslip_attr[1] = 0; // outlet
+    }
+
+
+    // Finite element spaces for solid and fluid domains
+    H1_FECollection xfec(order, dim); // Displacement field (solid and fluid domains)
+    H1_FECollection ufec(order, dim); // Velocity field (fluid domain)
+    H1_FECollection pfec(order-1, dim); // Pressure field (fluid domain)
+
+    ParFiniteElementSpace xs_fes(&solid_mesh, &xfec, dim, Ordering::byNODES);
+    ParFiniteElementSpace xf_fes(&fluid_mesh, &xfec, dim, Ordering::byNODES);
+    ParFiniteElementSpace u_fes(&fluid_mesh, &ufec, dim, Ordering::byNODES);
+    ParFiniteElementSpace p_fes(&fluid_mesh, &pfec);
+
+
+    Vector vzero(dim); vzero = 0.0;
+    VectorConstantCoefficient zerovec(vzero);
+
+    // Inlet velocity boundary condition
+    auto velocity_profile = [&Uavg, &t_dev](const Vector &x, double t, Vector &u) mutable
+    {
+        double xi = x(0), yi = x(1);
+        double U = Uavg;
+        double ramp_time = 2.0;
+        u = 0.0;
+        if (xi == 0.0)
+        {
+            u(0) = 1.5 * 4.0 * U * yi * (0.41 - yi) / (pow(0.41, 2.0));
+            u(1) = 0.0;
+        }
+        if(t < ramp_time) u(0) *= 0.5*(1.0 - cos(0.5*M_PI*t));
+    };
+
+
+    // Set material properties
+    // real_t E = 5.6e6; // Young's modulus
+    real_t E = 1.4e6; // Young's modulus
+    real_t nu = 4.0e-1; // Poisson's ratio
+    real_t solid_density = 1.0e4;
+    real_t viscosity = 1.0e-3;
+    real_t fluid_density = 1.0e3;
+    real_t mesh_diffusion = 1.0e0;
+    real_t lame_mu = E / (2.0 * (1.0 + nu));
+    real_t lame_lambda = 2.0*lame_mu*nu/(1.0-2.0*nu);
+    real_t compressibility = 0e-4; // Artificial Compressibility
+    bool scaled_pressure = true;
+
+    if(Mpi::Root())
+    {
+        std::cout << "Young's Modulus: " << E << std::endl;
+        std::cout << "Poisson's Ratio: " << nu << std::endl;
+        std::cout << "Lame's Mu: " << lame_mu << std::endl;
+        std::cout << "Lame's Lambda: " << lame_lambda << std::endl;
+        std::cout << "Solid Density: " << solid_density << std::endl;
+        std::cout << "Fluid Density: " << fluid_density << std::endl;
+        std::cout << "Viscosity: " << viscosity << std::endl;
+        std::cout << "Mesh Diffusion: " << mesh_diffusion << std::endl;
+        std::cout << "Artificial Compressibility: " << compressibility << std::endl;
+        std::cout << "Mean Velocity: " << Uavg << std::endl;
+        std::cout << "Reynolds Number: " << (Uavg*0.1)/viscosity << std::endl;
+        std::cout << "AE Number: " << (E)/(Uavg*Uavg*fluid_density) << std::endl;
+        std::cout << "Beta Number: " << solid_density/fluid_density << std::endl;
+    }
+
+    // Build individual applications (morpher, elasticity, navier-stokes)
+
+    // Mesh Morpher
+    MeshDiffusion morpher(xf_fes, xf_ess_attr, mesh_diffusion);
+
+    // Morphing grid functions
+    ParGridFunction &xf_gf  = *morpher.Fields().GetField("Displacement");
+    ParGridFunction &dxf_gf = *morpher.Fields().GetField("dxdt");
+    ParGridFunction &xf_gf_bc = *morpher.Fields().GetField("Displacement_BC");
+    ParGridFunction &dxf_gf_bc = *morpher.Fields().GetField("Velocity_BC");
+
+    /// Morphing solution
+    Vector xf(xf_fes.GetTrueVSize());
+    xf_gf.GetTrueDofs(xf);
+
+    ParGridFunction mesh_disp(&xf_fes), mesh_vel(&xf_fes);
+    mesh_disp = 0.0; mesh_vel = 0.0;
+
+
+    // Elasticity
+    Elasticity elasticity(xs_fes, xs_ess_attr, xs_nat_attr, solid_density, lame_mu, lame_lambda);
+
+    // Elasticity grid functions
+    ParGridFunction &xs_gf = *elasticity.Fields().GetField("Displacement");
+    ParGridFunction &us_gf = *elasticity.Fields().GetField("Velocity");
+    ParGridFunction &stress_gf = *elasticity.Fields().GetField("Traction");
+
+    /// Solid solution
+    Array<int> solid_offsets({0,xs_fes.GetTrueVSize(), xs_fes.GetTrueVSize()});
+    solid_offsets.PartialSum();
+    BlockVector xu(solid_offsets); // Solid displacement and velocity
+    us_gf.GetTrueDofs(xu.GetBlock(0));
+    xs_gf.GetTrueDofs(xu.GetBlock(1));
+
+
+    // Navier-Stokes
+    VectorGridFunctionCoefficient ale_uf(&dxf_gf);
+    ScalarVectorProductCoefficient neg_ale(-1.0, ale_uf);
+    VectorCoefficient *ale_velocity = &neg_ale;
+
+    NavierStokes nse(u_fes, p_fes, u_ess_attr, p_ess_attr,
+                     fluid_density, viscosity, compressibility,
+                     scaled_pressure, ale_velocity);
+
+    std::unique_ptr<ODESolver> nse_solver = ODESolver::Select(ode_solver_type);
+    nse_solver->Init(nse);
+
+    // Navier-Stokes grid functions
+    ParGridFunction &p_gf = *nse.Fields().GetField("Pressure");
+    ParGridFunction &uf_gf = *nse.Fields().GetField("Velocity");
+    ParGridFunction &tau_gf = *nse.Fields().GetField("Stress");
+    ParGridFunction &uf_gf_bc = *nse.Fields().GetField("Velocity_BC");
+
+    ConstantCoefficient p_coeff(0.0);
+    VectorFunctionCoefficient u_coeff(dim, velocity_profile);
+    
+    // Set initial conditions in fluid
+    p_gf.ProjectCoefficient(p_coeff);
+    uf_gf.ProjectCoefficient(u_coeff);
+    uf_gf.ProjectBdrCoefficient(zerovec,noslip_attr);
+
+    if(init) // Initialize from file
+    {
+        ReadGridFunctionFromFile(init_dir, "p-init.gf", fluid_mesh, p_gf);
+        ReadGridFunctionFromFile(init_dir, "u-init.gf", fluid_mesh, uf_gf);
+        u_coeff.SetTime(t_dev+2.0); // time beyond ramp-up
+        uf_gf.ProjectBdrCoefficient(u_coeff,u_ess_attr);
+    }
+    tau_gf.ProjectBdrCoefficient(nse.stress_coeff,u_ess_attr);
+
+    /// Fluid solution
+    Array<int> nse_offsets({0,u_fes.GetTrueVSize(), p_fes.GetTrueVSize()});
+    nse_offsets.PartialSum();
+    BlockVector up(nse_offsets);
+
+    uf_gf.GetTrueDofs(up.GetBlock(0));
+    p_gf.GetTrueDofs(up.GetBlock(1));
+
+
+    // Three coupled applications: navier, elasticity, and morpher
+    CoupledOperator multiapp(3);
+    std::unique_ptr<ODESolver> coupled_solver = ODESolver::Select(ode_solver_type);
+
+    Application* nse_app = multiapp.AddOperator(&nse);
+    Application* elasticity_app = multiapp.AddOperator(&elasticity);
+    Application* morpher_app = multiapp.AddOperator(&morpher);
+
+
+    // Set up field transfer
+    // NativeTransfer strsf_to_strss_map(u_fes, xs_fes); // default map if none provided
+    // GSLibTransfer strsf_to_strss_map(u_fes, xs_fes, xs_nat_attr);
+    nse_app->Fields().AddTargetField("Stress", &stress_gf);
+    elasticity_app->Fields().AddTargetField("Velocity_BC", &uf_gf_bc);
+    elasticity_app->Fields().AddTargetField("Displacement_BC", &xf_gf_bc);
+    elasticity_app->Fields().AddTargetField("Velocity_BC", &dxf_gf_bc);
+
+
+    // Assemble the true (offseted) dofs array for the coupled FSI system
+    Array<int> tdofs;
+    std::vector<Array<int>*> state_dofs = {&nse.u_ess_tdofs, &nse.p_ess_tdofs,
+                                           &elasticity.ess_tdofs, &elasticity.ess_tdofs,
+                                           &morpher.ess_tdofs};
+    CollectArrays(state_dofs, tdofs);
+
+
+    // Set up coupling parameters (schemes and solvers)
+    FPISolver fp_solver(MPI_COMM_WORLD); // For partitioned solves
+    AitkenRelaxation fp_relax;
+    ConstrainedInnerProduct constr_ipo(MPI_COMM_WORLD , tdofs);
+    SetSolverParameters(&fp_solver, 0.0, 5e-4, 100, 1, false);
+    fp_solver.SetRelaxation(relax_factor, &fp_relax);
+    // fp_solver.SetRelaxation(relax_factor, nullptr); // Use default relaxation method
+    // fp_relax.SetBounds(-1.0e0,1.0e0);
+    fp_relax.SetAbsoluteLowerBound(1e-1);
+    fp_solver.SetInnerProduct(&constr_ipo);
+    fp_relax.SetInnerProduct(&constr_ipo);
+
+
+    NewtonSolver newton_solver(MPI_COMM_WORLD); // For fully coupled
+    GMRESSolver gmres_solver(MPI_COMM_WORLD); // For fully coupled
+
+    if(couple_scheme == -1)
+    {
+        MFEM_ABORT("Monolithic coupling not supported for FSI.")
+        multiapp.SetCouplingScheme(CoupledOperator::Scheme::MONOLITHIC);
+
+        SetSolverParameters(&gmres_solver, 1e-3, 1e-3, 500, 1, false);
+        gmres_solver.SetKDim(300);
+
+        SetSolverParameters(&newton_solver, 0.0, 1e-4, 30, 1, false);
+        newton_solver.SetSolver(gmres_solver);
+
+        multiapp.SetSolver(&newton_solver); // Set the solver for the multiapp
+    }
+    else if(couple_scheme == 0)
+    {
+        multiapp.SetCouplingScheme(CoupledOperator::Scheme::ADDITIVE_SCHWARZ);
+        multiapp.SetSolver(&fp_solver); // Set the solver for the multiapp
+    }
+    else
+    {
+        multiapp.SetCouplingScheme(CoupledOperator::Scheme::ALTERNATING_SCHWARZ);
+        multiapp.SetSolver(&fp_solver); // Set the solver for the multiapp
+    }
+
+    multiapp.Assemble(false); // Assemble the multiapp (build OperatorCoupler)
+    multiapp.Finalize(false); // Finalize the multiapp (perform checks)
+
+    coupled_solver->Init(multiapp);
+
+    // Set up the initial conditions in block vector for the
+    // coupled application in the correct order
+    int fl_id = nse_app->GetOperatorIndex();
+    int el_id = elasticity_app->GetOperatorIndex();
+    int morph_id = morpher_app->GetOperatorIndex();
+
+    BlockVector xb(multiapp.GetBlockOffsets());
+    xb.GetBlock(el_id) = xu; // Solid displacement and velocity
+    xb.GetBlock(fl_id) = up; // Fluid velocity and pressure
+    xb.GetBlock(morph_id) = xf; // Fluid displacement
+
+
+
+    auto update_grid_functions  = [&](Vector &x) mutable
+    {
+        BlockVector xb(x.GetData(), multiapp.GetBlockOffsets());
+        BlockVector elas_x(xb.GetBlock(el_id).GetData(), solid_offsets);
+        BlockVector nse_x(xb.GetBlock(fl_id).GetData(), nse_offsets);
+        Vector morph_x(xb.GetBlock(morph_id).GetData(), xb.BlockSize(morph_id));
+
+        xf_gf.SetFromTrueDofs(morph_x);
+        us_gf.SetFromTrueDofs(elas_x.GetBlock(0));
+        xs_gf.SetFromTrueDofs(elas_x.GetBlock(1));
+        uf_gf.SetFromTrueDofs(nse_x.GetBlock(0));
+        p_gf.SetFromTrueDofs(nse_x.GetBlock(1));
+        tau_gf.ProjectBdrCoefficient(nse.stress_coeff,u_ess_attr);
+    };
+
+    auto update_nse_grid_functions  = [&](Vector &x) mutable
+    {
+        BlockVector nse_x(x.GetData(), nse_offsets);
+        uf_gf.SetFromTrueDofs(nse_x.GetBlock(0));
+        p_gf.SetFromTrueDofs(nse_x.GetBlock(1));
+        tau_gf.ProjectBdrCoefficient(nse.stress_coeff,u_ess_attr);
+    };    
+
+    // Set up visualization
+    ParaViewDataCollection *fluid_pv = nullptr;
+    ParaViewDataCollection *solid_pv = nullptr;
+
+    if(visualization)
+    {
+        fluid_pv = new ParaViewDataCollection("fsi-fluid", &fluid_mesh);
+        solid_pv = new ParaViewDataCollection("fsi-solid", &solid_mesh);
+
+        fluid_pv->SetLevelsOfDetail(order);
+        fluid_pv->SetDataFormat(VTKFormat::BINARY);
+        fluid_pv->SetHighOrderOutput(true);
+        fluid_pv->RegisterField("displacement",&mesh_disp);
+        fluid_pv->RegisterField("dxdt",&mesh_vel);
+
+        fluid_pv->RegisterField("pressure",&p_gf);
+        fluid_pv->RegisterField("velocity",&uf_gf);
+        fluid_pv->RegisterField("stress",&tau_gf);
+        fluid_pv->RegisterField("ale_velocity",&dxf_gf);
+
+        solid_pv->SetLevelsOfDetail(order);
+        solid_pv->SetDataFormat(VTKFormat::BINARY);
+        solid_pv->SetHighOrderOutput(true);
+        solid_pv->RegisterField("displacement",&xs_gf);
+        solid_pv->RegisterField("velocity",&us_gf);
+        solid_pv->RegisterField("stress",&stress_gf);
+    }
+
+    auto save_callback = [&](int cycle, double t)
+    {
+        if(fluid_pv)
+        {
+            fluid_pv->SetCycle(cycle);
+            fluid_pv->SetTime(t);
+            fluid_pv->Save();
+        }
+
+        if(solid_pv)
+        {
+            solid_pv->SetCycle(cycle);
+            solid_pv->SetTime(t);
+            solid_pv->Save();
+        }
+    };
+
+
+    if (Mpi::Root()) out << "Starting time integration..." << std::endl;
+
+    StopWatch timer;
+    timer.Start();
+
+    real_t t = 0.0;
+    bool last_step = false;
+    int tindex = 1;
+    save_callback(0, t);
+
+
+    // Solve the Navier-Stokes equations to fully developed flow time, t_dev
+    if(t_dev > 0.0 && !init)
+    {
+        nse_app->SetOperationID(Application::OperationID::STEP);
+        for (; !last_step; tindex++)
+        {
+            if (t + dt >= t_dev - dt/2){ last_step = true; }
+
+            u_coeff.SetTime(t); // Slowly ramp-up inlet velocity
+            uf_gf.ProjectBdrCoefficient(u_coeff,u_ess_attr);
+            nse_solver->Step(up,t,dt);
+
+            if (last_step || (tindex % vis_steps) == 0){
+                if (Mpi::Root()) { out << "step " << tindex << ", t = " << t << std::endl;}
+                update_nse_grid_functions(up);
+                save_callback(tindex, t);
+            }
+        }
+
+        if ((myid==0) && lsave)
+        {
+            std::filesystem::path dir_path = init_dir;
+            if(!std::filesystem::is_directory(dir_path)) std::filesystem::create_directory(dir_path);
+        }
+        if(lsave){
+            p_gf.Save((init_dir+"/p-init.gf").c_str());
+            uf_gf.Save((init_dir+"/u-init.gf").c_str());
+        }
+
+        xb.GetBlock(fl_id) = up; // Update nse block
+    }
+
+    vis_steps = 1;
+    last_step = false;
+
+    nse_app->SetCoupled(true);
+    morpher_app->SetCoupled(true);
+    elasticity_app->SetCoupled(true);
+    multiapp.Transfer(xb);
+
+    // Store original fluid nodes; morphing is done w.r.t. original configuration
+    GridFunction fluid_nodes_orig = *(fluid_mesh.GetNodes());
+
+    for (; !last_step; tindex++)
+    {
+        if (t + dt >= t_final - dt/2){ last_step = true; }
+
+        coupled_solver->Step(xb, t, dt);
+        // multiapp.Transfer(xb);
+        update_grid_functions(xb);
+
+
+        // Morph fluid mesh
+        mesh_disp.SetFromTrueDofs(xb.GetBlock(morph_id));
+        mesh_vel = dxf_gf;
+
+        GridFunction *fluid_nodes = fluid_mesh.GetNodes();
+        *fluid_nodes = fluid_nodes_orig;
+        *fluid_nodes += mesh_disp;
+        fluid_mesh.DeleteGeometricFactors();
+
+        // Update FE spaces, grid functions and forms
+        // after fluid mesh update
+        nse_app->Update();
+        morpher_app->Update();
+
+        if (last_step || (tindex % vis_steps) == 0){
+            if (Mpi::Root()){ out << "step " << tindex << ", t = " << t << std::endl;}
+            save_callback(tindex, t);
+        }
+    }
+
+    timer.Stop();
+    if (Mpi::Root()){
+        out << "Total time: " << timer.RealTime() << " seconds." << std::endl;
+    }
+
+    return 0;
+}
+
+void SetSolverParameters(IterativeSolver *solver, real_t rtol, real_t atol , int max_it, 
+                         int print_level, bool iterative_mode)
+{
+    solver->SetRelTol(rtol);
+    solver->SetAbsTol(atol);
+    solver->SetMaxIter(max_it);
+    solver->SetPrintLevel(print_level);
+    solver->iterative_mode = iterative_mode;
+}
+
+bool ReadGridFunctionFromFile(const string &dirname, const string &gf_name,
+                              ParMesh &mesh, ParGridFunction &gf)
+{
+    int myid = Mpi::WorldRank();
+    std::string mpirank = std::to_string(myid);
+    std::string filename = dirname+"/"+gf_name+"."+mpirank.insert(0,6-mpirank.length(),'0');
+    bool sucess = false;
+    if (std::filesystem::exists(filename))
+    {
+        istream *ifile;
+        ifile = new ifstream(filename);
+        gf = ParGridFunction(&mesh,*ifile);
+        delete ifile;
+        sucess = true;
+    }
+
+    return sucess;
+}
+
+
+void CollectArrays(std::vector<Array<int>*> &dof_arrays, Array<int> &tdof_array)
+{
+    int total_dofs = 0;
+    for (const auto& arr : dof_arrays)
+    {
+        total_dofs += arr->Size();
+    }
+
+    int n = 0;
+    tdof_array.SetSize(total_dofs);
+    for (auto arr : dof_arrays)
+    {
+        std::transform(arr->begin(), arr->end(), tdof_array.begin() + n,
+                       [&](int i) { return i+n; });
+        n += arr->Size();
+    }
+}
@@ -0,0 +1,538 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "multiapp.hpp"
+
+#ifdef MFEM_USE_MPI
+
+namespace mfem
+{
+
+FieldTransfer* FieldTransfer::Select(ParFiniteElementSpace *src,
+                                     ParFiniteElementSpace *tar, 
+                                     Type type)
+{
+    switch (type)
+    {
+        case Type::NATIVE:
+            return new NativeTransfer(src, tar);
+        case Type::GSLIB:
+            return new GSLibTransfer(src, tar);
+        default:
+            MFEM_ABORT("Unknown FieldTransfer scheme: " << static_cast<int>(type));
+    }
+}
+
+OperatorCoupler* OperatorCoupler::Select(CoupledOperator *op, 
+                                         Scheme scheme)
+{
+    switch (scheme)
+    {
+        case Scheme::MONOLITHIC:
+        {
+            real_t fd_eps = 1e-6;
+            return new JacobianFreeFullCoupler(op, fd_eps);
+        }
+        case Scheme::ADDITIVE_SCHWARZ:
+            return new AdditiveSchwarzCoupler(op);
+        case Scheme::ALTERNATING_SCHWARZ:
+            return new AlternatingSchwarzCoupler(op);
+        case Scheme::NONE:
+            return nullptr;
+        default:
+            MFEM_ABORT("Unknown coupling scheme: " << static_cast<int>(scheme));
+    }
+}
+
+CoupledOperator::~CoupledOperator()
+{
+    if(solver && own_solver) delete solver;
+    if(op_coupler && own_op_coupler) delete op_coupler;
+
+    for(int i=0; i < nops; i++)
+    {
+        if(operators_owned[i] && operators[i]) delete operators[i];
+    }
+}
+
+void CoupledOperator::SetOperatorCoupler(OperatorCoupler* op, bool own)
+{
+    if(op_coupler && own_op_coupler) delete op_coupler;
+    op_coupler = op;
+    own_op_coupler = own;
+    coupler_type = op_coupler->GetType();
+}
+
+void CoupledOperator::Initialize(bool do_ops)
+{
+    if (do_ops)
+    {
+        for (auto &op : operators)
+        {
+            op->Initialize();
+        }
+    }
+}
+
+void CoupledOperator::Assemble(bool do_ops)
+{
+    if (do_ops)
+    {
+        for (auto &op : operators)
+        {
+            op->Assemble();
+        }
+    }
+
+    // Check block offsets against operator size
+    Array<int> true_offsets(Size()+1);
+    bool offset_consistent = true;
+    true_offsets = 0;
+    int max_size = 0;
+
+    for (int i=0; i < nops; i++)
+    {
+        auto op = GetOperator(i);
+        int block_size = offsets[i+1]-offsets[i];
+        true_offsets[i+1] = true_offsets[i] + op->Width();
+        if (block_size != op->Width())
+        {
+            offset_consistent = false;
+        }
+    }
+    if (!offset_consistent)
+    {
+        MFEM_WARNING("Block offsets inconsistent with operator sizes."
+                        "Using default offsets.");
+        offsets = true_offsets;
+        max_op_size = max_size;
+    }
+
+    if(op_coupler && own_op_coupler) delete op_coupler;
+    op_coupler = OperatorCoupler::Select(this, coupler_type);
+    if(solver) solver->SetOperator(*op_coupler);
+}
+
+void CoupledOperator::Finalize(bool do_ops)
+{
+    if (do_ops)
+    {
+        for (auto &op : operators)
+        {
+            op->Finalize();
+        }
+    }
+}
+
+void CoupledOperator::PreProcess(Vector &x, bool do_ops) 
+{
+    if (do_ops)
+    {
+        BlockVector xb(x.GetData(), offsets);
+        for (int i=0; i < nops; i++)
+        {
+            Vector &xi = xb.GetBlock(i);
+            operators[i]->PreProcess(xi);
+        }
+    }
+}
+
+void CoupledOperator::PostProcess(Vector &x, bool do_ops) 
+{
+    if (do_ops)
+    {
+        BlockVector xb(x.GetData(), offsets);
+        for (int i=0; i < nops; i++)
+        {
+            Vector &xi = xb.GetBlock(i);
+            operators[i]->PostProcess(xi);
+        }
+    }
+}
+
+void CoupledOperator::SetOperationID(OperationID id, bool do_ops)
+{
+    Application::SetOperationID(id);
+    if (do_ops)
+    {
+        for (auto &op : operators)
+        {
+            op->SetOperationID(id);
+        }
+    }
+}
+
+void CoupledOperator::SetTime(const real_t t_) 
+{
+    TimeDependentOperator::SetTime(t_);
+    if(op_coupler) op_coupler->SetTime(t_);
+    for (auto &op : operators)
+    {
+        op->SetTime(t_);
+    }
+}
+
+void CoupledOperator::Transfer(const Vector &x)
+{
+    BlockVector xb(x.GetData(), offsets);
+    for (int i=0; i < nops; i++)
+    {
+        operators[i]->Transfer(xb.GetBlock(i));
+    }
+}
+
+void CoupledOperator::Transfer(const Vector &u, const Vector &k, real_t dt)
+{
+    BlockVector ub(u.GetData(), offsets);
+    BlockVector kb(k.GetData(), offsets);
+    for (int i=0; i < nops; i++)
+    {
+        operators[i]->Transfer(ub.GetBlock(i), kb.GetBlock(i), dt);
+    }
+}
+
+void CoupledOperator::Mult(const Vector &x, Vector &y) const
+{
+    if(op_coupler && coupler_type != Scheme::NONE)
+    {
+        if(solver) {
+            op_coupler->SetOperationID(OperationID::MULT);
+            op_coupler->SetInput(&x);
+            solver->Mult(b,y);
+        }
+        else {
+            op_coupler->Mult(x,y);
+        }
+    }
+    else
+    {
+        BlockVector xb(x.GetData(), offsets);
+        BlockVector yb(y.GetData(), offsets);
+
+        for (int i=0; i < nops; i++)
+        {
+            operators[i]->SetOperationID(OperationID::MULT);
+            operators[i]->Mult(xb.GetBlock(i), yb.GetBlock(i));
+        }
+    }
+}
+
+void CoupledOperator::ImplicitSolve(const real_t dt, const Vector &x, Vector &k ){
+
+    if(op_coupler && coupler_type != Scheme::NONE)
+    {
+        if(solver) {
+            op_coupler->SetOperationID(OperationID::IMPLICIT_SOLVE); ///< OperatorCoupler::Mult() -> OperatorCoupler::ImplicitSolve() 
+            op_coupler->SetTimeStep(dt);
+            op_coupler->SetInput(&x);
+            solver->Mult(b,k);
+        }
+        else {
+            op_coupler->ImplicitSolve(dt,x,k);
+        }                
+    }
+    else
+    {
+        BlockVector xb(x.GetData(), offsets);
+        BlockVector kb(k.GetData(), offsets);
+
+        for (int i=0; i < nops; i++)
+        {
+            Vector &xi = xb.GetBlock(i);
+            Vector &ki = kb.GetBlock(i);
+            operators[i]->SetOperationID(OperationID::IMPLICIT_SOLVE);
+            operators[i]->ImplicitSolve(dt,xi,ki); ///< Solve the implicit system for the application
+        }
+    }
+}
+
+void CoupledOperator::Step(Vector &x, real_t &t, real_t &dt)
+{
+    if(op_coupler && coupler_type != Scheme::NONE)
+    {
+        if(solver) {
+            op_coupler->SetOperationID(OperationID::STEP); ///< OperatorCoupler::Mult() -> OperatorCoupler::Mult() 
+            op_coupler->SetTimeStep(dt); ///< Set the time step for the ODE Solver
+            op_coupler->SetTime(t);
+            op_coupler->SetInput(&x);
+            solver->Mult(b,x);
+        }
+        else {
+            op_coupler->Step(x,t,dt);
+        }                
+    }
+    else
+    {
+        BlockVector xb(x.GetData(), offsets);
+        for (int i=0; i < nops; i++)
+        {
+            real_t t0 = t;   ///< Store the current time
+            real_t dt0 = dt; ///< Store the current time step
+            Vector &xi = xb.GetBlock(i);
+            operators[i]->SetOperationID(OperationID::STEP);
+            operators[i]->Step(xi,t0,dt0); ///< Advance the time step for application
+        }
+        t += dt; ///< Update the time after all applications have been stepped forward
+                 ///< NOTE: does not work for adaptive time-stepping
+    }
+}
+
+void CoupledOperator::ImplicitMult(const Vector &u, const Vector &k, Vector &v) const 
+{
+    BlockVector ub(u.GetData(), offsets);
+    BlockVector kb(k.GetData(), offsets);
+    BlockVector vb(v.GetData(), offsets);
+
+    for (int i=0; i < nops; i++)
+    {
+        Vector &ui = ub.GetBlock(i);
+        Vector &ki = kb.GetBlock(i);
+        Vector &vi = vb.GetBlock(i);
+        operators[i]->SetOperationID(OperationID::IMPLICIT_MULT);
+        operators[i]->ImplicitMult(ui,ki,vi); ///< Solve the implicit system for the application
+    }
+}
+
+void CoupledOperator::ExplicitMult(const Vector &u, Vector &v) const
+{
+    BlockVector ub(u.GetData(), offsets);
+    BlockVector vb(v.GetData(), offsets);
+
+    for (int i=0; i < nops; i++)
+    {
+        Vector &ui = ub.GetBlock(i);
+        Vector &vi = vb.GetBlock(i);
+        operators[i]->SetOperationID(OperationID::EXPLICIT_MULT);
+        operators[i]->ExplicitMult(ui,vi); ///< Solve the implicit system for the application
+    }
+}
+
+
+
+
+// AdditiveSchwarzCoupler methods
+void AdditiveSchwarzCoupler::Mult(const Vector &x, Vector &y) const
+{
+    /// This is use to call either ImplicitSolve or Step when Solver::Mult()
+    /// calls Solver.Operator::Mult()
+    if(GetOperationID() == OperationID::IMPLICIT_SOLVE)
+    {
+        y=x; // input vector passed as initial guess for k in ImpliicitSolve
+        ImplicitSolve(timestep,*input,y);
+        return;
+    }
+    else if(GetOperationID() == OperationID::STEP)
+    {
+        y=x; // input vector passed as initial condition in Step
+        real_t t_ = t, dt = timestep; 
+        Step(y,t_,dt);
+        return;
+    }    
+
+    int nops = coupled_op->Size();
+    const Array<int> offsets = coupled_op->GetBlockOffsets();
+
+    BlockVector xb(x.GetData(), offsets);
+    BlockVector yb(y.GetData(), offsets);
+
+    for (int i=0; i < nops; i++)
+    {
+        Vector &xi = xb.GetBlock(i);
+        Vector &yi = yb.GetBlock(i);
+        auto op = coupled_op->GetOperator(i);
+        op->Transfer(xi,yi,0.0);
+    }
+
+    for (int i=0; i < nops; i++)
+    {
+        Vector &xi = xb.GetBlock(i);
+        Vector &yi = yb.GetBlock(i);
+        auto op = coupled_op->GetOperator(i);
+        op->SetOperationID(OperationID::MULT);
+        
+        op->PreProcess(xi); ///< Postprocess the data for the application
+        op->Mult(xi,yi);
+        op->PostProcess(yi); ///< Postprocess the data for the application
+    }
+}
+
+void AdditiveSchwarzCoupler::ImplicitSolve(const real_t dt, const Vector &x, Vector &k ) const
+{
+    int nops = coupled_op->Size();
+    const Array<int> offsets = coupled_op->GetBlockOffsets();
+
+    BlockVector xb(x.GetData(), offsets);
+    BlockVector kb(k.GetData(), offsets);
+
+    for (int i=0; i < nops; i++)
+    {
+        Vector &xi = xb.GetBlock(i);
+        Vector &ki = kb.GetBlock(i);
+        auto op = coupled_op->GetOperator(i);
+        op->Transfer(xi,ki,dt);
+    }
+
+    for (int i=0; i < nops; i++)
+    {
+        Vector &xi = xb.GetBlock(i);
+        Vector &ki = kb.GetBlock(i);
+
+        auto op = coupled_op->GetOperator(i);
+        op->SetOperationID(OperationID::IMPLICIT_SOLVE);
+
+        op->PreProcess(xi);
+        op->ImplicitSolve(dt,xi,ki);
+        op->PostProcess(ki);
+    }
+}
+
+void AdditiveSchwarzCoupler::Step(Vector &x, real_t &t_, real_t &dt) const
+{
+    int nops = coupled_op->Size();
+    const Array<int> offsets = coupled_op->GetBlockOffsets();
+
+    BlockVector xb(x.GetData(), offsets);
+
+    for (int i=0; i < nops; i++)
+    {
+        Vector &xi = xb.GetBlock(i);
+        auto op = coupled_op->GetOperator(i);
+        op->Transfer(xi);
+    }
+
+    // TODO: Add time-interpolation to enable different time step for each operator; 
+    //       currently, all operators are stepped forward with the same time step
+    for (int i=0; i < nops; i++)
+    {
+        real_t ti = t_;  ///< Store the current time
+        real_t dti = dt; ///< Store the current time step
+
+        Vector &xi = xb.GetBlock(i);
+        auto op = coupled_op->GetOperator(i);
+        op->SetOperationID(OperationID::STEP);
+        
+        op->PreProcess(xi);
+        op->Step(xi,ti,dti);
+        op->PostProcess(xi);
+    }
+
+    t_ += dt; ///< Update the time after all applications have been stepped forward
+               ///< NOTE: does not work for adaptive time-stepping
+}
+
+
+// AlternatingSchwarzCoupler methods
+void AlternatingSchwarzCoupler::Mult(const Vector &x, Vector &y) const
+{
+    /// This is use to call either ImplicitSolve or Step when Solver::Mult()
+    /// calls Solver.Operator::Mult()
+    if(GetOperationID() == OperationID::IMPLICIT_SOLVE)
+    {
+        y=x; // input vector passed as initial guess for k in ImpliicitSolve
+        ImplicitSolve(timestep,*input,y);
+        return;
+    }
+    else if(GetOperationID() == OperationID::STEP)
+    {
+        y=x; // input vector passed as initial condition in Step
+        real_t t_ = t, dt = timestep; 
+        Step(y,t_,dt);
+        return;
+    }
+
+    int nops = coupled_op->Size();
+    const Array<int> offsets = coupled_op->GetBlockOffsets();
+
+    BlockVector xb(x.GetData(), offsets);
+    BlockVector yb(y.GetData(), offsets);
+
+    for (int i=0; i < nops; i++)
+    {
+        Vector &xi = xb.GetBlock(i);
+        Vector &yi = yb.GetBlock(i);
+        auto op = coupled_op->GetOperator(i);
+        op->SetOperationID(OperationID::MULT);
+
+        op->PreProcess(xi);
+        op->Mult(xi,yi);
+        op->PostProcess(yi);
+        op->Transfer(xi,yi,0.0);
+    }
+}
+
+void AlternatingSchwarzCoupler::ImplicitSolve(const real_t dt, const Vector &x, Vector &k ) const
+{
+    int nops = coupled_op->Size();
+    const Array<int> offsets = coupled_op->GetBlockOffsets();
+
+    BlockVector xb(x.GetData(), offsets);
+    BlockVector kb(k.GetData(), offsets);
+
+    for (int i=0; i < nops; i++)
+    {
+        Vector &xi = xb.GetBlock(i);
+        Vector &ki = kb.GetBlock(i);
+
+        auto op = coupled_op->GetOperator(i);
+        op->SetOperationID(OperationID::IMPLICIT_SOLVE);
+
+        op->PreProcess(xi);
+        op->ImplicitSolve(dt,xi,ki);
+        op->PostProcess(ki);
+        op->Transfer(xi,ki,dt);
+    }
+}
+
+void AlternatingSchwarzCoupler::Step(Vector &x, real_t &t_, real_t &dt) const
+{
+    int nops = coupled_op->Size();
+    const Array<int> offsets = coupled_op->GetBlockOffsets();
+
+    BlockVector xb(x.GetData(), offsets);
+
+    for (int i=0; i < nops; i++)
+    {
+        real_t ti = t_;   ///< Store the current time
+        real_t dti = dt;  ///< Store the current time step
+
+        Vector &xi = xb.GetBlock(i);
+        auto op = coupled_op->GetOperator(i);
+        op->SetOperationID(OperationID::STEP);
+
+        op->PreProcess(xi);
+        op->Step(xi,ti,dti);
+        op->PostProcess(xi);
+        op->Transfer(xi);
+    }
+
+    t_ += dt; ///< Update the time after all applications have been stepped forward
+               ///< NOTE: does not work for adaptive time-stepping
+}
+
+
+
+// JacobianFreeFullCoupler methods
+void JacobianFreeFullCoupler::Mult(const Vector &k, Vector &y) const
+{
+    add(1.0,*input,timestep,k,u); // u = u + dt*k
+    coupled_op->Transfer(u);
+    coupled_op->ImplicitMult(u,k,y); //compute residual y = f(u,k,t)
+}
+
+Operator& JacobianFreeFullCoupler::GetGradient(const Vector &k) const
+{
+    grad.Update(k);
+    return const_cast<future::FDJacobian&>(grad);
+}
+
+}
+
+#endif // MFEM_USE_MPI
Author	SHA1	Message	Date
Sohail Reddy	b7753cb4db	fixed cmake	2026-01-06 18:59:44 -08:00
Sohail Reddy	842fec483f	Merge branch 'master' into multiapp	2026-01-06 18:53:05 -08:00
Sohail Reddy	4b0eafd047	Fixed const correctness	2026-01-06 18:49:38 -08:00
Sohail Reddy	e793c70058	Updated fsi app for constrained inner product on ess_dof	2025-11-19 14:28:38 -08:00
Sohail Reddy	1856912d63	Fixed style	2025-11-19 14:24:23 -08:00
Sohail Reddy	3b6817e36b	Added custom inner product to fixed point relaxation methods	2025-11-19 14:24:04 -08:00
Sohail Reddy	e4dcfb684f	Merge branch 'ode-stage-solve' into multiapp	2025-11-18 23:26:55 -08:00
Sohail Reddy	7dcf77da54	Merge branch 'master' into multiapp	2025-11-18 23:25:23 -08:00
Sohail Reddy	3e823f1da7	updated morpher coupled bc	2025-11-18 23:16:06 -08:00
Sohail Reddy	9a2fee8b6f	moved FPI relaxation member functions to cpp	2025-10-19 21:44:25 -07:00
Sohail Reddy	4b38bd79e9	clean up fsi example	2025-10-16 18:05:41 -07:00
Sohail Reddy	d1378b8a7b	bug fix fpisolver; reset relaxation each call	2025-10-16 18:04:17 -07:00
Sohail Reddy	3972b24f3b	fixed ale velocity calculation	2025-10-15 22:19:59 -07:00
Sohail Reddy	48bd66640d	added command line option for relaxation factor and fixed mpi_c_bool	2025-10-15 22:19:35 -07:00
Sohail Reddy	a230b87a20	Fixed pure virtual function overide if GSLib not built	2025-10-14 18:31:08 -07:00
Sohail Reddy	5f9f782ad2	fixed style	2025-10-14 17:01:51 -07:00
Sohail Reddy	159a259f30	fixed redeclared variable in FPISolver	2025-10-14 16:56:20 -07:00
Sohail Reddy	67da324c26	Merge branch 'master' into multiapp	2025-10-14 16:28:08 -07:00
Sohail Reddy	b5cf20e76c	added cmake for multiapp miniapps	2025-10-14 16:26:00 -07:00
Sohail Reddy	1510492dfd	fixed style	2025-10-14 16:24:23 -07:00
Sohail Reddy	360e662cf8	FSI miniapp and updated elasticity and navier-stokes apps	2025-10-14 16:22:53 -07:00
Sohail Reddy	c1ec6a8c72	bug fix in navier miniapp for Step() function	2025-10-14 13:27:04 -07:00
Sohail Reddy	805499d799	bugfix for updating grid functions for transfer	2025-10-14 13:26:15 -07:00
Sohail Reddy	7b0330ffca	Removed repetitive code in cht miniapp	2025-10-14 13:23:55 -07:00
Sohail Reddy	68fd37de1d	Added completed conjugate heat transfer test case	2025-10-11 17:44:41 -07:00
Sohail Reddy	9f1d2ac249	reverted changes to navier miniapp	2025-10-11 17:43:24 -07:00
Sohail Reddy	a594ab5a3a	fixed namespace std use	2025-10-11 17:42:57 -07:00
Sohail Reddy	343a8b473e	moved fixed point solver to solvers.cpp	2025-10-11 17:13:31 -07:00
Sohail Reddy	911709c619	fixed use of MFEM_USE_MPI in multiapp	2025-10-11 16:10:17 -07:00
Sohail Reddy	c303956b9f	fixed source field ownership and cleaning	2025-10-11 01:29:50 -07:00
Sohail Reddy	c031ec4466	Fixed const correctness and added operator ownership in multiapp	2025-10-10 17:37:12 -07:00
Sohail Reddy	58c56ec9dc	added LinkedFieldsCollection to manage linked fields for each application	2025-10-10 00:01:12 -07:00
Sohail Reddy	7db740ede3	Added documentation, updated class names for coupling schemes	2025-09-28 23:36:06 -07:00
Sohail Reddy	0436adfb1c	updated fixed point solver	2025-09-28 23:31:33 -07:00
Sohail Reddy	a6f2d8981a	Merge branch 'master' into multiapp	2025-09-08 18:02:54 -07:00
Sohail Reddy	90edc198e3	update fsi app	2025-09-08 17:58:24 -07:00
Sohail Reddy	c691c7c854	Formatted multiapp files and added documentation	2025-09-08 10:21:16 -07:00
Sohail Reddy	7c9105bd0e	fixed style	2025-09-08 01:41:26 -07:00
Sohail Reddy	773bac0114	updated navier-stokes miniapp for ALE	2025-09-08 01:29:33 -07:00
Sohail Reddy	440db8189c	Added option for ALE to VectorConvectionNLFIntegrator	2025-09-08 01:28:26 -07:00
Sohail Reddy	94efb4cd81	Updated multiapp interface for type-erased operator, field transfer and coupling scheme	2025-09-08 01:27:47 -07:00
Sohail Reddy	ce3bdc0dfc	Simplified mesh-morphing/diffusion app	2025-09-08 01:24:53 -07:00
Sohail Reddy	aaf59d3a83	updated elasticity miniapp to use Schur complement	2025-09-08 01:22:39 -07:00
Sohail Reddy	3b3ffa7d7a	updated native navier miniapp to use standard ODESolve::Step signature	2025-09-05 13:53:12 -07:00
Sohail Reddy	2898ec669c	Updated native navier miniapp to own/not-own the boundary conditions	2025-09-04 17:19:45 -07:00
Sohail Reddy	f34bfbbd8d	Added field transfer based on FindPoints	2025-09-04 10:35:50 -07:00
Sohail Reddy	3cec2cba78	updated navier-stokes and elasticity examples	2025-09-03 16:34:45 -07:00
Sohail Reddy	f769e2a7d7	updated morpher to diffuse displacement or velocity	2025-09-03 16:34:08 -07:00
Sohail Reddy	493674e650	updated navier-stokes miniapp to work with scaled pressure	2025-09-03 16:31:58 -07:00
Sohail Reddy	13632cde51	updated elasticity miniapp	2025-09-03 16:31:03 -07:00
Sohail Reddy	ae63fda0fd	Updated multiapp object for multi-stage coupling	2025-08-19 11:54:02 -07:00
Sohail Reddy	48b3e428fb	clean up	2025-08-19 11:53:32 -07:00
Sohail Reddy	89dd191866	bug-fix: fixed point solver; updated relaxation schemes to clamp relaxation factors	2025-08-19 11:47:14 -07:00
Sohail Reddy	2a4f9184be	Updated Navier-Stokes miniapp for multi-stage coupling and artificial compressibility	2025-08-19 11:40:47 -07:00
Sohail Reddy	ff66dce92e	updated elasticity and morpher miniapps for multi-stage coupling	2025-08-19 11:38:25 -07:00
Sohail Reddy	01ef886faf	Added parallelized DistanceSquared for Vectors	2025-08-19 11:09:47 -07:00
Sohail Reddy	4be67907c4	Working FSI miniapp	2025-08-13 16:12:04 -07:00
Sohail Reddy	9443523dec	Updated the multiapp (CoupledApplication) framework	2025-08-13 16:10:19 -07:00
Sohail Reddy	fac6dfe283	Added mesh morphing app	2025-08-13 16:09:50 -07:00
Sohail Reddy	7776766ee7	updated FPISolver	2025-08-13 16:09:09 -07:00
Sohail Reddy	8e24224a18	updated cht example...unused atm	2025-08-13 16:08:33 -07:00
Sohail Reddy	fa4c607622	updated Navier-Stokes miniapp for partial assembly	2025-08-13 16:07:09 -07:00
Sohail Reddy	29bd541a3e	updated elasticity miniapp	2025-08-13 16:03:23 -07:00
Sohail Reddy	369f82adad	Added function to fetch TimeDependantOperator from ODESolver	2025-08-13 16:01:39 -07:00
Sohail Reddy	e2bc1aefd8	Working elasticity miniapp	2025-08-01 19:09:05 -07:00
Sohail Reddy	033814fe70	fixed code style	2025-07-31 10:14:48 -07:00
Sohail Reddy	f4d4a5c4da	Working fully implicit incompressible navier-stokes solver	2025-07-31 10:13:33 -07:00
Sohail Reddy	6b8d341478	Renamed some multiapp member functions	2025-07-31 09:56:27 -07:00
Sohail Reddy	7a02905d32	Updated constructors for some bilinear and nonlinearform integrator to now take integration rule as input	2025-07-31 09:53:16 -07:00
Sohail Reddy	a4c2ad2f3c	Updated type erasure to support non-mfem apps	2025-07-21 10:27:35 -07:00
Sohail Reddy	4a11684323	Working monolithic (linear and nonlinear) solves for multiapp	2025-07-17 18:22:41 -07:00
Sohail Reddy	7b4b2de092	Merge branch 'master' into multiapp	2025-07-15 16:37:35 -07:00
Sohail Reddy	285aa03049	conv-diff mulitapp example	2025-07-15 16:35:56 -07:00
Sohail Reddy	27af36adc5	Changed class naming and fixed point solver and relaxation schemes	2025-07-15 16:35:03 -07:00
Sohail Reddy	de781c05b8	Fixed function pointer for relaxation method in FPISolver	2025-07-10 17:04:58 -07:00
Sohail Reddy	125f5c9e48	Added relaxation methods to fixed point solvers and updated multiapp for explicit coupling	2025-07-10 16:59:03 -07:00
Sohail Reddy	0d6c8fbbf1	Updated multiapp partition implicit solve to use fixed point Solver object	2025-07-09 19:48:04 -07:00
Sohail Reddy	34ebc0f158	Changed scheme naming convention and added/separated Addtive and Alternating Schwarz coupling schemes	2025-07-08 18:16:32 -07:00
Sohail Reddy	8db7813d86	Working paritioned implicit solves	2025-07-08 11:06:27 -07:00
Sohail Reddy	d31e5c8af0	Merge branch 'master' into multiapp	2025-07-02 20:46:57 -07:00
Sohail Reddy	d7e43cfc74	added test case for two-way partitioned coupled solve	2025-07-02 20:45:51 -07:00
Sohail Reddy	82193e0404	Merge branch 'master' into multiapp	2025-06-16 15:09:01 -07:00
Sohail Reddy	5796633307	updated example for implicit time-integration	2025-06-16 15:05:22 -07:00
Sohail Reddy	3d55c021e3	updated header copyright	2025-06-12 14:47:11 -07:00
Sohail Reddy	a30303ecf9	fixed type checking with SFINAE	2025-06-12 14:31:05 -07:00
Sohail Reddy	9287488024	Merge branch 'master' into multiapp	2025-06-11 14:11:32 -07:00
Sohail Reddy	9aa2866b0d	Merge branch 'master' into multiapp	2025-05-29 14:00:03 -07:00
Sohail Reddy	ad16d2f5ae	added field transfer, multiapp abstractions and working one-way coupling	2025-05-29 13:58:58 -07:00
Sohail Reddy	c552925a54	Merge branch 'master' into multiapp	2025-03-13 15:41:44 -07:00
Sohail Reddy	c3b9b62336	updated mfem_application to perform index controlled operators	2025-02-20 11:28:02 -08:00
Sohail Reddy	c36e5facd1	initial multiapp skeleton with MFEMApplication class and type-erased multiapp	2024-10-02 17:11:14 -07:00