Implements partitioned Fluid-Structure Interaction (FSI) coupling in svMultiPhysics

Code/Source/solver/CMakeLists.txt

@@ -178,6 +178,8 @@ set(CSRCS
   eq_assem.h eq_assem.cpp
   fluid.h fluid.cpp
   fsi.h fsi.cpp
+  fsi_coupling.h fsi_coupling.cpp
+  PartitionedFSI.h PartitionedFSI.cpp
   fs.h fs.cpp
   fft.h fft.cpp
   heatf.h heatf.cpp
@@ -357,7 +359,13 @@ if(ENABLE_UNIT_TEST)
 
   # include source files (same as what svMultiPhysics does except for main.cpp)
   add_executable(run_all_unit_tests ${CSRCS})
-
+
+  # Define test data directory for integration tests
+  # CMAKE_SOURCE_DIR points to Code/, so go up one level to reach the repo root
+  target_compile_definitions(run_all_unit_tests PRIVATE
+    TEST_DATA_DIR="${CMAKE_SOURCE_DIR}/../tests/cases"
+  )
+
   if(USE_TRILINOS)
     target_link_libraries(run_all_unit_tests ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES})
   endif()

Code/Source/solver/ComMod.h

@@ -1571,9 +1571,14 @@ class ComMod {
     /// @brief Whether there is a requirement to update mesh and Dn-Do variables
     bool dFlag = false;
 
-    /// @brief Whether mesh is moving
+    /// @brief Whether mesh is moving (monolithic FSI: mesh velocity in DOFs nsd+1..2*nsd)
     bool mvMsh = false;
 
+    /// @brief ALE mesh velocity for partitioned FSI (nsd, tnNo).
+    /// When non-empty, the fluid assembly subtracts this from the convective
+    /// velocity, providing the ALE correction without expanding tDof.
+    Array<double> ale_mesh_velocity;
+
     /// @brief Whether to averaged results
     bool saveAve = false;
 

Code/Source/solver/Integrator.cpp

@@ -17,6 +17,7 @@
 #include "utils.h"
 
 #include <algorithm>
+#include <cmath>
 #include <iostream>
 #include <set>
 
@@ -171,6 +172,94 @@ bool Integrator::step() {
   } // End of Newton iteration loop
 }
 
+//------------------------
+// step_equation
+//------------------------
+/// @brief Solve a single equation to convergence without cycling to other equations
+bool Integrator::step_equation(int iEq, std::function<void()> post_assembly) {
+  using namespace consts;
+
+  auto& com_mod = simulation_->com_mod;
+  auto& cm_mod = simulation_->cm_mod;
+
+  int& cTS = com_mod.cTS;
+
+  // Set up for this equation
+  com_mod.cEq = iEq;
+  auto& eq = com_mod.eq[iEq];
+  eq.itr = 0;
+  eq.ok = false;
+  eq.iNorm = 0.0;
+
+  newton_count_ = 1;
+
+  while (true) {
+    istr_ = "_" + std::to_string(cTS) + "_" + std::to_string(newton_count_);
+    auto& eq = com_mod.eq[iEq];
+
+    // Coupled BC handling (same as step())
+    if (com_mod.cplBC.coupled && iEq == 0) {
+      set_bc::set_bc_cpl(com_mod, cm_mod, solutions_);
+      set_bc::set_bc_dir(com_mod, solutions_);
+    }
+
+    // Initiator step for Generalized α-Method
+    initiator_step();
+
+    if (com_mod.Rd.size() != 0) {
+      com_mod.Rd = 0.0;
+      com_mod.Kd = 0.0;
+    }
+
+    // Allocate, assemble, apply BCs
+    allocate_linear_system(eq);
+    set_body_forces();
+    assemble_equations();
+    apply_boundary_conditions();
+
+    // Optional post-assembly callback (e.g., for injecting interface traction)
+    if (post_assembly) {
+      post_assembly();
+    }
+
+    // Synchronize R across processes
+    if (!eq.assmTLS) {
+      all_fun::commu(com_mod, com_mod.R);
+    }
+
+    // Update residual in displacement equation for USTRUCT phys
+    if (com_mod.sstEq) {
+      ustruct::ustruct_r(com_mod, solutions_);
+    }
+
+    // Set the residual of the continuity equation to 0 on edge nodes
+    if (std::set<EquationType>{Equation_stokes, Equation_fluid, Equation_ustruct, Equation_FSI}.count(eq.phys) != 0) {
+      fs::thood_val_rc(com_mod);
+    }
+
+    set_bc::set_bc_undef_neu(com_mod);
+    update_residual_arrays(eq);
+
+    // Solve linear system
+    solve_linear_system();
+
+    // Update solution and check convergence (no equation cycling)
+    update_solution();
+
+    if (eq.ok) {
+      return true;
+    }
+
+    // Abort on NaN in residual norm (indicates divergence).
+    if (newton_count_ > 1 && std::isnan(eq.FSILS.RI.iNorm)) {
+      return false;
+    }
+
+    output::output_result(simulation_, com_mod.timeP, 2, iEq);
+    newton_count_ += 1;
+  }
+}
+
 //------------------------
 // initiator_step
 //------------------------
@@ -374,8 +463,8 @@ void Integrator::update_residual_arrays(eqType& eq) {
 //  1.  Bazilevs, et al. "Isogeometric fluid-structure interaction:
 //      theory, algorithms, and computations.", Computational Mechanics,
 //      43 (2008): 3-37. doi: 10.1007/s00466-008-0315-x
-//  2. Bazilevs, et al. "Variational multiscale residual-based 
-//      turbulence modeling for large eddy simulation of incompressible 
+//  2. Bazilevs, et al. "Variational multiscale residual-based
+//      turbulence modeling for large eddy simulation of incompressible
 //      flows.", CMAME (2007)
 //------------------------
 // predictor (picp)
@@ -621,11 +710,13 @@ void Integrator::initiator(SolutionStates& solutions)
   }
 }
 //------------------------
-// corrector
+// update_solution
 //------------------------
-/// @brief Corrector with convergence check
+/// @brief Update solution from linear solve result and check convergence
 ///
-/// Decision for next eqn is also made here (modifies cEq global).
+/// Performs the corrector update of An, Yn, Dn from the Newton solve result
+/// and checks convergence norms. Sets eq.ok if converged.
+/// Does NOT handle equation cycling (that is done by corrector()).
 ///
 /// Modifies:
 /// \code {.cpp}
@@ -637,13 +728,12 @@ void Integrator::initiator(SolutionStates& solutions)
 ///   com_mod.pS0
 ///   com_mod.pSa
 ///   com_mod.pSn
-///
-///   com_mod.cEq
 ///   eq.FSILS.RI.iNorm
 ///   eq.pNorm
+///   eq.ok
 /// \endcode
 //
-void Integrator::corrector()
+void Integrator::update_solution()
 {
   using namespace consts;
 
@@ -853,10 +943,42 @@ void Integrator::corrector()
     dmsg << "com_mod.eq[1].ok: " << com_mod.eq[1].ok;
     #endif
   }
+}
+
+//------------------------
+// corrector
+//------------------------
+/// @brief Corrector with convergence check and equation cycling
+///
+/// Calls update_solution() to update the solution and check convergence,
+/// then handles equation switching for coupled problems (modifies cEq global).
+//
+void Integrator::corrector()
+{
+  using namespace consts;
 
+  auto& com_mod = simulation_->com_mod;
+  const int tnNo = com_mod.tnNo;
+  auto& cEq = com_mod.cEq;
+  auto& eq = com_mod.eq[cEq];
+
+  auto& An = solutions_.current.get_acceleration();
+  auto& Dn = solutions_.current.get_displacement();
+  auto& Yn = solutions_.current.get_velocity();
+
+  #define n_debug_corrector_cycling
+  #ifdef debug_corrector_cycling
+  DebugMsg dmsg(__func__, com_mod.cm.idcm());
+  dmsg.banner();
+  #endif
+
+  // Update solution and check convergence
+  update_solution();
+
+  // Check if all equations converged - if so, skip equation cycling
   auto& eqs = com_mod.eq;
   if (std::count_if(eqs.begin(),eqs.end(),[](eqType& eq){return eq.ok;}) == eqs.size()) {
-    #ifdef debug_corrector
+    #ifdef debug_corrector_cycling
     dmsg << "all ok";
     #endif
     return;
@@ -868,7 +990,7 @@ void Integrator::corrector()
       // For coupled equations, if explicit geometric coupling is not used,
       // increment the equation counter after each linear solve
       cEq = cEq + 1;
-      #ifdef debug_corrector
+      #ifdef debug_corrector_cycling
       dmsg << "eq " << " coupled ";
       dmsg << "1st update cEq: " << cEq;
       #endif
@@ -929,7 +1051,7 @@ void Integrator::corrector()
       cEq = cEq + 1;
     }
   }
- #ifdef debug_corrector
+ #ifdef debug_corrector_cycling
  dmsg << "eq " << " coupled ";
  dmsg << "2nd update cEq: " << cEq;
  #endif

Code/Source/solver/Integrator.h

@@ -9,6 +9,38 @@
 #include "Vector.h"
 #include "Simulation.h"
 
+#include <functional>
+
+/// @brief Newmark time integration utilities.
+///
+/// Compute consistent state variables from a prescribed displacement or
+/// velocity using the Newmark-beta / generalized-alpha relationships:
+///   Dn = Do + dt*Yo + dt^2*((0.5-beta)*Ao + beta*An)
+///   Yn = Yo + dt*((1-gamma)*Ao + gamma*An)
+namespace newmark {
+
+/// Compute acceleration and velocity from prescribed displacement.
+inline void state_from_displacement(
+    double d_new, double d_old, double v_old, double a_old,
+    double dt, double beta, double gam,
+    double& a_new, double& v_new)
+{
+  a_new = (d_new - d_old - dt * v_old) / (beta * dt * dt)
+        - (0.5 - beta) / beta * a_old;
+  v_new = v_old + dt * ((1.0 - gam) * a_old + gam * a_new);
+}
+
+/// Compute acceleration from prescribed velocity.
+inline void state_from_velocity(
+    double v_new, double v_old, double a_old,
+    double dt, double gam,
+    double& a_new)
+{
+  a_new = (v_new - v_old) / (gam * dt) - (1.0 - gam) / gam * a_old;
+}
+
+} // namespace newmark
+
 /**
  * @brief Integrator class encapsulates the Newton iteration loop for time integration
  *
@@ -42,6 +74,21 @@ class Integrator {
    */
   bool step();
 
+  /**
+   * @brief Solve a single equation to convergence (for partitioned coupling)
+   *
+   * Runs Newton iterations for only the specified equation, without cycling
+   * to other equations. Used by partitioned FSI to solve fluid, solid, and
+   * mesh equations independently.
+   *
+   * @param iEq Index of the equation to solve
+   * @param post_assembly Optional callback invoked after boundary condition
+   *        application but before the linear solve. Used by partitioned FSI
+   *        to inject interface traction into the residual.
+   * @return True if the equation converged, false if max iterations reached
+   */
+  bool step_equation(int iEq, std::function<void()> post_assembly = nullptr);
+
   /**
    * @brief Perform predictor step for next time step
    *
@@ -173,11 +220,20 @@ class Integrator {
    */
   void initiator(SolutionStates& solutions);
 
+  /**
+   * @brief Update solution and check convergence of current equation
+   *
+   * Performs the corrector step: updates An, Yn, Dn from the linear solve
+   * result and checks convergence norms. Sets eq.ok if converged.
+   * Does NOT handle equation cycling -- that is done separately in corrector().
+   */
+  void update_solution();
+
   /**
    * @brief Corrector function with convergence check (corrector)
    *
-   * Updates solution at n+1 time level and checks convergence of Newton
-   * iterations. Also handles equation switching for coupled problems.
+   * Calls update_solution(), then handles equation switching for coupled
+   * problems (modifies cEq global).
    */
   void corrector();