diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index cee095791..e30cfa28a 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -182,6 +182,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 heatf.h heatf.cpp heats.h heats.cpp @@ -382,7 +384,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() diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index c7f0dc56e..6e661557b 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -1609,9 +1609,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 ale_mesh_velocity; + /// @brief Whether to averaged results bool saveAve = false; diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index 90dddb1d1..dcb936c7d 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -17,6 +17,7 @@ #include "utils.h" #include +#include #include #include @@ -171,6 +172,212 @@ 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 post_assembly, + std::function pre_boundary_conditions, + std::function post_boundary_conditions) { + 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(); + if (pre_boundary_conditions) { + pre_boundary_conditions(); + } + try { + apply_boundary_conditions(); + } catch (...) { + if (post_boundary_conditions) { + post_boundary_conditions(); + } + throw; + } + if (post_boundary_conditions) { + post_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{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; + } +} + +//------------------------ +// assemble_equation_residual +//------------------------ +/// @brief Assemble a single equation residual without solving or correcting. +void Integrator::assemble_equation_residual( + int iEq, std::function post_assembly, + std::function pre_boundary_conditions, + std::function post_boundary_conditions) +{ + using namespace consts; + + auto& com_mod = simulation_->com_mod; + auto& cm_mod = simulation_->cm_mod; + + const int saved_cEq = com_mod.cEq; + const int saved_dof = com_mod.dof; + const int saved_newton_count = newton_count_; + const std::string saved_istr = istr_; + + com_mod.cEq = iEq; + auto& eq = com_mod.eq[iEq]; + + const int saved_itr = eq.itr; + const bool saved_ok = eq.ok; + const double saved_iNorm = eq.iNorm; + const double saved_pNorm = eq.pNorm; + const auto saved_RI = eq.FSILS.RI; + + auto restore_state = [&]() { + eq.itr = saved_itr; + eq.ok = saved_ok; + eq.iNorm = saved_iNorm; + eq.pNorm = saved_pNorm; + eq.FSILS.RI = saved_RI; + com_mod.cEq = saved_cEq; + com_mod.dof = saved_dof; + newton_count_ = saved_newton_count; + istr_ = saved_istr; + }; + + try { + const int& cTS = com_mod.cTS; + istr_ = "_" + std::to_string(cTS) + "_residual"; + + 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(); + eq.itr = saved_itr; + + if (com_mod.Rd.size() != 0) { + com_mod.Rd = 0.0; + com_mod.Kd = 0.0; + } + + allocate_linear_system(eq); + set_body_forces(); + assemble_equations(); + if (pre_boundary_conditions) { + pre_boundary_conditions(); + } + try { + apply_boundary_conditions(); + } catch (...) { + if (post_boundary_conditions) { + post_boundary_conditions(); + } + throw; + } + if (post_boundary_conditions) { + post_boundary_conditions(); + } + + if (post_assembly) { + post_assembly(); + } + + if (!eq.assmTLS) { + all_fun::commu(com_mod, com_mod.R); + } + + if (com_mod.sstEq) { + ustruct::ustruct_r(com_mod, solutions_); + } + + if (std::set{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); + } catch (...) { + restore_state(); + throw; + } + + restore_state(); +} + //------------------------ // initiator_step //------------------------ @@ -374,8 +581,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 +828,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 +846,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 +1061,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 +1108,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 +1169,7 @@ void Integrator::corrector() cEq = cEq + 1; } } - #ifdef debug_corrector + #ifdef debug_corrector_cycling dmsg << "eq " << " coupled "; dmsg << "2nd update cEq: " << cEq; #endif diff --git a/Code/Source/solver/Integrator.h b/Code/Source/solver/Integrator.h index 790e17eff..a9a71d3c8 100644 --- a/Code/Source/solver/Integrator.h +++ b/Code/Source/solver/Integrator.h @@ -9,6 +9,38 @@ #include "Vector.h" #include "Simulation.h" +#include + +/// @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,48 @@ 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. + * @param pre_boundary_conditions Optional callback invoked after volume + * assembly and before boundary-condition assembly. + * @param post_boundary_conditions Optional callback invoked immediately after + * boundary-condition assembly. + * @return True if the equation converged, false if max iterations reached + */ + bool step_equation(int iEq, std::function post_assembly = nullptr, + std::function pre_boundary_conditions = nullptr, + std::function post_boundary_conditions = nullptr); + + /** + * @brief Assemble one equation residual at the current solution state. + * + * Reuses the normal assembly and boundary-condition path without solving the + * linear system or applying a Newton correction. The optional callback is + * invoked after boundary-condition assembly and before the standard residual + * communication/row filtering, matching step_equation(). + * + * @param iEq Index of the equation to assemble + * @param post_assembly Optional callback invoked after boundary condition + * application but before residual communication + * @param pre_boundary_conditions Optional callback invoked after volume + * assembly and before boundary-condition assembly. + * @param post_boundary_conditions Optional callback invoked immediately after + * boundary-condition assembly. + */ + void assemble_equation_residual(int iEq, + std::function post_assembly = nullptr, + std::function pre_boundary_conditions = nullptr, + std::function post_boundary_conditions = nullptr); + /** * @brief Perform predictor step for next time step * @@ -173,11 +247,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(); diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index fdc5ae50c..4a4f92495 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -227,14 +227,29 @@ void Parameters::print_parameters() /// @brief Set the simulation parameter values given in an XML format file. void Parameters::read_xml(std::string file_name) -{ +{ tinyxml2::XMLDocument doc; - auto error = doc.LoadFile(file_name.c_str()); - + set_values_from_doc(doc, error, "XML file '" + file_name + "'"); +} + +// Parse parameters from an in-memory XML string (no file on disk). Used by the +// partitioned FSI driver to initialize its sub-simulations directly from the +// single solver.xml without writing temporary per-role XML files. +void Parameters::read_xml_from_string(const std::string& xml_content) +{ + tinyxml2::XMLDocument doc; + auto error = doc.Parse(xml_content.c_str()); + set_values_from_doc(doc, error, "in-memory sub-simulation XML"); +} + +void Parameters::set_values_from_doc(tinyxml2::XMLDocument& doc, + tinyxml2::XMLError error, + const std::string& source_desc) +{ auto root_element = doc.FirstChildElement(FSI_FILE.c_str()); if (error != tinyxml2::XML_SUCCESS || root_element == nullptr) { - svmp::raise("The following error occurred while reading the XML file '" + file_name + "'.\n" + + svmp::raise("The following error occurred while reading the " + source_desc + ".\n" + "[svMultiPhysics] ERROR " + std::string(doc.ErrorStr())); } @@ -253,6 +268,9 @@ void Parameters::read_xml(std::string file_name) // Set mesh projection parameters. set_projection_values(root_element); + // Set partitioned coupling parameters. + set_partitioned_coupling_values(root_element); + // Set Add_equation values. set_equation_values(root_element); @@ -283,6 +301,8 @@ void Parameters::set_equation_values(tinyxml2::XMLElement* root_element) auto eq_params = new EquationParameters(); eq_params->type.set(std::string(eq_type)); + const char* eq_role = add_eq_item->Attribute("role"); + if (eq_role) eq_params->role.set(std::string(eq_role)); eq_params->set_values(add_eq_item); equation_parameters.push_back(eq_params); @@ -3043,6 +3063,50 @@ void ProjectionParameters::set_values(tinyxml2::XMLElement* xml_elem) xml_util_set_parameters(ftpr, xml_elem, error_msg); } +////////////////////////////////////////////////////////// +// PartitionedCouplingParameters // +////////////////////////////////////////////////////////// + +const std::string PartitionedCouplingParameters::xml_element_name_ = "Partitioned_coupling"; + +PartitionedCouplingParameters::PartitionedCouplingParameters() +{ + bool required = true; + + set_parameter("Max_coupling_iterations", 50, !required, max_coupling_iterations); + set_parameter("Coupling_tolerance", 1e-6, !required, coupling_tolerance); + set_parameter("Initial_relaxation", 1.0, !required, initial_relaxation); + set_parameter("Omega_max", 1.0, !required, omega_max); + set_parameter("Coupling_method", "aitken", !required, coupling_method); + set_parameter("Fluid_interface_face", "", required, fluid_interface_face); + set_parameter("Solid_interface_face", "", required, solid_interface_face); + +} + +void PartitionedCouplingParameters::set_values(tinyxml2::XMLElement* xml_elem) +{ + using namespace tinyxml2; + std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; + + using std::placeholders::_1; + using std::placeholders::_2; + + std::function ftpr = + std::bind(&PartitionedCouplingParameters::set_parameter_value, *this, _1, _2); + + xml_util_set_parameters(ftpr, xml_elem, error_msg); + value_set = true; +} + +void Parameters::set_partitioned_coupling_values(tinyxml2::XMLElement* root_element) +{ + auto item = root_element->FirstChildElement(PartitionedCouplingParameters::xml_element_name_.c_str()); + if (item == nullptr) { + return; + } + partitioned_coupling_parameters.set_values(item); +} + ////////////////////////////////////////////////////////// // RISProjectionParameters // ////////////////////////////////////////////////////////// diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 23c377b10..3925d77de 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1637,6 +1637,7 @@ class EquationParameters : public ParameterLists Parameter tolerance; Parameter type; + Parameter role; // "partitioned_fluid", "partitioned_solid", or "partitioned_mesh" Parameter use_taylor_hood_type_basis; // Explicit geometric coupling for FSI simulations: the fluid-structure equations @@ -1923,6 +1924,34 @@ class URISMeshParameters : public ParameterLists +////////////////////////////////////////////////////////// +// PartitionedCouplingParameters // +////////////////////////////////////////////////////////// + +/// @brief Parameters for the 'Partitioned_coupling' XML element. +/// +/// Configures partitioned FSI coupling between separately solved +/// fluid and solid equations with Aitken relaxation. +class PartitionedCouplingParameters : public ParameterLists +{ + public: + PartitionedCouplingParameters(); + static const std::string xml_element_name_; + void set_values(tinyxml2::XMLElement* xml_elem); + bool defined() const { return value_set; } + + Parameter max_coupling_iterations; + Parameter coupling_tolerance; + Parameter initial_relaxation; + Parameter omega_max; + Parameter coupling_method; // "constant" or "aitken" + Parameter fluid_interface_face; + Parameter solid_interface_face; + + bool value_set = false; +}; + + /// @brief The Parameters class stores parameter values read in from a solver input file. class Parameters { @@ -1936,6 +1965,9 @@ class Parameters { void get_logging_levels(int& verbose, int& warning, int& debug); void print_parameters(); void read_xml(std::string file_name); + void read_xml_from_string(const std::string& xml_content); + void set_values_from_doc(tinyxml2::XMLDocument& doc, tinyxml2::XMLError error, + const std::string& source_desc); void set_contact_values(tinyxml2::XMLElement* root_element); void set_equation_values(tinyxml2::XMLElement* root_element); @@ -1946,6 +1978,7 @@ class Parameters { void set_RIS_projection_values(tinyxml2::XMLElement* root_element); void set_URIS_mesh_values(tinyxml2::XMLElement* root_element); + void set_partitioned_coupling_values(tinyxml2::XMLElement* root_element); // Objects representing each parameter section of XML file. ContactParameters contact_parameters; @@ -1958,6 +1991,8 @@ class Parameters { std::vector RIS_projection_parameters; std::vector URIS_mesh_parameters; + PartitionedCouplingParameters partitioned_coupling_parameters; + }; #endif diff --git a/Code/Source/solver/PartitionedFSI.cpp b/Code/Source/solver/PartitionedFSI.cpp new file mode 100644 index 000000000..686031971 --- /dev/null +++ b/Code/Source/solver/PartitionedFSI.cpp @@ -0,0 +1,1257 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +#include "PartitionedFSI.h" +#include "Integrator.h" +#include "all_fun.h" +#include "consts.h" +#include "fsi_coupling.h" +#include "set_bc.h" +#include "distribute.h" +#include "initialize.h" +#include "output.h" +#include "vtk_xml.h" +#include "read_files.h" +#include "utils.h" + +#include +#include +#include +#include +#include +#include + +/// Check if any value in the solution arrays is NaN +static bool has_nan(const SolutionStates& sol) { + const Array* arrays[] = { + &sol.current.get_velocity(), + &sol.current.get_acceleration(), + &sol.current.get_displacement() + }; + for (auto* arr : arrays) + for (int a = 0; a < arr->ncols(); a++) + for (int i = 0; i < arr->nrows(); i++) + if (std::isnan((*arr)(i, a))) return true; + return false; +} + +//---------------------------------------------------------------------- +// Helper: initialize one sub-simulation through the standard pipeline. +// xml_content is the in-memory sub-XML produced by build_sub_xml (no file +// is written to disk). +//---------------------------------------------------------------------- +static void init_sub_sim(Simulation* sim, const std::string& xml_content) +{ + read_files_ns::read_files(sim, xml_content, /*from_string=*/true); + sim->logger.set_cout_write(false); + + // The mesh sub-sim includes a stub so that read_files + // accepts the mesh equation, but this sim has only the mesh equation (tDof=3). + // baf_ini assumes FSI DOF layout (Do(i+nsd+1,Ac) reaches row 4) when mvMsh=true, + // which would go out-of-bounds for the standalone mesh sub-sim. + if (sim->com_mod.nEq == 1 && + sim->com_mod.eq[0].phys == consts::EquationType::phys_mesh) { + sim->com_mod.mvMsh = false; + } + + distribute(sim); + Vector init_time(3); + initialize(sim, init_time); + for (int iEq = 0; iEq < sim->com_mod.nEq; iEq++) { + add_eq_linear_algebra(sim->com_mod, sim->com_mod.eq[iEq]); + } +} + +//---------------------------------------------------------------------- +// build_sub_xml — extract one role's equation + its meshes from the +// main XML and return a minimal standalone sub-simulation XML as an +// in-memory string (consumed by read_files with from_string=true). +// +// Mesh association: the meshes included in the sub-XML are those whose +// value matches any child of the target +// equation. The partitioned_mesh role falls back to the fluid meshes +// when the mesh equation carries no explicit domain block. +//---------------------------------------------------------------------- +std::string PartitionedFSI::build_sub_xml(const std::string& main_xml_path, + const std::string& role) +{ + using namespace tinyxml2; + + XMLDocument doc; + if (doc.LoadFile(main_xml_path.c_str()) != XML_SUCCESS) + throw std::runtime_error("[PartitionedFSI] Cannot parse main XML: " + main_xml_path); + + XMLElement* root = doc.FirstChildElement("svMultiPhysicsFile"); + if (!root) + throw std::runtime_error("[PartitionedFSI] Missing root in " + main_xml_path); + + // Find the equation element with the requested role attribute. + XMLElement* target_eq = nullptr; + for (XMLElement* eq = root->FirstChildElement("Add_equation"); eq; + eq = eq->NextSiblingElement("Add_equation")) { + const char* r = eq->Attribute("role"); + if (r && std::string(r) == role) { target_eq = eq; break; } + } + if (!target_eq) + throw std::runtime_error("[PartitionedFSI] No found in " + main_xml_path); + + // Collect domain IDs used by the target equation. + std::set domain_ids; + for (XMLElement* d = target_eq->FirstChildElement("Domain"); d; + d = d->NextSiblingElement("Domain")) + domain_ids.insert(d->IntAttribute("id", -1)); + + // Mesh role with no domain block: use the same domains as the fluid equation. + if (domain_ids.empty() && role == "partitioned_mesh") { + for (XMLElement* eq = root->FirstChildElement("Add_equation"); eq; + eq = eq->NextSiblingElement("Add_equation")) { + const char* r = eq->Attribute("role"); + if (r && std::string(r) == "partitioned_fluid") { + for (XMLElement* d = eq->FirstChildElement("Domain"); d; + d = d->NextSiblingElement("Domain")) + domain_ids.insert(d->IntAttribute("id", -1)); + break; + } + } + } + + // Build sub-document. + XMLDocument sub; + XMLElement* sub_root = sub.NewElement("svMultiPhysicsFile"); + sub_root->SetAttribute("version", "0.1"); + sub.InsertFirstChild(sub_root); + + // Copy GeneralSimulationParameters, overriding the VTK save prefix + // so each sub-sim writes to distinct files (result_fluid_*, result_solid_*, result_mesh_*). + XMLElement* gen = root->FirstChildElement("GeneralSimulationParameters"); + if (gen) { + XMLElement* gen_clone = gen->DeepClone(&sub)->ToElement(); + // Map role → short suffix used in the prefix name + std::string suffix; + if (role == "partitioned_fluid") suffix = "fluid"; + else if (role == "partitioned_solid") suffix = "solid"; + else if (role == "partitioned_mesh") suffix = "mesh"; + if (!suffix.empty()) { + XMLElement* name_elem = gen_clone->FirstChildElement("Name_prefix_of_saved_VTK_files"); + if (name_elem) { + std::string base_prefix = name_elem->GetText() ? name_elem->GetText() : "result"; + // trim whitespace + base_prefix.erase(0, base_prefix.find_first_not_of(" \t\n\r")); + base_prefix.erase(base_prefix.find_last_not_of(" \t\n\r") + 1); + name_elem->SetText((" " + base_prefix + "_" + suffix + " ").c_str()); + } + } + sub_root->InsertEndChild(gen_clone); + } + + // Copy matching Add_mesh elements. + for (XMLElement* mesh = root->FirstChildElement("Add_mesh"); mesh; + mesh = mesh->NextSiblingElement("Add_mesh")) { + XMLElement* dom_elem = mesh->FirstChildElement("Domain"); + if (!dom_elem) continue; + int mesh_dom = dom_elem->IntText(-1); + if (domain_ids.empty() || domain_ids.count(mesh_dom)) + sub_root->InsertEndChild(mesh->DeepClone(&sub)); + } + + // Clone the equation, stripping the role attribute. + XMLElement* eq_clone = target_eq->DeepClone(&sub)->ToElement(); + eq_clone->DeleteAttribute("role"); + sub_root->InsertEndChild(eq_clone); + + // The mesh sub-sim needs a minimal block so that + // read_files sets mvMsh=true (required for the mesh equation to be valid). + if (role == "partitioned_mesh") { + XMLElement* pcp = sub.NewElement("Partitioned_coupling"); + XMLElement* fface = sub.NewElement("Fluid_interface_face"); + fface->SetText("dummy"); + XMLElement* sface = sub.NewElement("Solid_interface_face"); + sface->SetText("dummy"); + pcp->InsertEndChild(fface); + pcp->InsertEndChild(sface); + sub_root->InsertEndChild(pcp); + } + + // Serialize the sub-document to an in-memory string instead of writing a + // temp file. The string is handed directly to read_files(..., from_string). + tinyxml2::XMLPrinter printer; + sub.Print(&printer); + return std::string(printer.CStr()); +} + +//---------------------------------------------------------------------- +// Constructor +//---------------------------------------------------------------------- +PartitionedFSI::PartitionedFSI(Simulation* main_simulation, + const PartitionedFSIConfig& config, + const std::string& xml_file_path) + : main_sim_(main_simulation), config_(config), + xml_file_path_(xml_file_path), omega_(config.initial_relaxation) +{ + auto& cm = main_sim_->com_mod.cm; + auto& cm_mod = main_sim_->cm_mod; + + // Build each role's sub-XML in memory from the single main solver.xml and + // initialize a standalone sub-simulation from it. No temp files are written: + // every rank parses the same main XML (read-only) and passes the resulting + // sub-XML string straight to read_files, so no rank-0 write / barrier is needed. + fluid_sim_ = std::make_unique(); + init_sub_sim(fluid_sim_.get(), build_sub_xml(xml_file_path_, "partitioned_fluid")); + + solid_sim_ = std::make_unique(); + init_sub_sim(solid_sim_.get(), build_sub_xml(xml_file_path_, "partitioned_solid")); + // Monolithic FSI advances the solid domain with the FSI equation's + // generalized-alpha coefficients, not standalone structural coefficients. + fsi_coupling::copy_time_integration_parameters( + fluid_sim_->com_mod.eq[0], solid_sim_->com_mod.eq[0]); + + mesh_sim_ = std::make_unique(); + init_sub_sim(mesh_sim_.get(), build_sub_xml(xml_file_path_, "partitioned_mesh")); + + if (cm.mas(cm_mod)) { + // Open log files + std::string log_dir = fluid_sim_->get_chnl_mod().appPath; + coupling_log_.open(log_dir + "coupling.dat"); + char hdr[256]; + snprintf(hdr, sizeof(hdr), "# %4s %3s %10s %5s %10s %10s %10s %10s", + "cTS", "cp", "time", "dB", "Ri/R1", "Ri/R0", "omega", "|disp|"); + coupling_log_ << hdr << std::endl; + + histor_log_.open(log_dir + "histor.dat"); + } + + resolve_faces(); + build_node_maps(); + build_solid_interface_force_mask(); +} + +PartitionedFSI::~PartitionedFSI() = default; + +//---------------------------------------------------------------------- +// resolve_faces +//---------------------------------------------------------------------- +void PartitionedFSI::resolve_faces() +{ + auto find_face = [](Simulation* sim, const std::string& face_name, + const faceType*& face_out, const mshType*& mesh_out) { + for (int iM = 0; iM < sim->com_mod.nMsh; iM++) { + auto& msh = sim->com_mod.msh[iM]; + for (int iFa = 0; iFa < msh.nFa; iFa++) { + if (msh.fa[iFa].name == face_name) { + face_out = &msh.fa[iFa]; + mesh_out = &msh; + return; + } + } + } + throw std::runtime_error("[PartitionedFSI] Face '" + face_name + "' not found."); + }; + + find_face(fluid_sim_.get(), config_.fluid_interface_face, fluid_face_, fluid_mesh_); + find_face(solid_sim_.get(), config_.solid_interface_face, solid_face_, solid_mesh_); + find_face(mesh_sim_.get(), config_.fluid_interface_face, mesh_face_, mesh_mesh_); +} + +//---------------------------------------------------------------------- +// compute_face_global_info — gather per-rank nNo to compute global nNo +// and this rank's offset within the global face node ordering. +//---------------------------------------------------------------------- +void PartitionedFSI::compute_face_global_info( + const faceType& face, cmType& cm, const CmMod& cm_mod, + int& global_nNo, int& local_offset) +{ + int np = cm.np(); + int my_rank = cm.id(); + int local_nNo = face.nNo; + + std::vector all_nNo(np); + MPI_Allgather(&local_nNo, 1, MPI_INT, + all_nNo.data(), 1, MPI_INT, cm.com()); + + global_nNo = 0; + local_offset = 0; + for (int p = 0; p < np; p++) { + if (p < my_rank) local_offset += all_nNo[p]; + global_nNo += all_nNo[p]; + } +} + +//---------------------------------------------------------------------- +// gather_face_data — MPI_Allgatherv local face data to all ranks. +// local_data is (nrows, local_nNo); returns (nrows, global_nNo). +// Rank ordering in global array matches the local_offset convention. +//---------------------------------------------------------------------- +Array PartitionedFSI::gather_face_data( + const Array& local_data, + int global_nNo, int /*local_offset*/, + cmType& cm, const CmMod& cm_mod) +{ + const int nrows = local_data.nrows(); + const int local_nNo = local_data.ncols(); + const int np = cm.np(); + + // Pack as flat row-major: [node0_row0, node0_row1, ..., node1_row0, ...] + std::vector local_flat(nrows * local_nNo); + for (int a = 0; a < local_nNo; a++) + for (int i = 0; i < nrows; i++) + local_flat[a * nrows + i] = local_data(i, a); + + int send_count = nrows * local_nNo; + std::vector recv_counts(np), displs(np); + MPI_Allgather(&send_count, 1, MPI_INT, + recv_counts.data(), 1, MPI_INT, cm.com()); + int total = 0; + for (int p = 0; p < np; p++) { displs[p] = total; total += recv_counts[p]; } + + std::vector global_flat(total); + MPI_Allgatherv(local_flat.data(), send_count, MPI_DOUBLE, + global_flat.data(), recv_counts.data(), displs.data(), + MPI_DOUBLE, cm.com()); + + Array result(nrows, global_nNo); + int node_offset = 0; + for (int p = 0; p < np; p++) { + int p_nNo = recv_counts[p] / nrows; + for (int la = 0; la < p_nNo; la++) + for (int i = 0; i < nrows; i++) + result(i, node_offset + la) = global_flat[displs[p] + la * nrows + i]; + node_offset += p_nNo; + } + return result; +} + +//---------------------------------------------------------------------- +// gather_global_map — all-gather a local (local_src → global_tgt) map +// into a global (global_src → global_tgt) map. +//---------------------------------------------------------------------- +void PartitionedFSI::gather_global_map( + const std::vector& local_map, + int global_src_nNo, + cmType& cm, const CmMod& cm_mod, + std::vector& global_map) +{ + int np = cm.np(); + int send_count = static_cast(local_map.size()); + + std::vector recv_counts(np), displs(np); + MPI_Allgather(&send_count, 1, MPI_INT, + recv_counts.data(), 1, MPI_INT, cm.com()); + int total = 0; + for (int p = 0; p < np; p++) { displs[p] = total; total += recv_counts[p]; } + + std::vector global_flat(total); + MPI_Allgatherv(local_map.data(), send_count, MPI_INT, + global_flat.data(), recv_counts.data(), displs.data(), + MPI_INT, cm.com()); + + global_map.assign(global_src_nNo, -1); + int offset = 0; + for (int p = 0; p < np; p++) { + for (int la = 0; la < recv_counts[p]; la++) + global_map[offset + la] = global_flat[displs[p] + la]; + offset += recv_counts[p]; + } +} + +//---------------------------------------------------------------------- +// build_face_node_map — match each LOCAL face_a node to its nearest +// node in the PRE-GATHERED global face_b coordinates. +// Returns local_src_idx → global_tgt_idx map. +//---------------------------------------------------------------------- +void PartitionedFSI::build_face_node_map( + const faceType& face_a, const ComMod& com_a, + int global_b_nNo, const Array& global_b_coords, + std::vector& a_to_global_b) +{ + const int nsd = com_a.nsd; + const double tol = 1e-8; + a_to_global_b.assign(face_a.nNo, -1); + + for (int a = 0; a < face_a.nNo; a++) { + int Ac = face_a.gN(a); + double best = 1e30; + int best_b = -1; + for (int bg = 0; bg < global_b_nNo; bg++) { + double d2 = 0.0; + for (int i = 0; i < nsd; i++) { + double d = com_a.x(i, Ac) - global_b_coords(i, bg); + d2 += d * d; + } + if (d2 < best) { best = d2; best_b = bg; } + } + if (best < tol * tol) a_to_global_b[a] = best_b; + } +} + +//---------------------------------------------------------------------- +// build_node_maps — build global→global interface node maps. +// Each sub-mesh is independently distributed, so we gather all face +// coordinates from all ranks before performing the nearest-neighbor +// search; then gather the local partial maps into global maps. +//---------------------------------------------------------------------- +void PartitionedFSI::build_node_maps() +{ + auto& cm = main_sim_->com_mod.cm; + auto& cm_mod = main_sim_->cm_mod; + const int nsd = main_sim_->com_mod.nsd; + + // Compute global nNo and per-rank offsets for each interface face + compute_face_global_info(*solid_face_, cm, cm_mod, + solid_face_global_nNo_, solid_face_local_offset_); + compute_face_global_info(*fluid_face_, cm, cm_mod, + fluid_face_global_nNo_, fluid_face_local_offset_); + compute_face_global_info(*mesh_face_, cm, cm_mod, + mesh_face_global_nNo_, mesh_face_local_offset_); + + // Gather face coordinates globally for each sub-mesh face + auto pack_coords = [&](const faceType& face, const ComMod& com) { + Array local(nsd, face.nNo); + for (int a = 0; a < face.nNo; a++) { + int Ac = face.gN(a); + for (int i = 0; i < nsd; i++) + local(i, a) = com.x(i, Ac); + } + return local; + }; + + auto global_solid_coords = gather_face_data( + pack_coords(*solid_face_, solid_sim_->com_mod), + solid_face_global_nNo_, solid_face_local_offset_, cm, cm_mod); + auto global_fluid_coords = gather_face_data( + pack_coords(*fluid_face_, fluid_sim_->com_mod), + fluid_face_global_nNo_, fluid_face_local_offset_, cm, cm_mod); + auto global_mesh_coords = gather_face_data( + pack_coords(*mesh_face_, mesh_sim_->com_mod), + mesh_face_global_nNo_, mesh_face_local_offset_, cm, cm_mod); + + // Build local (local_src → global_tgt) maps, then gather to global maps + std::vector local_s2f, local_f2s, local_s2m; + build_face_node_map(*solid_face_, solid_sim_->com_mod, + fluid_face_global_nNo_, global_fluid_coords, local_s2f); + build_face_node_map(*fluid_face_, fluid_sim_->com_mod, + solid_face_global_nNo_, global_solid_coords, local_f2s); + build_face_node_map(*solid_face_, solid_sim_->com_mod, + mesh_face_global_nNo_, global_mesh_coords, local_s2m); + + gather_global_map(local_s2f, solid_face_global_nNo_, cm, cm_mod, solid_to_fluid_map_); + gather_global_map(local_f2s, fluid_face_global_nNo_, cm, cm_mod, fluid_to_solid_map_); + gather_global_map(local_s2m, solid_face_global_nNo_, cm, cm_mod, solid_to_mesh_map_); + + solid_face_canonical_ = build_face_dedup_map(*solid_face_, solid_sim_->com_mod, + solid_face_global_nNo_, cm, cm_mod); + fluid_face_canonical_ = build_face_dedup_map(*fluid_face_, fluid_sim_->com_mod, + fluid_face_global_nNo_, cm, cm_mod); + mesh_face_canonical_ = build_face_dedup_map(*mesh_face_, mesh_sim_->com_mod, + mesh_face_global_nNo_, cm, cm_mod); + + // Owner flags for the solid interface face: the recovered interface force is + // a precomputed nodal load that must be injected exactly once per physical + // node, since all_fun::commu sums it across ranks after the callback. + solid_face_owner_ = build_face_owner_flags(*solid_face_, solid_sim_->com_mod, cm); + + if (cm.mas(cm_mod)) { + auto count_dups = [](const std::vector& c) { + int n = 0; for (int i = 0; i < (int)c.size(); i++) if (c[i] != i) n++; return n; + }; + std::cout << "[PartitionedFSI] face dedup: solid " + << count_dups(solid_face_canonical_) << "/" << solid_face_global_nNo_ + << " dups, fluid " << count_dups(fluid_face_canonical_) << "/" << fluid_face_global_nNo_ + << " dups, mesh " << count_dups(mesh_face_canonical_) << "/" << mesh_face_global_nNo_ + << " dups" << std::endl; + } +} + +//---------------------------------------------------------------------- +// build_solid_interface_force_mask — mask interface force components that +// are prescribed by strong solid Dirichlet boundary conditions. +//---------------------------------------------------------------------- +void PartitionedFSI::build_solid_interface_force_mask() +{ + auto& solid_com = solid_sim_->com_mod; + const int nsd = main_sim_->com_mod.nsd; + + solid_interface_force_mask_.resize(nsd, solid_face_->nNo); + solid_interface_force_mask_ = 1.0; + + std::unordered_map interface_node_to_face_index; + interface_node_to_face_index.reserve(solid_face_->nNo); + for (int a = 0; a < solid_face_->nNo; a++) { + interface_node_to_face_index.emplace(solid_face_->gN(a), a); + } + + if (solid_com.nEq == 0) { + return; + } + + const auto& solid_eq = solid_com.eq[0]; + for (int iBc = 0; iBc < solid_eq.nBc; iBc++) { + const auto& bc = solid_eq.bc[iBc]; + if (!utils::btest(bc.bType, consts::iBC_Dir) || bc.weakDir) { + continue; + } + if (bc.iM < 0 || bc.iM >= solid_com.nMsh) { + continue; + } + const auto& mesh = solid_com.msh[bc.iM]; + if (bc.iFa < 0 || bc.iFa >= mesh.nFa) { + continue; + } + + bool has_effective_direction = false; + for (int i = 0; i < nsd; i++) { + if (i < bc.eDrn.size() && bc.eDrn(i) != 0) { + has_effective_direction = true; + break; + } + } + + const auto& bc_face = mesh.fa[bc.iFa]; + for (int a = 0; a < bc_face.nNo; a++) { + const auto it = interface_node_to_face_index.find(bc_face.gN(a)); + if (it == interface_node_to_face_index.end()) { + continue; + } + + const int interface_a = it->second; + for (int i = 0; i < nsd; i++) { + if (!has_effective_direction || (i < bc.eDrn.size() && bc.eDrn(i) != 0)) { + solid_interface_force_mask_(i, interface_a) = 0.0; + } + } + } + } +} + +//---------------------------------------------------------------------- +// build_face_dedup_map — canonical[i] = first global slot index that +// shares the same physical node (global node ID) as slot i. +// Slots with no duplicate have canonical[i] == i. +//---------------------------------------------------------------------- +std::vector PartitionedFSI::build_face_dedup_map( + const faceType& face, const ComMod& com, + int global_nNo, cmType& cm, const CmMod& cm_mod) +{ + std::vector local_ids(face.nNo); + for (int a = 0; a < face.nNo; a++) + local_ids[a] = com.ltg[face.gN(a)]; + + int np = cm.np(); + std::vector counts(np), displs(np); + MPI_Allgather(&face.nNo, 1, MPI_INT, counts.data(), 1, MPI_INT, cm.com()); + int total = 0; + for (int p = 0; p < np; p++) { displs[p] = total; total += counts[p]; } + + std::vector all_ids(total); + MPI_Allgatherv(local_ids.data(), face.nNo, MPI_INT, + all_ids.data(), counts.data(), displs.data(), MPI_INT, cm.com()); + + std::unordered_map first; + std::vector canonical(total); + for (int i = 0; i < total; i++) { + auto [it, inserted] = first.emplace(all_ids[i], i); + canonical[i] = it->second; + } + return canonical; +} + +//---------------------------------------------------------------------- +// apply_dedup — copy canonical slot value to every duplicate slot. +//---------------------------------------------------------------------- +void PartitionedFSI::apply_dedup(Array& arr, const std::vector& canonical) +{ + for (int i = 0; i < (int)canonical.size(); i++) + if (canonical[i] != i) + for (int row = 0; row < arr.nrows(); row++) + arr(row, i) = arr(row, canonical[i]); +} + +//---------------------------------------------------------------------- +// build_face_owner_flags — flags[a] = 1 if THIS rank is the lowest-ID rank +// that contains local face node a, else 0. A node shared across ranks +// (a partition-ring node) is owned by exactly one rank under this rule. +//---------------------------------------------------------------------- +std::vector PartitionedFSI::build_face_owner_flags( + const faceType& face, const ComMod& com, cmType& cm) +{ + const int np = cm.np(); + const int myrank = cm.id(); + + std::vector local_ids(face.nNo); + for (int a = 0; a < face.nNo; a++) + local_ids[a] = com.ltg[face.gN(a)]; + + std::vector counts(np), displs(np); + MPI_Allgather(&face.nNo, 1, MPI_INT, counts.data(), 1, MPI_INT, cm.com()); + int total = 0; + for (int p = 0; p < np; p++) { displs[p] = total; total += counts[p]; } + + std::vector all_ids(total); + MPI_Allgatherv(local_ids.data(), face.nNo, MPI_INT, + all_ids.data(), counts.data(), displs.data(), MPI_INT, cm.com()); + + // owner[globalID] = lowest rank containing it (ranks scanned in ascending order) + std::unordered_map owner; + for (int p = 0; p < np; p++) + for (int k = 0; k < counts[p]; k++) + owner.emplace(all_ids[displs[p] + k], p); + + std::vector flags(face.nNo); + for (int a = 0; a < face.nNo; a++) + flags[a] = (owner[local_ids[a]] == myrank) ? 1 : 0; + return flags; +} + +//---------------------------------------------------------------------- +// transfer_data +//---------------------------------------------------------------------- +Array PartitionedFSI::transfer_data( + const std::vector& src_to_tgt_map, + const Array& src_data, int tgt_nNo) +{ + int nrows = src_data.nrows(); + Array result(nrows, tgt_nNo); + for (int a = 0; a < static_cast(src_to_tgt_map.size()); a++) { + int b = src_to_tgt_map[a]; + if (b >= 0) { + for (int i = 0; i < nrows; i++) result(i, b) = src_data(i, a); + } + } + return result; +} + +//---------------------------------------------------------------------- +// relax_interface — updates disp_prev_ and vel_prev_ +//---------------------------------------------------------------------- +void PartitionedFSI::relax_interface(int cp, int nsd, + const Array& disp_current) +{ + switch (config_.coupling_method) { + case CouplingMethod::constant: + relax_constant(cp, nsd, disp_current); + break; + case CouplingMethod::aitken: + relax_aitken(cp, nsd, disp_current); + break; + } +} + +//---------------------------------------------------------------------- +// relax_constant — fixed relaxation (operates on global face arrays) +//---------------------------------------------------------------------- +void PartitionedFSI::relax_constant(int cp, int nsd, + const Array& disp_current) +{ + omega_ = config_.initial_relaxation; + for (int a = 0; a < solid_face_global_nNo_; a++) + for (int i = 0; i < nsd; i++) + disp_prev_(i, a) += omega_ * (disp_current(i, a) - disp_prev_(i, a)); +} + +//---------------------------------------------------------------------- +// relax_aitken — Aitken Delta^2 (Küttler & Wall 2008, Eq. 44) +// Operates on global face arrays; all ranks have identical data so no +// MPI reduction is needed here. +//---------------------------------------------------------------------- +void PartitionedFSI::relax_aitken(int cp, int nsd, + const Array& disp_current) +{ + const int u = nsd * solid_face_global_nNo_; + + // Build residual r = x_tilde - x + std::vector r(u); + for (int a = 0; a < solid_face_global_nNo_; a++) + for (int i = 0; i < nsd; i++) + r[a * nsd + i] = disp_current(i, a) - disp_prev_(i, a); + + // Aitken update: omega = -omega * r^T (r_new - r_old) / |r_new - r_old|^2 + // Negative omega allowed (corrects overshoot) + if (cp > 0 && !r_prev_.empty()) { + double num = 0, den = 0; + for (int j = 0; j < u; j++) { + // Count each physical node once so omega is independent of partition + // count (duplicate ring slots hold identical residual entries). + if (solid_face_canonical_[j / nsd] != j / nsd) continue; + double dr = r[j] - r_prev_[j]; + num += r_prev_[j] * dr; + den += dr * dr; + } + if (den > 1e-30) { + omega_ = -omega_ * num / den; + if (std::abs(omega_) > config_.omega_max) + omega_ = (omega_ > 0) ? config_.omega_max : -config_.omega_max; + } + } + r_prev_ = r; + + // Apply: x_{k+1} = x_k + omega * r + for (int a = 0; a < solid_face_global_nNo_; a++) + for (int i = 0; i < nsd; i++) + disp_prev_(i, a) += omega_ * (disp_current(i, a) - disp_prev_(i, a)); +} + +//====================================================================== +// run — full time-stepping loop with Dirichlet-Neumann coupling +//====================================================================== +void PartitionedFSI::run() +{ + auto& main_com = main_sim_->com_mod; + auto& cm_mod = main_sim_->cm_mod; + auto& cm = main_com.cm; + + int nTS = main_com.nTS; + int& cTS = main_com.cTS; + double& dt = main_com.dt; + double& time = main_com.time; + int nITs = main_com.nITs; + + if (cTS <= nITs) dt = dt / 10.0; + + Simulation* sims[3] = {fluid_sim_.get(), solid_sim_.get(), mesh_sim_.get()}; + + while (true) { + if (cTS == nITs) dt = 10.0 * dt; + cTS = cTS + 1; + time = time + dt; + + // Sync time to sub-sims + for (auto* sim : sims) { + sim->com_mod.cTS = cTS; + sim->com_mod.time = time; + sim->com_mod.dt = dt; + for (auto& eq : sim->com_mod.eq) { eq.itr = 0; eq.ok = false; } + } + + if (cm.mas(cm_mod)) { + if (histor_log_.is_open()) { + histor_log_ << std::string(70, '=') << std::endl; + histor_log_ << " TIME STEP " << cTS << " t=" << time << " dt=" << dt << std::endl; + histor_log_ << std::string(70, '=') << std::endl; + } + } + + // Predictor + Dirichlet BCs for each sub-sim + for (auto* sim : sims) { + sim->get_integrator().predictor(); + set_bc::set_bc_dir(sim->com_mod, sim->get_integrator().get_solutions()); + } + + // Coupling loop + bool converged = step(); + + if (!converged && cm.mas(cm_mod)) { + std::cout << " TIME STEP " << cTS << " FAILED (NaN or no convergence)" << std::endl; + if (histor_log_.is_open()) + histor_log_ << " TIME STEP " << cTS << " FAILED (NaN or no convergence)" << std::endl; + } + + // Stop on failure + if (!converged) break; + + // Save results + save_results(); + + // Copy current -> old + for (auto* sim : sims) { + auto& sol = sim->get_integrator().get_solutions(); + sol.old.get_acceleration() = sol.current.get_acceleration(); + sol.old.get_velocity() = sol.current.get_velocity(); + if (sim->com_mod.dFlag) + sol.old.get_displacement() = sol.current.get_displacement(); + } + + // Stop condition + int stopTS = nTS; + if (cm.mas(cm_mod)) { + if (FILE* fp = fopen(main_com.stopTrigName.c_str(), "r")) { + int count = fscanf(fp, "%d", &stopTS); + if (count == 0) stopTS = cTS; + fclose(fp); + } + } + cm.bcast(cm_mod, &stopTS); + if (cTS >= stopTS) break; + } +} + +//---------------------------------------------------------------------- +// compute_interface_velocity — Newmark-consistent velocity from disp_prev_. +// Each rank computes its local solid face nodes, then all-gather to +// produce the global vel_prev_ replicated on all ranks. +//---------------------------------------------------------------------- +void PartitionedFSI::compute_interface_velocity() +{ + auto& solid_com = solid_sim_->com_mod; + auto& solid_sol = solid_sim_->get_integrator().get_solutions(); + const auto& eq = solid_com.eq[0]; + const int s = eq.s; + const int nsd = main_sim_->com_mod.nsd; + const double dt = solid_com.dt; + const auto& Do = solid_sol.old.get_displacement(); + const auto& Yo = solid_sol.old.get_velocity(); + const auto& Ao = solid_sol.old.get_acceleration(); + auto& cm = main_sim_->com_mod.cm; + auto& cm_mod = main_sim_->cm_mod; + + // Compute velocity for this rank's local solid face nodes + Array local_vel(nsd, solid_face_->nNo); + for (int a = 0; a < solid_face_->nNo; a++) { + int Ac = solid_face_->gN(a); + for (int i = 0; i < nsd; i++) { + double disp_a = disp_prev_(i, solid_face_local_offset_ + a); + double a_new, v_new; + newmark::state_from_displacement( + disp_a, Do(i + s, Ac), Yo(i + s, Ac), Ao(i + s, Ac), + dt, eq.beta, eq.gam, a_new, v_new); + local_vel(i, a) = v_new; + } + } + + // All-gather to global. Duplicate ring slots are filled identically by + // their owning ranks (velocity is computed from the replicated disp_prev_), + // so no dedup is needed. + vel_prev_ = gather_face_data(local_vel, solid_face_global_nNo_, + solid_face_local_offset_, cm, cm_mod); +} + +//---------------------------------------------------------------------- +// solve_fluid — fluid equation with interface velocity and ALE. +// vel_prev_ is global; extract this rank's local fluid face portion. +//---------------------------------------------------------------------- +bool PartitionedFSI::solve_fluid( + const Array& fluid_x_old, + const Array& mesh_vel_Yo, const Array& mesh_vel_Yn) +{ + auto& fluid_com = fluid_sim_->com_mod; + auto& fluid_int = fluid_sim_->get_integrator(); + auto& fluid_sol = fluid_int.get_solutions(); + const int nsd = main_sim_->com_mod.nsd; + + // Transfer global solid velocity → global fluid velocity, then extract local. + // transfer_data only writes the canonical target slot of each node; dedup + // fills the duplicate ring slots so every rank reads a correct Dirichlet + // value below (the velocity is applied on all face nodes, not just owners). + auto global_fluid_vel = transfer_data(solid_to_fluid_map_, vel_prev_, + fluid_face_global_nNo_); + apply_dedup(global_fluid_vel, fluid_face_canonical_); + Array local_fluid_vel(nsd, fluid_face_->nNo); + for (int a = 0; a < fluid_face_->nNo; a++) + for (int i = 0; i < nsd; i++) + local_fluid_vel(i, a) = global_fluid_vel(i, fluid_face_local_offset_ + a); + + set_bc::set_bc_dir(fluid_com, fluid_sol); + fsi_coupling::apply_velocity_on_fluid( + fluid_com, fluid_com.eq[0], *fluid_face_, local_fluid_vel, fluid_sol); + + // ALE mesh velocity at generalized-alpha intermediate time + double af = fluid_com.eq[0].af; + fluid_com.ale_mesh_velocity.resize(nsd, fluid_com.tnNo); + for (int a = 0; a < fluid_com.tnNo; a++) + for (int i = 0; i < nsd; i++) + fluid_com.ale_mesh_velocity(i, a) = (1.0 - af) * mesh_vel_Yo(i, a) + + af * mesh_vel_Yn(i, a); + + // On the FSI interface, monolithic residual assembly has ALE velocity at + // the same generalized-alpha stage as the imposed fluid velocity. + auto sync_interface_ale_to_fluid_stage = [&]() { + const auto& Yo = fluid_sol.old.get_velocity(); + const auto& Yn = fluid_sol.current.get_velocity(); + const int s = fluid_com.eq[0].s; + for (int a = 0; a < fluid_face_->nNo; a++) { + const int Ac = fluid_face_->gN(a); + for (int i = 0; i < nsd; i++) { + fluid_com.ale_mesh_velocity(i, Ac) = + (1.0 - af) * Yo(s + i, Ac) + af * Yn(s + i, Ac); + } + } + }; + sync_interface_ale_to_fluid_stage(); + + // Monolithic FSI assembles the fluid volume on the generalized-alpha ALE + // geometry, but its face-based Neumann BC path uses the old mesh geometry + // through gnnb(). Match that split while boundary conditions are assembled. + Array saved_fluid_x_for_bc; + auto use_old_geometry_for_boundary_conditions = [&]() { + if (fluid_x_old.nrows() != fluid_com.x.nrows() || + fluid_x_old.ncols() != fluid_com.x.ncols()) { + throw std::runtime_error( + "[PartitionedFSI] Old fluid coordinates are incompatible with fluid coordinates."); + } + saved_fluid_x_for_bc = fluid_com.x; + fluid_com.x = fluid_x_old; + }; + auto restore_staged_geometry_after_boundary_conditions = [&]() { + fluid_com.x = saved_fluid_x_for_bc; + }; + + bool fluid_ok = fluid_int.step_equation(0, [&]() { + set_bc::enforce_dirichlet_dofs_on_face(fluid_com, *fluid_face_, 0, nsd); + }, use_old_geometry_for_boundary_conditions, + restore_staged_geometry_after_boundary_conditions); + + if (fluid_ok && !has_nan(fluid_sol)) { + set_bc::set_bc_dir(fluid_com, fluid_sol); + fsi_coupling::apply_velocity_on_fluid( + fluid_com, fluid_com.eq[0], *fluid_face_, local_fluid_vel, fluid_sol); + sync_interface_ale_to_fluid_stage(); + + auto recover_interface_force_before_boundary_conditions = [&]() { + Array fluid_residual = fluid_com.R; + if (!fluid_com.eq[fluid_com.cEq].assmTLS) { + all_fun::commu(fluid_com, fluid_residual); + } + local_fluid_interface_force_ = + fsi_coupling::extract_fluid_residual_force( + fluid_com, *fluid_face_, fluid_residual); + + use_old_geometry_for_boundary_conditions(); + }; + + fluid_int.assemble_equation_residual(0, [&]() { + set_bc::enforce_dirichlet_dofs_on_face(fluid_com, *fluid_face_, 0, nsd); + }, recover_interface_force_before_boundary_conditions, + restore_staged_geometry_after_boundary_conditions); + } + + return fluid_ok && !has_nan(fluid_sol); +} + +//---------------------------------------------------------------------- +// solve_solid — transfer recovered fluid force, solve solid. +// All-gathers local fluid force to global, transfers to global solid, +// then extracts this rank's local solid portion. +//---------------------------------------------------------------------- +bool PartitionedFSI::solve_solid() +{ + auto& solid_com = solid_sim_->com_mod; + auto& solid_int = solid_sim_->get_integrator(); + auto& solid_sol = solid_int.get_solutions(); + auto& cm = main_sim_->com_mod.cm; + auto& cm_mod = main_sim_->cm_mod; + + const int nsd = main_sim_->com_mod.nsd; + if (local_fluid_interface_force_.nrows() != nsd || + local_fluid_interface_force_.ncols() != fluid_face_->nNo) { + throw std::runtime_error( + "[PartitionedFSI] Fluid interface force is unavailable or has wrong shape."); + } + + // All-gather the residual-recovered fluid force to the global fluid face. + // The residual copy was already communicated before Dirichlet row enforcement, + // so duplicate partition-ring slots carry identical physical-node values. + auto global_fluid_force = gather_face_data(local_fluid_interface_force_, + fluid_face_global_nNo_, + fluid_face_local_offset_, + cm, cm_mod); + + // Transfer global fluid → global solid, then extract local solid portion. + // Only the owning (min-rank) slot of each node is read below, and that slot + // is the canonical one that transfer_data writes, so no dedup is needed. + auto global_solid_force = transfer_data(fluid_to_solid_map_, + global_fluid_force, + solid_face_global_nNo_); + const int nrows = global_solid_force.nrows(); + Array local_solid_force(nrows, solid_face_->nNo); + for (int a = 0; a < solid_face_->nNo; a++) { + // Apply the nodal force on the owning rank only; partition-ring nodes are + // shared by several ranks and all_fun::commu (called after the force + // callback in step_equation) SUMS R across them. Subtracting the full + // force on every sharing rank would multiply it by the node's multiplicity. + double w = solid_face_owner_[a] ? 1.0 : 0.0; + for (int i = 0; i < nrows; i++) + local_solid_force(i, a) = + w * global_solid_force(i, solid_face_local_offset_ + a); + } + + if (solid_interface_force_mask_.nrows() == nrows && + solid_interface_force_mask_.ncols() == solid_face_->nNo) { + for (int a = 0; a < solid_face_->nNo; a++) + for (int i = 0; i < nrows; i++) + local_solid_force(i, a) *= solid_interface_force_mask_(i, a); + } + + set_bc::set_bc_dir(solid_com, solid_sol); + solid_int.step_equation(0, [&]() { + fsi_coupling::apply_traction_on_solid( + solid_com, solid_com.eq[0], *solid_face_, local_solid_force); + }); + return !has_nan(solid_sol); +} + +//---------------------------------------------------------------------- +// solve_mesh — mesh equation with relaxed displacement, deform fluid mesh. +// disp_prev_ is global; extract this rank's local mesh face portion. +//---------------------------------------------------------------------- +bool PartitionedFSI::solve_mesh(const Array& x_ref, int mesh_s) +{ + auto& mesh_com = mesh_sim_->com_mod; + auto& mesh_int = mesh_sim_->get_integrator(); + auto& mesh_sol = mesh_int.get_solutions(); + const int nsd = main_sim_->com_mod.nsd; + + // Transfer global solid displacement → global mesh, extract local portion. + // transfer_data only writes the canonical target slot of each node; dedup + // fills the duplicate ring slots so every rank reads a correct Dirichlet + // value below (the displacement is applied on all face nodes, not just owners). + auto global_mesh_disp = transfer_data(solid_to_mesh_map_, disp_prev_, + mesh_face_global_nNo_); + apply_dedup(global_mesh_disp, mesh_face_canonical_); + Array local_mesh_disp(nsd, mesh_face_->nNo); + for (int a = 0; a < mesh_face_->nNo; a++) + for (int i = 0; i < nsd; i++) + local_mesh_disp(i, a) = global_mesh_disp(i, mesh_face_local_offset_ + a); + + set_bc::set_bc_dir(mesh_com, mesh_sol); + fsi_coupling::apply_displacement_on_mesh( + mesh_com, mesh_com.eq[0], *mesh_face_, local_mesh_disp, mesh_sol); + mesh_int.step_equation(0, [&]() { + set_bc::enforce_dirichlet_on_face(mesh_com, *mesh_face_, nsd); + }); + if (has_nan(mesh_sol)) return false; + + // Match the generalized-alpha fluid residual geometry used by monolithic + // FSI: fluid assembly sees x_ref + alpha_f * (Dn - Do). + update_fluid_mesh_coordinates(x_ref, mesh_s, fluid_sim_->com_mod.eq[0].af); + return true; +} + +//---------------------------------------------------------------------- +// update_fluid_mesh_coordinates +//---------------------------------------------------------------------- +void PartitionedFSI::update_fluid_mesh_coordinates( + const Array& x_ref, int mesh_s, double theta) +{ + auto& fluid_com = fluid_sim_->com_mod; + auto& mesh_sol = mesh_sim_->get_integrator().get_solutions(); + + fluid_com.x = fsi_coupling::staged_fluid_mesh_coordinates( + x_ref, + mesh_sol.current.get_displacement(), + mesh_sol.old.get_displacement(), + mesh_s, + main_sim_->com_mod.nsd, + theta); +} + +//====================================================================== +// step — one coupling iteration loop for one time step +//====================================================================== +bool PartitionedFSI::step() +{ + auto& fluid_com = fluid_sim_->com_mod; + auto& solid_com = solid_sim_->com_mod; + auto& mesh_com = mesh_sim_->com_mod; + auto& cm_mod = main_sim_->cm_mod; + auto& cm = main_sim_->com_mod.cm; + const int nsd = main_sim_->com_mod.nsd; + const int cTS = main_sim_->com_mod.cTS; + + auto& fluid_sol = fluid_sim_->get_integrator().get_solutions(); + auto& solid_sol = solid_sim_->get_integrator().get_solutions(); + auto& mesh_sol = mesh_sim_->get_integrator().get_solutions(); + + omega_ = config_.initial_relaxation; + r_prev_.clear(); + + // Save predictor state + struct SavedState { Array An, Yn, Dn; }; + auto save_state = [](SolutionStates& s) -> SavedState { + return {s.current.get_acceleration(), s.current.get_velocity(), s.current.get_displacement()}; + }; + auto restore_state = [](SolutionStates& s, const SavedState& st) { + s.current.get_acceleration() = st.An; + s.current.get_velocity() = st.Yn; + s.current.get_displacement() = st.Dn; + }; + SavedState fluid_pred = save_state(fluid_sol); + SavedState solid_pred = save_state(solid_sol); + SavedState mesh_pred = save_state(mesh_sol); + + // Save mesh coordinates at start of time step = x_original + Do + Array x_ref(fluid_com.x); + + // ALE mesh velocity from predictor (updated after each mesh solve) + const int mesh_s = mesh_com.eq[0].s; + Array mesh_vel_Yn(nsd, mesh_com.tnNo); + Array mesh_vel_Yo(nsd, mesh_com.tnNo); + { + auto& mYn = mesh_sol.current.get_velocity(); + auto& mYo = mesh_sol.old.get_velocity(); + for (int a = 0; a < mesh_com.tnNo; a++) + for (int i = 0; i < nsd; i++) { + mesh_vel_Yn(i, a) = mYn(mesh_s + i, a); + mesh_vel_Yo(i, a) = mYo(mesh_s + i, a); + } + } + + // Initial interface state from predictor — extract local, all-gather to global + Array disp_current; + { + auto local_disp = fsi_coupling::extract_solid_displacement( + solid_com, solid_com.eq[0], *solid_face_, solid_sol); + disp_prev_ = gather_face_data(local_disp, solid_face_global_nNo_, + solid_face_local_offset_, cm, cm_mod); + } + compute_interface_velocity(); + + bool converged = false; + + for (int cp = 0; cp < config_.max_coupling_iterations; cp++) { + + // Restore all sub-sims to predictor state + restore_state(fluid_sol, fluid_pred); + restore_state(solid_sol, solid_pred); + restore_state(mesh_sol, mesh_pred); + + // ---- 1. Mesh solve + deform fluid mesh ---- + // Use latest disp_prev_ (relaxed from previous iter, or predictor on iter 0). + // Writes fluid_com.x = x_ref + (Dn - Do) so the fluid solves on the deformed mesh. + if (!solve_mesh(x_ref, mesh_s)) { + if (cm.mas(cm_mod)) std::cout << " ABORT: NaN in mesh solve" << std::endl; + return false; + } + + // Update ALE mesh velocity from this iteration's mesh solve + { + auto& mYn = mesh_sol.current.get_velocity(); + for (int a = 0; a < mesh_com.tnNo; a++) + for (int i = 0; i < nsd; i++) + mesh_vel_Yn(i, a) = mYn(mesh_s + i, a); + } + + // ---- 2. Fluid solve ---- + if (!solve_fluid(x_ref, mesh_vel_Yo, mesh_vel_Yn)) { + if (cm.mas(cm_mod)) std::cout << " ABORT: NaN in fluid solve" << std::endl; + return false; + } + + // ---- 3. Solid solve ---- + if (!solve_solid()) { + if (cm.mas(cm_mod)) std::cout << " ABORT: NaN in solid solve" << std::endl; + return false; + } + + // ---- 4. Extract displacement (global), check convergence ---- + // Extract local solid displacement and all-gather to global so all ranks + // have identical arrays — no MPI reduction needed for norms. + { + auto local_disp = fsi_coupling::extract_solid_displacement( + solid_com, solid_com.eq[0], *solid_face_, solid_sol); + disp_current = gather_face_data(local_disp, solid_face_global_nNo_, + solid_face_local_offset_, cm, cm_mod); + } + + double res_norm = 0.0, disp_norm = 0.0; + for (int a = 0; a < solid_face_global_nNo_; a++) { + // Count each physical node once: duplicate ring slots (canonical != a) + // hold identical values and would otherwise be weighted by multiplicity, + // making the convergence metric depend on the partition count. + if (solid_face_canonical_[a] != a) continue; + for (int i = 0; i < nsd; i++) { + double res = disp_current(i, a) - disp_prev_(i, a); + res_norm += res * res; + disp_norm += disp_current(i, a) * disp_current(i, a); + } + } + res_norm = sqrt(res_norm); + disp_norm = sqrt(disp_norm); + double rel = (disp_norm > 1e-30) ? res_norm / disp_norm : res_norm; + + // ---- 5. Relaxation ---- + relax_interface(cp, nsd, disp_current); + compute_interface_velocity(); + + // Check for NaN/divergence (global arrays — consistent on all ranks) + { + bool bad = false; + double max_disp = 0; + for (int a = 0; a < solid_face_global_nNo_ && !bad; a++) + for (int i = 0; i < nsd; i++) { + if (std::isnan(disp_prev_(i, a)) || std::isinf(disp_prev_(i, a))) + { bad = true; break; } + max_disp = std::max(max_disp, std::abs(disp_prev_(i, a))); + } + if (bad || max_disp > 1e10) { + if (cm.mas(cm_mod)) std::cout << " ABORT: NaN/divergence after relaxation" << std::endl; + return false; + } + } + + // ---- 6. Output ---- + if (cp == 0) first_res_norm_ = res_norm; + int dB_val = 0; + double ri_r1 = 1.0; + if (first_res_norm_ > 1e-30 && res_norm > 0) { + ri_r1 = res_norm / first_res_norm_; + dB_val = static_cast(20.0 * log10(ri_r1)); + } + + if (cm.mas(cm_mod)) { + bool conv = rel < config_.coupling_tolerance; + bool saved = conv + && (cTS % fluid_sim_->com_mod.saveIncr == 0) + && (cTS >= fluid_sim_->com_mod.saveATS); + char buf[256]; + snprintf(buf, sizeof(buf), " CP %d-%d%s %10.3e %5d %10.3e %10.3e %10.3e %10.3e", + cTS, cp + 1, saved ? "s" : " ", + main_sim_->com_mod.timer.get_elapsed_time(), + dB_val, ri_r1, rel, omega_, disp_norm); + std::cout << buf << std::endl; + if (coupling_log_.is_open()) coupling_log_ << buf << std::endl; + if (histor_log_.is_open()) histor_log_ << buf << std::endl; + } + + if (rel < config_.coupling_tolerance) { converged = true; break; } + } + if (converged) { + // Keep saved results and the next time step's reference coordinates at + // the full n+1 mesh position after solving the fluid on n+alpha_f geometry. + update_fluid_mesh_coordinates(x_ref, mesh_s, 1.0); + } + return converged; +} + +//---------------------------------------------------------------------- +// save_results +//---------------------------------------------------------------------- +void PartitionedFSI::save_results() +{ + int cTS = main_sim_->com_mod.cTS; + Simulation* sims[3] = {fluid_sim_.get(), solid_sim_.get(), mesh_sim_.get()}; + + for (auto* sim : sims) { + auto& com = sim->com_mod; + auto& sol = sim->get_integrator().get_solutions(); + if (com.saveVTK) { + bool l2 = ((cTS % com.saveIncr) == 0); + bool l3 = (cTS >= com.saveATS); + if (l2 && l3) { + Array saved_fluid_x; + bool restore_fluid_x = false; + if (sim == fluid_sim_.get()) { + auto& mesh_x = mesh_sim_->com_mod.x; + if (mesh_x.nrows() != com.x.nrows() || + mesh_x.ncols() != com.x.ncols()) { + throw std::runtime_error( + "[PartitionedFSI] Fluid and mesh coordinate arrays are incompatible for output."); + } + saved_fluid_x = com.x; + com.x = mesh_x; + restore_fluid_x = true; + } + + try { + output::output_result(sim, com.timeP, 3, 0); + vtk_xml::write_vtus(sim, sol, false); + } catch (...) { + if (restore_fluid_x) { + com.x = saved_fluid_x; + } + throw; + } + if (restore_fluid_x) { + com.x = saved_fluid_x; + } + } + } + } +} diff --git a/Code/Source/solver/PartitionedFSI.h b/Code/Source/solver/PartitionedFSI.h new file mode 100644 index 000000000..8cd612e1f --- /dev/null +++ b/Code/Source/solver/PartitionedFSI.h @@ -0,0 +1,187 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef PARTITIONED_FSI_H +#define PARTITIONED_FSI_H + +#include "Simulation.h" +#include "Integrator.h" +#include "Array.h" +#include "CmMod.h" + +#include +#include +#include +#include + +/// @brief Coupling method for interface relaxation +enum class CouplingMethod { constant, aitken }; + +/// @brief Configuration for partitioned FSI coupling, read from XML input. +struct PartitionedFSIConfig { + int max_coupling_iterations = 50; + double coupling_tolerance = 1e-6; + double initial_relaxation = 1.0; + double omega_max = 1.0; + CouplingMethod coupling_method = CouplingMethod::aitken; + + // Face names for the FSI interface + std::string fluid_interface_face; + std::string solid_interface_face; +}; + +/// @brief Partitioned FSI coupling with 3 independent sub-Simulations. +/// +/// Each sub-field (fluid, struct, mesh) has its own Simulation object with +/// independent mesh, solution arrays, and linear system. No shared global +/// arrays, no DOF offsets (each eq.s=0), no regularization of inactive nodes. +/// +/// Implements Dirichlet-Neumann coupling with Aitken relaxation: +/// 1. Transfer solid displacement to mesh interface, solve mesh equation +/// 2. Deform fluid mesh using mesh displacement, solve fluid equation +/// 3. Recover fluid interface residual force, apply to solid, solve solid equation +/// 4. Extract solid displacement, apply Aitken relaxation +/// 5. Check coupling convergence +/// +/// Related to GitHub issue #431: Implement partitioned FSI in svMultiPhysics +class PartitionedFSI { +public: + PartitionedFSI(Simulation* main_simulation, const PartitionedFSIConfig& config, + const std::string& xml_file_path); + + ~PartitionedFSI(); + + void run(); + bool step(); + +private: + Simulation* main_sim_; + PartitionedFSIConfig config_; + std::string xml_file_path_; + + // Sub-simulations (owned) + std::unique_ptr fluid_sim_; + std::unique_ptr solid_sim_; + std::unique_ptr mesh_sim_; + + // Interface face pointers within each sub-sim + const faceType* fluid_face_ = nullptr; + const faceType* solid_face_ = nullptr; + const faceType* mesh_face_ = nullptr; + + // Mesh pointers within each sub-sim + const mshType* fluid_mesh_ = nullptr; + const mshType* solid_mesh_ = nullptr; + const mshType* mesh_mesh_ = nullptr; + + // Global face node counts and per-rank offsets (for MPI distribution) + int solid_face_global_nNo_ = 0; + int solid_face_local_offset_ = 0; + int fluid_face_global_nNo_ = 0; + int fluid_face_local_offset_ = 0; + int mesh_face_global_nNo_ = 0; + int mesh_face_local_offset_ = 0; + + // Node maps between interface faces: global_src_idx → global_tgt_idx + std::vector solid_to_fluid_map_; + std::vector fluid_to_solid_map_; + std::vector solid_to_mesh_map_; + + // Canonical slot maps: canonical_[i] = first global slot with same physical node as i + std::vector solid_face_canonical_; + std::vector fluid_face_canonical_; + std::vector mesh_face_canonical_; + + // Per-local-node owner flag for the solid interface face: 1 if THIS rank is + // the min-rank owner of the node, else 0. Used to apply the recovered + // interface force exactly once per physical node (it is summed by + // all_fun::commu). + std::vector solid_face_owner_; + + // Coupling state — indexed by GLOBAL solid face nodes (replicated on all ranks) + Array disp_prev_; + Array vel_prev_; + Array local_fluid_interface_force_; + Array solid_interface_force_mask_; + double omega_; + double first_res_norm_ = 0.0; + + // Aitken state — sized for global solid face + std::vector r_prev_; + + // Output files for coupling convergence history + std::ofstream coupling_log_; + std::ofstream histor_log_; + + void resolve_faces(); + void build_node_maps(); + void build_solid_interface_force_mask(); + + /// Solve fluid equation with current interface velocity and ALE mesh velocity + bool solve_fluid(const Array& fluid_x_old, + const Array& mesh_vel_Yo, + const Array& mesh_vel_Yn); + + /// Transfer recovered fluid interface force to solid, solve solid equation + bool solve_solid(); + + /// Solve mesh equation with relaxed displacement, deform fluid mesh + bool solve_mesh(const Array& x_ref, int mesh_s); + + /// Update fluid mesh coordinates with a staged mesh displacement increment + void update_fluid_mesh_coordinates(const Array& x_ref, int mesh_s, + double theta); + + /// Compute vel_prev_ (global) from disp_prev_ (global) using Newmark relationship + void compute_interface_velocity(); + + void relax_interface(int cp, int nsd, const Array& disp_current); + void relax_constant(int cp, int nsd, const Array& disp_current); + void relax_aitken(int cp, int nsd, const Array& disp_current); + + /// Compute global_nNo and local_offset for a distributed face + static void compute_face_global_info(const faceType& face, cmType& cm, + const CmMod& cm_mod, + int& global_nNo, int& local_offset); + + /// All-gather local face data (nrows, local_nNo) to global (nrows, global_nNo) + static Array gather_face_data(const Array& local_data, + int global_nNo, int local_offset, + cmType& cm, const CmMod& cm_mod); + + /// All-gather local (local_src → global_tgt) map to global (global_src → global_tgt) map + static void gather_global_map(const std::vector& local_map, + int global_src_nNo, + cmType& cm, const CmMod& cm_mod, + std::vector& global_map); + + /// Build local face_a → global face_b node map using pre-gathered global face_b coords + static void build_face_node_map(const faceType& face_a, const ComMod& com_a, + int global_b_nNo, const Array& global_b_coords, + std::vector& a_to_global_b); + + /// Build a minimal sub-simulation XML for the given role by extracting + /// the tagged equation and its meshes from the main XML. Returns temp file path. + static std::string build_sub_xml(const std::string& main_xml_path, + const std::string& role); + + /// Transfer data from global src face to global tgt face using global map + static Array transfer_data(const std::vector& src_to_tgt_map, + const Array& src_data, int tgt_nNo); + + /// Build canonical[i] = index of first occurrence of global node ID i across ranks + std::vector build_face_dedup_map(const faceType& face, const ComMod& com, + int global_nNo, cmType& cm, const CmMod& cm_mod); + + /// Propagate canonical slot values to duplicate slots (canonical[i] != i) + static void apply_dedup(Array& arr, const std::vector& canonical); + + /// Build per-local-node flags marking whether THIS rank is the min-rank + /// owner of each face node (1 = owner, 0 = shared but owned by a lower rank). + static std::vector build_face_owner_flags(const faceType& face, + const ComMod& com, cmType& cm); + + void save_results(); +}; + +#endif // PARTITIONED_FSI_H diff --git a/Code/Source/solver/Simulation.cpp b/Code/Source/solver/Simulation.cpp index 7e6d56523..c1ce70250 100644 --- a/Code/Source/solver/Simulation.cpp +++ b/Code/Source/solver/Simulation.cpp @@ -3,6 +3,7 @@ #include "Simulation.h" #include "Integrator.h" +#include "PartitionedFSI.h" #include "all_fun.h" #include "load_msh.h" @@ -12,6 +13,17 @@ #include #include +void add_eq_linear_algebra(ComMod& com_mod, eqType& lEq) +{ + lEq.linear_algebra = LinearAlgebraFactory::create_interface(lEq.linear_algebra_type); + lEq.linear_algebra->set_preconditioner(lEq.linear_algebra_preconditioner); + lEq.linear_algebra->initialize(com_mod, lEq); + + if (lEq.linear_algebra_assembly_type != consts::LinearAlgebraType::none) { + lEq.linear_algebra->set_assembly(lEq.linear_algebra_assembly_type); + } +} + Simulation::Simulation() { roInf = 0.2; @@ -40,6 +52,13 @@ void Simulation::read_parameters(const std::string& file_name) parameters.read_xml(file_name); } +/// @brief Read solver parameters from an in-memory XML string (no file). +// +void Simulation::read_parameters_from_string(const std::string& xml_content) +{ + parameters.read_xml_from_string(xml_content); +} + /// @brief Set the simulation and module member data. /// /// Replicates the README subroutine lines to set COMMOD module varliables @@ -112,3 +131,81 @@ Integrator& Simulation::get_integrator() } return *integrator_; } + +/// @brief Get pointer to PartitionedFSI object (null if not configured) +PartitionedFSI* Simulation::get_partitioned_fsi() +{ + return partitioned_fsi_.get(); +} + +/// @brief Initialize partitioned FSI if configured in parameters. +/// +/// Parameters are only parsed on rank 0 (slaves skip read_files), so we +/// broadcast the active flag and config to all ranks before branching. +void Simulation::initialize_partitioned_fsi(const std::string& xml_file_path) +{ + auto& cm = com_mod.cm; + auto& cm_mod_ref = cm_mod; + + // Rank 0 determines whether partitioned FSI is active and builds the config. + // Broadcast the decision so all ranks take the same path. + int active = 0; + PartitionedFSIConfig config; + + if (cm.mas(cm_mod_ref)) { + auto& pcp = parameters.partitioned_coupling_parameters; + // Active when is present and at least one equation + // carries a partitioned role attribute. + if (pcp.defined()) { + for (auto* ep : parameters.equation_parameters) { + if (ep->role.defined() && !ep->role.value().empty()) { active = 1; break; } + } + } + if (active) { + config.max_coupling_iterations = pcp.max_coupling_iterations.value(); + config.coupling_tolerance = pcp.coupling_tolerance.value(); + config.initial_relaxation = pcp.initial_relaxation.value(); + config.omega_max = pcp.omega_max.value(); + + std::string method = pcp.coupling_method.value(); + if (method == "constant") config.coupling_method = CouplingMethod::constant; + else if (method == "aitken") config.coupling_method = CouplingMethod::aitken; + else throw std::runtime_error("[PartitionedFSI] Unknown Coupling_method: " + method); + + config.fluid_interface_face = pcp.fluid_interface_face.value(); + config.solid_interface_face = pcp.solid_interface_face.value(); + } + } + + // Broadcast the active flag and config fields to all ranks. + MPI_Bcast(&active, 1, MPI_INT, 0, cm.com()); + if (!active) return; + + int max_iter = config.max_coupling_iterations; + double tol = config.coupling_tolerance; + double relax = config.initial_relaxation; + double omax = config.omega_max; + int method_i = static_cast(config.coupling_method); + MPI_Bcast(&max_iter, 1, MPI_INT, 0, cm.com()); + MPI_Bcast(&tol, 1, MPI_DOUBLE, 0, cm.com()); + MPI_Bcast(&relax, 1, MPI_DOUBLE, 0, cm.com()); + MPI_Bcast(&omax, 1, MPI_DOUBLE, 0, cm.com()); + MPI_Bcast(&method_i, 1, MPI_INT, 0, cm.com()); + + auto bcast_str = [&](std::string& s) { + int len = static_cast(s.size()); + MPI_Bcast(&len, 1, MPI_INT, 0, cm.com()); + s.resize(len); + MPI_Bcast(s.data(), len, MPI_CHAR, 0, cm.com()); + }; + bcast_str(config.fluid_interface_face); + bcast_str(config.solid_interface_face); + + config.max_coupling_iterations = max_iter; + config.coupling_tolerance = tol; + config.initial_relaxation = relax; + config.omega_max = omax; + config.coupling_method = static_cast(method_i); + + partitioned_fsi_ = std::make_unique(this, config, xml_file_path); +} diff --git a/Code/Source/solver/Simulation.h b/Code/Source/solver/Simulation.h index a8e092158..f38d7751f 100644 --- a/Code/Source/solver/Simulation.h +++ b/Code/Source/solver/Simulation.h @@ -13,8 +13,11 @@ #include #include -// Forward declaration +// Forward declarations class Integrator; +class PartitionedFSI; + +void add_eq_linear_algebra(ComMod& com_mod, eqType& lEq); class Simulation { @@ -28,6 +31,8 @@ class Simulation { ChnlMod& get_chnl_mod() { return chnl_mod; }; ComMod& get_com_mod() { return com_mod; }; Integrator& get_integrator(); + PartitionedFSI* get_partitioned_fsi(); + void initialize_partitioned_fsi(const std::string& xml_file_path); // Initialize the Integrator object after simulation setup is complete // Takes ownership of solution states via move semantics @@ -36,6 +41,9 @@ class Simulation { // Read a solver paramerer input XML file. void read_parameters(const std::string& fileName); + // Read solver parameters from an in-memory XML string (no file on disk). + void read_parameters_from_string(const std::string& xml_content); + // Set simulation and module member data from Parameters. void set_module_parameters(); @@ -75,6 +83,9 @@ class Simulation { private: // Time integrator for Newton iteration loop std::unique_ptr integrator_; + + // Partitioned FSI coupling (null if not configured) + std::unique_ptr partitioned_fsi_; }; #endif diff --git a/Code/Source/solver/SimulationLogger.h b/Code/Source/solver/SimulationLogger.h index 4023079bb..03dd96b44 100755 --- a/Code/Source/solver/SimulationLogger.h +++ b/Code/Source/solver/SimulationLogger.h @@ -33,6 +33,8 @@ class SimulationLogger { bool is_initialized() const { return log_file_.is_open(); } + void set_cout_write(bool v) const { cout_write_ = v; } + ~SimulationLogger() { log_file_.close(); diff --git a/Code/Source/solver/eq_assem.cpp b/Code/Source/solver/eq_assem.cpp index 9cea03b30..39c69719f 100644 --- a/Code/Source/solver/eq_assem.cpp +++ b/Code/Source/solver/eq_assem.cpp @@ -58,9 +58,13 @@ void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& cDmn = all_fun::domain(com_mod, msh, cEq, Ec); auto cPhys = eq.dmn[cDmn].phys; - Vector ptr(eNoN); - Vector N(eNoN), hl(eNoN); - Array yl(tDof,eNoN), lR(dof,eNoN); + // For ALE partitioned FSI, extend yl to hold mesh velocity + const bool has_ale = (com_mod.ale_mesh_velocity.size() > 0); + const int yl_nrows = has_ale ? tDof + nsd : tDof; + + Vector ptr(eNoN); + Vector N(eNoN), hl(eNoN); + Array yl(yl_nrows,eNoN), lR(dof,eNoN); Array3 lK(dof*dof,eNoN,eNoN); for (int a = 0; a < eNoN; a++) { @@ -70,6 +74,11 @@ void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& for (int i = 0; i < tDof; i++) { yl(i,a) = Yg(i,Ac); } + if (has_ale) { + for (int i = 0; i < nsd; i++) { + yl(nsd+1+i, a) = com_mod.ale_mesh_velocity(i, Ac); + } + } } // Updating the shape functions, if neccessary @@ -87,7 +96,7 @@ void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& N = lFa.N.col(g); double h = 0.0; - Vector y(tDof); + Vector y(yl_nrows); for (int a = 0; a < eNoN; a++) { h = h + N(a)*hl(a); @@ -162,6 +171,9 @@ void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& /// @param Dg void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const Vector& hg, const SolutionStates& solutions) { + // Local alias for old displacement + const auto& Do = solutions.old.get_displacement(); + using namespace consts; using namespace utils; const auto& Dg = solutions.intermediate.get_displacement(); @@ -296,6 +308,10 @@ void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const /// an arbitrary vector. void fsi_ls_upd(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const SolutionStates& solutions) { + // Local aliases for displacement arrays + const auto& Dn = solutions.current.get_displacement(); + const auto& Do = solutions.old.get_displacement(); + using namespace consts; using namespace utils; using namespace fsi_linear_solver; diff --git a/Code/Source/solver/fluid.cpp b/Code/Source/solver/fluid.cpp index 0d58a2725..1a21d0449 100644 --- a/Code/Source/solver/fluid.cpp +++ b/Code/Source/solver/fluid.cpp @@ -51,7 +51,7 @@ void b_fluid(ComMod& com_mod, const int eNoN, const double w, const Vector u(nsd); - if (com_mod.mvMsh) { + if (com_mod.mvMsh || com_mod.ale_mesh_velocity.size() > 0) { for (int i = 0; i < nsd; i++) { int j = i + nsd + 1; u(i) = y(i) - y(j); @@ -163,7 +163,7 @@ void bw_fluid_2d(ComMod& com_mod, const int eNoNw, const int eNoNq, const double Vector uh(2); - if (com_mod.mvMsh) { + if (com_mod.mvMsh || com_mod.ale_mesh_velocity.size() > 0) { for (int a = 0; a < eNoNw; a++) { uh(0) = uh(0) + Nw(a)*yl(3,a); uh(1) = uh(1) + Nw(a)*yl(4,a); @@ -319,7 +319,7 @@ void bw_fluid_3d(ComMod& com_mod, const int eNoNw, const int eNoNq, const double Vector uh(3); - if (com_mod.mvMsh) { + if (com_mod.mvMsh || com_mod.ale_mesh_velocity.size() > 0) { for (int a = 0; a < eNoNw; a++) { uh(0) = uh(0) + Nw(a)*yl(4,a); uh(1) = uh(1) + Nw(a)*yl(5,a); @@ -522,14 +522,18 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const SolutionStates& s #endif // FLUID: dof = nsd+1 - Vector ptr(eNoN); - Array xl(nsd,eNoN); - + // For ALE partitioned FSI, extend yl to hold mesh velocity at DOFs nsd+1..2*nsd + const bool has_ale = (com_mod.ale_mesh_velocity.size() > 0); + const int yl_nrows = has_ale ? tDof + nsd : tDof; + + Vector ptr(eNoN); + Array xl(nsd,eNoN); + // local acceleration vector (for a single element) Array al(tDof,eNoN); - + // local velocity vector (for a single element) - Array yl(tDof,eNoN); + Array yl(yl_nrows,eNoN); Array bfl(nsd,eNoN); // local (weak form) residual vector (for a single element) @@ -569,10 +573,16 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const SolutionStates& s xl(i,a) = com_mod.x(i,Ac); bfl(i,a) = com_mod.Bf(i,Ac); } - for (int i = 0; i < al.nrows(); i++) { + for (int i = 0; i < tDof; i++) { al(i,a) = Ag(i,Ac); yl(i,a) = Yg(i,Ac); } + // Append ALE mesh velocity at DOFs nsd+1..2*nsd (same layout as mvMsh) + if (has_ale) { + for (int i = 0; i < nsd; i++) { + yl(nsd+1+i, a) = com_mod.ale_mesh_velocity(i, Ac); + } + } } // Initialize residual and tangents @@ -895,7 +905,7 @@ void fluid_2d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e } // Update convection velocity relative to mesh velocity - if (com_mod.mvMsh) { + if (com_mod.mvMsh || com_mod.ale_mesh_velocity.size() > 0) { for (int a = 0; a < eNoNw; a++) { u(0) = u(0) - Nw(a)*yl(3,a); u(1) = u(1) - Nw(a)*yl(4,a); @@ -1212,7 +1222,7 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e } // Update convection velocity relative to mesh velocity - if (com_mod.mvMsh) { + if (com_mod.mvMsh || com_mod.ale_mesh_velocity.size() > 0) { for (int a = 0; a < eNoNw; a++) { u(0) = u(0) - Nw(a)*yl(3,a); u(1) = u(1) - Nw(a)*yl(4,a); @@ -1562,7 +1572,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // Update convection velocity relative to mesh velocity // - if (com_mod.mvMsh) { + if (com_mod.mvMsh || com_mod.ale_mesh_velocity.size() > 0) { for (int a = 0; a < eNoNw; a++) { u[0] = u[0] - Nw(a)*yl(4,a); u[1] = u[1] - Nw(a)*yl(5,a); @@ -1906,7 +1916,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // Update convection velocity relative to mesh velocity // - if (com_mod.mvMsh) { + if (com_mod.mvMsh || com_mod.ale_mesh_velocity.size() > 0) { for (int a = 0; a < eNoNw; a++) { u[0] = u[0] - Nw(a)*yl(4,a); u[1] = u[1] - Nw(a)*yl(5,a); diff --git a/Code/Source/solver/fsi.h b/Code/Source/solver/fsi.h index def8ed8fb..49fd25c87 100644 --- a/Code/Source/solver/fsi.h +++ b/Code/Source/solver/fsi.h @@ -15,4 +15,3 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So }; #endif - diff --git a/Code/Source/solver/fsi_coupling.cpp b/Code/Source/solver/fsi_coupling.cpp new file mode 100644 index 000000000..8f6c01fcf --- /dev/null +++ b/Code/Source/solver/fsi_coupling.cpp @@ -0,0 +1,191 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +#include "fsi_coupling.h" +#include "Integrator.h" + +#include +#include + +namespace fsi_coupling { + +//---------------------------------------------------------------------- +// extract_solid_displacement +//---------------------------------------------------------------------- +Array extract_solid_displacement( + const ComMod& com_mod, const eqType& solid_eq, + const faceType& lFa, const SolutionStates& solutions) +{ + const int nsd = com_mod.nsd; + const int s = solid_eq.s; // DOF offset for the solid equation + const auto& Dn = solutions.current.get_displacement(); + + Array result(nsd, lFa.nNo); + for (int a = 0; a < lFa.nNo; a++) { + int Ac = lFa.gN(a); + for (int i = 0; i < nsd; i++) { + result(i, a) = Dn(i + s, Ac); + } + } + return result; +} + +//---------------------------------------------------------------------- +// copy_time_integration_parameters +//---------------------------------------------------------------------- +void copy_time_integration_parameters(const eqType& source_eq, + eqType& target_eq) +{ + target_eq.roInf = source_eq.roInf; + target_eq.am = source_eq.am; + target_eq.af = source_eq.af; + target_eq.beta = source_eq.beta; + target_eq.gam = source_eq.gam; +} + +//---------------------------------------------------------------------- +// apply_velocity_on_fluid +//---------------------------------------------------------------------- +void apply_velocity_on_fluid( + ComMod& com_mod, const eqType& fluid_eq, + const faceType& lFa, + const Array& velocity, + SolutionStates& solutions) +{ + const int nsd = com_mod.nsd; + const int s = fluid_eq.s; + const double dt = com_mod.dt; + const double gam = fluid_eq.gam; + + auto& An = solutions.current.get_acceleration(); + auto& Yn = solutions.current.get_velocity(); + const auto& Yo = solutions.old.get_velocity(); + const auto& Ao = solutions.old.get_acceleration(); + + for (int a = 0; a < lFa.nNo; a++) { + int Ac = lFa.gN(a); + for (int i = 0; i < nsd; i++) { + Yn(i + s, Ac) = velocity(i, a); + double a_new; + newmark::state_from_velocity( + velocity(i, a), Yo(i + s, Ac), Ao(i + s, Ac), dt, gam, a_new); + An(i + s, Ac) = a_new; + } + } +} + +//---------------------------------------------------------------------- +// extract_fluid_residual_force +//---------------------------------------------------------------------- +Array extract_fluid_residual_force( + const ComMod& com_mod, + const faceType& fluid_face, + const Array& residual) +{ + const int nsd = com_mod.nsd; + if (residual.nrows() < nsd || residual.ncols() != com_mod.tnNo) { + throw std::runtime_error( + "[fsi_coupling] Fluid residual shape does not match fluid simulation."); + } + + Array force(nsd, fluid_face.nNo); + for (int a = 0; a < fluid_face.nNo; a++) { + int Ac = fluid_face.gN(a); + for (int i = 0; i < nsd; i++) { + force(i, a) = -residual(i, Ac); + } + } + return force; +} + +//---------------------------------------------------------------------- +// apply_traction_on_solid +//---------------------------------------------------------------------- +void apply_traction_on_solid( + ComMod& com_mod, const eqType& solid_eq, + const faceType& lFa, + const Array& force) +{ + // The array contains consistent nodal forces (external force on solid). + // In svMultiPhysics, external forces are SUBTRACTED from R (for example: + // lR -= w*N*h). So R -= force. + for (int a = 0; a < lFa.nNo; a++) { + int Ac = lFa.gN(a); + for (int i = 0; i < force.nrows(); i++) { + com_mod.R(i, Ac) -= force(i, a); + } + } +} + +//---------------------------------------------------------------------- +// apply_displacement_on_mesh +//---------------------------------------------------------------------- +void apply_displacement_on_mesh( + ComMod& com_mod, const eqType& mesh_eq, + const faceType& lFa, + const Array& displacement, + SolutionStates& solutions) +{ + const int nsd = com_mod.nsd; + const int s = mesh_eq.s; + const double dt = com_mod.dt; + const double gam = mesh_eq.gam; + const double beta = mesh_eq.beta; + + auto& An = solutions.current.get_acceleration(); + auto& Yn = solutions.current.get_velocity(); + auto& Dn = solutions.current.get_displacement(); + const auto& Do = solutions.old.get_displacement(); + const auto& Yo = solutions.old.get_velocity(); + const auto& Ao = solutions.old.get_acceleration(); + + for (int a = 0; a < lFa.nNo; a++) { + int Ac = lFa.gN(a); + for (int i = 0; i < nsd; i++) { + Dn(i + s, Ac) = displacement(i, a); + double a_new, v_new; + newmark::state_from_displacement( + displacement(i, a), Do(i + s, Ac), Yo(i + s, Ac), Ao(i + s, Ac), + dt, beta, gam, a_new, v_new); + An(i + s, Ac) = a_new; + Yn(i + s, Ac) = v_new; + } + } +} + +//---------------------------------------------------------------------- +// staged_fluid_mesh_coordinates +//---------------------------------------------------------------------- +Array staged_fluid_mesh_coordinates( + const Array& x_ref, + const Array& mesh_Dn, + const Array& mesh_Do, + int mesh_s, + int nsd, + double theta) +{ + if (nsd <= 0 || mesh_s < 0 || !std::isfinite(theta)) { + throw std::runtime_error( + "[fsi_coupling] Invalid inputs for staged fluid mesh coordinates."); + } + if (x_ref.nrows() != nsd || + mesh_Dn.ncols() != x_ref.ncols() || + mesh_Do.ncols() != x_ref.ncols() || + mesh_Dn.nrows() < mesh_s + nsd || + mesh_Do.nrows() < mesh_s + nsd) { + throw std::runtime_error( + "[fsi_coupling] Mesh displacement shapes do not match fluid coordinates."); + } + + Array x_stage(nsd, x_ref.ncols()); + for (int a = 0; a < x_ref.ncols(); a++) { + for (int i = 0; i < nsd; i++) { + x_stage(i, a) = x_ref(i, a) + + theta * (mesh_Dn(i + mesh_s, a) + - mesh_Do(i + mesh_s, a)); + } + } + return x_stage; +} + +} // namespace fsi_coupling diff --git a/Code/Source/solver/fsi_coupling.h b/Code/Source/solver/fsi_coupling.h new file mode 100644 index 000000000..e7014a444 --- /dev/null +++ b/Code/Source/solver/fsi_coupling.h @@ -0,0 +1,112 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef FSI_COUPLING_H +#define FSI_COUPLING_H + +#include "ComMod.h" +#include "SolutionStates.h" + +/// @brief FSI interface data exchange functions for partitioned coupling. +/// +/// These functions apply residual-recovered fluid nodal forces and transfer +/// solid displacement/velocity data at the FSI interface, enabling partitioned +/// (Dirichlet-Neumann) coupling between separately solved fluid and solid +/// equations. +/// +/// Related to GitHub issue #431: Implement partitioned FSI in svMultiPhysics + +namespace fsi_coupling { + +/// @brief Extract solid displacement at interface face nodes. +/// +/// @param com_mod Common module +/// @param solid_eq The solid equation (for DOF offset) +/// @param solid_face The solid-side FSI interface face +/// @param solutions Solution states +/// @return Array(nsd, solid_face.nNo) of displacement values +Array extract_solid_displacement( + const ComMod& com_mod, const eqType& solid_eq, + const faceType& solid_face, const SolutionStates& solutions); + +/// @brief Copy generalized-alpha/Newmark coefficients from one equation to another. +/// +/// Partitioned FSI uses a standalone structural subproblem, but to match the +/// monolithic FSI structural domain it must use the same time-integration +/// parameters as the FSI/fluid equation. +void copy_time_integration_parameters(const eqType& source_eq, + eqType& target_eq); + +/// @brief Apply velocity as strong Dirichlet BC on fluid interface nodes. +/// Directly sets Yn at the fluid equation DOF range for the face nodes. +void apply_velocity_on_fluid( + ComMod& com_mod, const eqType& fluid_eq, + const faceType& fluid_face, + const Array& velocity, + SolutionStates& solutions); + +/// @brief Extract force ON the solid from the unconstrained fluid residual. +/// +/// The returned array contains the velocity-row residual at the fluid interface +/// converted to the external-force-on-solid convention expected by +/// apply_traction_on_solid(). +/// +/// @param com_mod Common module for the fluid subproblem +/// @param fluid_face The fluid-side FSI interface face +/// @param residual Unconstrained assembled fluid residual +/// @return Array(nsd, fluid_face.nNo) of consistent nodal forces +Array extract_fluid_residual_force( + const ComMod& com_mod, + const faceType& fluid_face, + const Array& residual); + +/// @brief Apply pre-computed consistent nodal forces to the solid residual. +/// +/// The force array is interpreted as the external force ON the solid. It is +/// subtracted from com_mod.R at the global node locations corresponding to the +/// solid face. This should be called during the post-assembly callback of +/// step_equation() for the solid equation. +/// +/// @param com_mod Common module (R is modified) +/// @param solid_eq The solid equation (for DOF offset) +/// @param solid_face The solid-side FSI interface face +/// @param force Array(nsd, solid_face.nNo) of consistent nodal forces +void apply_traction_on_solid( + ComMod& com_mod, const eqType& solid_eq, + const faceType& solid_face, + const Array& force); + +/// @brief Apply displacement as strong Dirichlet BC on mesh interface nodes. +/// +/// Directly sets the displacement in the solution arrays (An, Yn, Dn) for +/// the mesh equation DOF range at the interface face nodes. +/// +/// @param com_mod Common module +/// @param mesh_eq The mesh equation (for DOF offset) +/// @param mesh_face The mesh-side FSI interface face +/// @param displacement Array(nsd, mesh_face.nNo) of displacement values +/// @param solutions Solution states (modified) +void apply_displacement_on_mesh( + ComMod& com_mod, const eqType& mesh_eq, + const faceType& mesh_face, + const Array& displacement, + SolutionStates& solutions); + +/// @brief Compute staged ALE fluid mesh coordinates from mesh displacement. +/// +/// The returned coordinates are +/// x_stage = x_ref + theta * (mesh_Dn - mesh_Do), using the mesh equation +/// displacement block starting at mesh_s. For a generalized-alpha fluid solve, +/// theta should be the fluid equation alpha_f; before output/old-state update, +/// theta should be 1.0. +Array staged_fluid_mesh_coordinates( + const Array& x_ref, + const Array& mesh_Dn, + const Array& mesh_Do, + int mesh_s, + int nsd, + double theta); + +} // namespace fsi_coupling + +#endif // FSI_COUPLING_H diff --git a/Code/Source/solver/main.cpp b/Code/Source/solver/main.cpp index 4cfb35fb4..19bcea7dc 100644 --- a/Code/Source/solver/main.cpp +++ b/Code/Source/solver/main.cpp @@ -9,6 +9,7 @@ // #include "Simulation.h" #include "Integrator.h" +#include "PartitionedFSI.h" #include "all_fun.h" #include "bf.h" @@ -34,23 +35,7 @@ #include #include #include -#include -//------------------------ -// add_eq_linear_algebra -//------------------------ -// Create a LinearAlgebra object for an equation. -// -void add_eq_linear_algebra(ComMod& com_mod, eqType& lEq) -{ - lEq.linear_algebra = LinearAlgebraFactory::create_interface(lEq.linear_algebra_type); - lEq.linear_algebra->set_preconditioner(lEq.linear_algebra_preconditioner); - lEq.linear_algebra->initialize(com_mod, lEq); - - if (lEq.linear_algebra_assembly_type != consts::LinearAlgebraType::none) { - lEq.linear_algebra->set_assembly(lEq.linear_algebra_assembly_type); - } -} void finalize_linear_algebra(eqType& lEq) { @@ -350,6 +335,8 @@ void iterate_solution(Simulation* simulation) #endif int iEqOld = cEq; + + // Monolithic Newton iteration loop integrator.step(); #ifdef debug_iterate_solution @@ -572,7 +559,12 @@ void iterate_solution(Simulation* simulation) void run_simulation(Simulation* simulation) { - iterate_solution(simulation); + auto* partitioned_fsi = simulation->get_partitioned_fsi(); + if (partitioned_fsi) { + partitioned_fsi->run(); + } else { + iterate_solution(simulation); + } } @@ -658,6 +650,9 @@ int main(int argc, char *argv[]) add_eq_linear_algebra(simulation->com_mod, eq); } + // Initialize partitioned FSI coupling if configured + simulation->initialize_partitioned_fsi(file_name); + #ifdef debug_main for (int iM = 0; iM < simulation->com_mod.nMsh; iM++) { dmsg << "---------- iM " << iM; diff --git a/Code/Source/solver/mesh.cpp b/Code/Source/solver/mesh.cpp index b7b8cdfb9..c3d6c38cf 100644 --- a/Code/Source/solver/mesh.cpp +++ b/Code/Source/solver/mesh.cpp @@ -45,9 +45,12 @@ void construct_mesh(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const S auto& pSa = com_mod.pSa; bool pstEq = com_mod.pstEq; - // Start and end DOF - int is = nsd + 1; - int ie = 2*nsd; + // Start and end DOF for mesh equation in the global solution arrays. + // Use eq.s (DOF offset) instead of hardcoded nsd+1, so this works both + // in monolithic FSI (where mesh is the 3rd equation) and in partitioned + // FSI (where mesh is the only equation with eq.s=0). + int is = eq.s; + int ie = eq.s + nsd - 1; int eNoN = lM.eNoN; #ifdef debug_construct_mesh dmsg << "cEq: " << cEq; diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 63d9ba918..8fb56fdc7 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1829,7 +1829,7 @@ void read_eq(Simulation* simulation, EquationParameters* eq_params, eqType& lEq) // The many global COMMOD values set in the Fortan subroutine are set in the // Simulation::set_module_parameters() method. // -void read_files(Simulation* simulation, const std::string& file_name) +void read_files(Simulation* simulation, const std::string& file_name, bool from_string) { using namespace consts; @@ -1841,12 +1841,18 @@ void read_files(Simulation* simulation, const std::string& file_name) dmsg.banner(); #endif - // Read the solver XML file. + // Read the solver parameters. 'file_name' is either a path to an XML file or, + // when from_string is true, the XML content itself (used by the partitioned + // FSI driver to build sub-simulations without temp files). #ifdef debug_read_files - dmsg << "Read the solver XML file " << " ... "; + dmsg << "Read the solver XML " << " ... "; #endif if (!com_mod.resetSim) { - simulation->read_parameters(std::string(file_name)); + if (from_string) { + simulation->read_parameters_from_string(file_name); + } else { + simulation->read_parameters(std::string(file_name)); + } } auto& chnl_mod = simulation->get_chnl_mod(); @@ -1994,12 +2000,18 @@ void read_files(Simulation* simulation, const std::string& file_name) } } - if (eq.phys == EquationType::phys_mesh) { + if (eq.phys == EquationType::phys_mesh) { + // For partitioned FSI, mvMsh is set when Partitioned_coupling is configured + if (!com_mod.mvMsh && simulation->parameters.partitioned_coupling_parameters.defined()) { + com_mod.mvMsh = true; + } if (!com_mod.mvMsh) { - throw std::runtime_error("mesh equation can only be specified after FSI equation"); + throw std::runtime_error("mesh equation can only be specified after FSI or with Partitioned_coupling"); + } + if (com_mod.nEq > 0 && com_mod.eq[0].phys == EquationType::phys_FSI) { + // Use the explicit geometry coupling flag of the FSI equation. + eq.expl_geom_cpl = com_mod.eq[0].expl_geom_cpl; } - // Use the explicit geometry coupling flag of the FSI equation. - eq.expl_geom_cpl = com_mod.eq[0].expl_geom_cpl; } } #ifdef debug_read_files diff --git a/Code/Source/solver/read_files.h b/Code/Source/solver/read_files.h index 197f83fdb..8b698c86d 100644 --- a/Code/Source/solver/read_files.h +++ b/Code/Source/solver/read_files.h @@ -33,7 +33,7 @@ namespace read_files_ns { void read_eq(Simulation* simulation, EquationParameters* params, eqType& eq); - void read_files(Simulation* simulation, const std::string& file_name); + void read_files(Simulation* simulation, const std::string& file_name, bool from_string = false); void read_ls(Simulation* simulation, EquationParameters* eq_params, consts::SolverType solver_type, eqType& lEq); diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 65f744727..78a95bf85 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -1185,6 +1185,9 @@ void set_bc_dir_l(ComMod& com_mod, const bcType& lBc, const faceType& lFa, Array // void set_bc_dir_w(ComMod& com_mod, const SolutionStates& solutions) { + // Local alias for old displacement + const auto& Do = solutions.old.get_displacement(); + using namespace consts; const int cEq = com_mod.cEq; @@ -1307,8 +1310,12 @@ void set_bc_dir_wl(ComMod& com_mod, const bcType& lBc, const mshType& lM, const } } - Vector ptr(eNoN); - Array xl(nsd,eNoN), yl(tDof,eNoN), lR(dof,eNoN); + // For ALE partitioned FSI, extend yl to hold mesh velocity + const bool has_ale = (com_mod.ale_mesh_velocity.size() > 0); + const int yl_nrows = has_ale ? tDof + nsd : tDof; + + Vector ptr(eNoN); + Array xl(nsd,eNoN), yl(yl_nrows,eNoN), lR(dof,eNoN); Array3 lK(dof*dof,eNoN,eNoN); Array xbl(nsd,eNoNb), ubl(nsd,eNoNb); @@ -1337,6 +1344,11 @@ void set_bc_dir_wl(ComMod& com_mod, const bcType& lBc, const mshType& lM, const for (int i = 0; i < tDof; i++) { yl(i,a) = Yg(i,Ac); } + if (has_ale) { + for (int i = 0; i < nsd; i++) { + yl(nsd+1+i, a) = com_mod.ale_mesh_velocity(i, Ac); + } + } for (int i = 0; i < nsd; i++) { xl(i,a) = com_mod.x(i,Ac); @@ -2101,6 +2113,71 @@ void set_bc_undef_neu_l(ComMod& com_mod, const bcType& lBc, const faceType& lFa) } } +//---------------------------------------------------------------------- +// enforce_dirichlet_on_face +//---------------------------------------------------------------------- +void enforce_dirichlet_on_face(ComMod& com_mod, const faceType& lFa, int nsd) +{ + const auto& eq = com_mod.eq[com_mod.cEq]; + const int dof = eq.dof; + const auto& rowPtr = com_mod.rowPtr; + const auto& colPtr = com_mod.colPtr; + auto& R = com_mod.R; + auto& Val = com_mod.Val; + + for (int a = 0; a < lFa.nNo; a++) { + int rowN = lFa.gN(a); + for (int i = 0; i < dof; i++) { + R(i, rowN) = 0.0; + } + for (int j = rowPtr(rowN); j <= rowPtr(rowN + 1) - 1; j++) { + int colN = colPtr(j); + for (int iDof = 0; iDof < dof * dof; iDof++) { + Val(iDof, j) = 0.0; + } + if (colN == rowN) { + for (int i = 0; i < dof; i++) { + Val(i * dof + i, j) = 1.0; + } + } + } + } +} + +//---------------------------------------------------------------------- +// enforce_dirichlet_dofs_on_face +//---------------------------------------------------------------------- +void enforce_dirichlet_dofs_on_face(ComMod& com_mod, const faceType& lFa, + int dof_start, int num_dofs) +{ + const auto& eq = com_mod.eq[com_mod.cEq]; + const int dof = eq.dof; + const auto& rowPtr = com_mod.rowPtr; + const auto& colPtr = com_mod.colPtr; + auto& R = com_mod.R; + auto& Val = com_mod.Val; + + for (int a = 0; a < lFa.nNo; a++) { + int rowN = lFa.gN(a); + for (int i = dof_start; i < dof_start + num_dofs; i++) { + R(i, rowN) = 0.0; + } + for (int j = rowPtr(rowN); j <= rowPtr(rowN + 1) - 1; j++) { + int colN = colPtr(j); + for (int i = dof_start; i < dof_start + num_dofs; i++) { + for (int k = 0; k < dof; k++) { + Val(i * dof + k, j) = 0.0; + } + } + if (colN == rowN) { + for (int i = dof_start; i < dof_start + num_dofs; i++) { + Val(i * dof + i, j) = 1.0; + } + } + } + } +} + }; diff --git a/Code/Source/solver/set_bc.h b/Code/Source/solver/set_bc.h index b84eadba7..229788d0a 100644 --- a/Code/Source/solver/set_bc.h +++ b/Code/Source/solver/set_bc.h @@ -45,6 +45,13 @@ void set_bc_undef_neu(ComMod& com_mod); void set_bc_undef_neu_l(ComMod& com_mod, const bcType& lBc, const faceType& lFa); +/// @brief Enforce Dirichlet BC at all DOFs of face nodes in assembled system. +void enforce_dirichlet_on_face(ComMod& com_mod, const faceType& lFa, int nsd); + +/// @brief Enforce Dirichlet BC for DOFs [dof_start, dof_start+num_dofs) at face nodes. +void enforce_dirichlet_dofs_on_face(ComMod& com_mod, const faceType& lFa, + int dof_start, int num_dofs); + }; #endif diff --git a/tests/cases/fsi/compare_fsi.py b/tests/cases/fsi/compare_fsi.py new file mode 100644 index 000000000..2a03dafb4 --- /dev/null +++ b/tests/cases/fsi/compare_fsi.py @@ -0,0 +1,215 @@ +""" +Compare monolithic vs partitioned FSI results and generate: +1. Field comparison at final time step +2. Videos of both simulations +3. Coupling performance plot for partitioned FSI +""" + +import numpy as np +import meshio +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from matplotlib.gridspec import GridSpec +from matplotlib.colors import Normalize +import matplotlib.animation as animation +import os, re + +mono_dir = "pipe_3d/1-procs" +part_dir = "pipe_3d_partitioned/1-procs" +coupling_log = "pipe_3d_partitioned/coupling_log.txt" +out_dir = "pipe_3d_partitioned" + +n_steps = 50 + +# ============================================================ +# 1. Parse coupling log +# ============================================================ +print("Parsing coupling log...") +coupling_data = {} # ts -> [(outer, omega, rel_disp), ...] + +with open(coupling_log) as f: + for line in f: + # Format: CP TS-ITER TIME [dB REL_DISP OMEGA] + m = re.match(r'\s*CP\s+(\d+)-(\d+)\s+([\d.e+-]+)\s+\[(-?\d+)\s+([\d.e+-]+)\s+([\d.e+-]+)\]', line) + if m: + ts = int(m.group(1)) + outer = int(m.group(2)) + rel_disp = float(m.group(5)) + omega = float(m.group(6)) + if ts not in coupling_data: + coupling_data[ts] = [] + coupling_data[ts].append((outer, omega, rel_disp)) + +# ============================================================ +# 2. Coupling performance plot +# ============================================================ +print("Creating coupling performance plot...") + +fig, axes = plt.subplots(2, 1, figsize=(10, 7), sharex=True) + +# Top: number of coupling iterations per time step +ts_list = sorted(coupling_data.keys()) +n_iters = [len(coupling_data[ts]) for ts in ts_list] +axes[0].bar(ts_list, n_iters, color='steelblue', width=0.8) +axes[0].set_ylabel('Coupling iterations') +axes[0].set_title('Partitioned FSI Coupling Performance (Aitken relaxation)') +axes[0].set_ylim(0, max(n_iters) + 1) + +# Bottom: convergence history (all time steps overlaid) +cmap = plt.cm.viridis +for i, ts in enumerate(ts_list): + iters = coupling_data[ts] + x = [it[0] for it in iters] + y = [it[2] for it in iters] + color = cmap(i / max(len(ts_list) - 1, 1)) + axes[1].semilogy(x, y, 'o-', color=color, alpha=0.6, lw=1, markersize=4) + +# Highlight first and last +for ts, color, label in [(ts_list[0], 'blue', f'TS {ts_list[0]}'), + (ts_list[-1], 'red', f'TS {ts_list[-1]}')]: + iters = coupling_data[ts] + x = [it[0] for it in iters] + y = [it[2] for it in iters] + axes[1].semilogy(x, y, 'o-', color=color, lw=2, markersize=6, label=label) + +axes[1].axhline(1e-8, color='green', linestyle='--', linewidth=1.5, label='Tolerance (1e-8)') +axes[1].set_xlabel('Coupling iteration within time step') +axes[1].set_ylabel('Relative displacement change') +axes[1].legend(loc='upper right') +axes[1].set_ylim(1e-12, 10) +axes[1].set_xlim(0.5, max(n_iters) + 0.5) + +plt.tight_layout() +plt.savefig(os.path.join(out_dir, 'coupling_performance.png'), dpi=150) +plt.close() +print(f" Saved {out_dir}/coupling_performance.png") + +# ============================================================ +# 3. Final time step comparison +# ============================================================ +print("Comparing final time step fields...") + +mono = meshio.read(os.path.join(mono_dir, f"result_{n_steps:03d}.vtu")) +part = meshio.read(os.path.join(part_dir, f"result_{n_steps:03d}.vtu")) + +print(f"\n{'Field':<20} {'Mono max':>12} {'Part max':>12} {'Rel diff':>12}") +print("-" * 60) +for f in mono.point_data: + m = np.max(np.abs(mono.point_data[f])) + if f in part.point_data: + p = np.max(np.abs(part.point_data[f])) + else: + p = 0.0 + rd = abs(m - p) / (m + 1e-30) + print(f"{f:<20} {m:12.4e} {p:12.4e} {rd:12.4e}") + +# ============================================================ +# 4. Videos (side-by-side velocity magnitude on z=0 slice) +# ============================================================ +print("\nCreating videos...") + +def get_velocity_mag(fname): + """Read VTU and return (points, velocity_magnitude)""" + m = meshio.read(fname) + pts = m.points + vel = m.point_data.get('Velocity', np.zeros((len(pts), 3))) + vmag = np.sqrt(vel[:, 0]**2 + vel[:, 1]**2 + vel[:, 2]**2) + return pts, vmag + +def get_displacement_mag(fname): + """Read VTU and return (points, displacement_magnitude)""" + m = meshio.read(fname) + pts = m.points + # Try FS_Displacement first (monolithic), then Displacement + for key in ['FS_Displacement', 'Displacement']: + if key in m.point_data: + d = m.point_data[key] + if np.max(np.abs(d)) > 1e-15: + return pts, np.sqrt(d[:, 0]**2 + d[:, 1]**2 + d[:, 2]**2) + return pts, np.zeros(len(pts)) + +# Get all data for velocity video +print(" Reading velocity data...") +mono_vel = [] +part_vel = [] +for ts in range(1, n_steps + 1): + mono_pts, mono_vmag = get_velocity_mag(os.path.join(mono_dir, f"result_{ts:03d}.vtu")) + part_pts, part_vmag = get_velocity_mag(os.path.join(part_dir, f"result_{ts:03d}.vtu")) + mono_vel.append((mono_pts, mono_vmag)) + part_vel.append((part_pts, part_vmag)) + +# Use monolithic velocity range for color scaling (partitioned may differ) +vmax_mono = max(np.max(v[1]) for v in mono_vel) + +# Create velocity animation +print(" Rendering velocity video...") +fig, axes = plt.subplots(1, 2, figsize=(14, 5)) + +def plot_frame(frame_idx): + for ax in axes: + ax.clear() + + pts_m, vmag_m = mono_vel[frame_idx] + pts_p, vmag_p = part_vel[frame_idx] + + # Each subplot uses its own normalization for best visibility + sc1 = axes[0].scatter(pts_m[:, 2], pts_m[:, 0], c=vmag_m, s=2, + vmin=0, vmax=vmax_mono, cmap='coolwarm') + axes[0].set_title(f'Monolithic FSI') + axes[0].set_xlabel('z') + axes[0].set_ylabel('x') + axes[0].set_aspect('equal') + + sc2 = axes[1].scatter(pts_p[:, 2], pts_p[:, 0], c=vmag_p, s=2, + vmin=0, vmax=vmax_mono, cmap='coolwarm') + axes[1].set_title(f'Partitioned FSI') + axes[1].set_xlabel('z') + axes[1].set_ylabel('x') + axes[1].set_aspect('equal') + + fig.suptitle(f'Velocity Magnitude (step {frame_idx+1}/{n_steps}, dt=1e-4)', fontsize=14) + return sc1, sc2 + +plot_frame(0) +plt.tight_layout() + +ani = animation.FuncAnimation(fig, plot_frame, frames=n_steps, interval=200, blit=False) +ani.save(os.path.join(out_dir, 'velocity_comparison.gif'), writer='pillow', fps=5) +plt.close() +print(f" Saved {out_dir}/velocity_comparison.gif") + +# Create pressure comparison at final step +print(" Creating final pressure comparison...") +fig, axes = plt.subplots(1, 2, figsize=(14, 5)) + +mono_final = meshio.read(os.path.join(mono_dir, f"result_{n_steps:03d}.vtu")) +part_final = meshio.read(os.path.join(part_dir, f"result_{n_steps:03d}.vtu")) + +p_mono = mono_final.point_data.get('Pressure', np.zeros(len(mono_final.points))) +p_part = part_final.point_data.get('Pressure', np.zeros(len(part_final.points))) +pmax = max(np.max(np.abs(p_mono)), np.max(np.abs(p_part))) + +sc1 = axes[0].scatter(mono_final.points[:, 2], mono_final.points[:, 0], + c=p_mono, s=2, vmin=-pmax, vmax=pmax, cmap='RdBu_r') +axes[0].set_title('Monolithic FSI') +axes[0].set_xlabel('z') +axes[0].set_ylabel('x') +axes[0].set_aspect('equal') +plt.colorbar(sc1, ax=axes[0], label='Pressure') + +sc2 = axes[1].scatter(part_final.points[:, 2], part_final.points[:, 0], + c=p_part, s=2, vmin=-pmax, vmax=pmax, cmap='RdBu_r') +axes[1].set_title('Partitioned FSI') +axes[1].set_xlabel('z') +axes[1].set_ylabel('x') +axes[1].set_aspect('equal') +plt.colorbar(sc2, ax=axes[1], label='Pressure') + +fig.suptitle(f'Pressure at step {n_steps}', fontsize=14) +plt.tight_layout() +plt.savefig(os.path.join(out_dir, 'pressure_comparison.png'), dpi=150) +plt.close() +print(f" Saved {out_dir}/pressure_comparison.png") + +print("\nDone!") diff --git a/tests/cases/fsi/pipe_3d_partitioned/mesh/fluid/mesh-complete.mesh.vtu b/tests/cases/fsi/pipe_3d_partitioned/mesh/fluid/mesh-complete.mesh.vtu new file mode 100644 index 000000000..8ddb1119b --- /dev/null +++ b/tests/cases/fsi/pipe_3d_partitioned/mesh/fluid/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bff21b54094fe528e6f7cd9b57dafc5863b530e112e08f8a62a04431ec14a7a3 +size 67969 diff --git a/tests/cases/fsi/pipe_3d_partitioned/mesh/fluid/mesh-surfaces/end.vtp b/tests/cases/fsi/pipe_3d_partitioned/mesh/fluid/mesh-surfaces/end.vtp new file mode 100644 index 000000000..e7c6c313c --- /dev/null +++ b/tests/cases/fsi/pipe_3d_partitioned/mesh/fluid/mesh-surfaces/end.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:57511367796bab30c07b3b5dbd3b6ec596e433cb1aae2afe763f05ab357655bf +size 13661 diff --git a/tests/cases/fsi/pipe_3d_partitioned/mesh/fluid/mesh-surfaces/interface.vtp b/tests/cases/fsi/pipe_3d_partitioned/mesh/fluid/mesh-surfaces/interface.vtp new file mode 100644 index 000000000..c52d70d2a --- /dev/null +++ b/tests/cases/fsi/pipe_3d_partitioned/mesh/fluid/mesh-surfaces/interface.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4c4a25e331d2c3ccbfa2ffcdd52e0ebbd37a8045989bec9abd7b2c0ed775e2b8 +size 28724 diff --git a/tests/cases/fsi/pipe_3d_partitioned/mesh/fluid/mesh-surfaces/start.vtp b/tests/cases/fsi/pipe_3d_partitioned/mesh/fluid/mesh-surfaces/start.vtp new file mode 100644 index 000000000..55f41e4f7 --- /dev/null +++ b/tests/cases/fsi/pipe_3d_partitioned/mesh/fluid/mesh-surfaces/start.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d90d5fb5ce21d0889b27fd0a592aafa6121d19047bebd801710caf5cba6fb48a +size 13685 diff --git a/tests/cases/fsi/pipe_3d_partitioned/mesh/solid/mesh-complete.mesh.vtu b/tests/cases/fsi/pipe_3d_partitioned/mesh/solid/mesh-complete.mesh.vtu new file mode 100644 index 000000000..08402c8eb --- /dev/null +++ b/tests/cases/fsi/pipe_3d_partitioned/mesh/solid/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2527472b08ece18c95175c0661ca838f90ac8b97406d0c1c77459e94df620fe7 +size 36980 diff --git a/tests/cases/fsi/pipe_3d_partitioned/mesh/solid/mesh-surfaces/end.vtp b/tests/cases/fsi/pipe_3d_partitioned/mesh/solid/mesh-surfaces/end.vtp new file mode 100644 index 000000000..4f2a81290 --- /dev/null +++ b/tests/cases/fsi/pipe_3d_partitioned/mesh/solid/mesh-surfaces/end.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:abb10cc3c374742b481031050fd6ceb2af30b8a14bd5627d89f930a4df01708a +size 12351 diff --git a/tests/cases/fsi/pipe_3d_partitioned/mesh/solid/mesh-surfaces/interface.vtp b/tests/cases/fsi/pipe_3d_partitioned/mesh/solid/mesh-surfaces/interface.vtp new file mode 100644 index 000000000..9a8bafd82 --- /dev/null +++ b/tests/cases/fsi/pipe_3d_partitioned/mesh/solid/mesh-surfaces/interface.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9302a80f392f30ae07b98d0048cb7ad85f2fd3fcb0a86dbe50def69ba5a6cafc +size 28900 diff --git a/tests/cases/fsi/pipe_3d_partitioned/mesh/solid/mesh-surfaces/outside.vtp b/tests/cases/fsi/pipe_3d_partitioned/mesh/solid/mesh-surfaces/outside.vtp new file mode 100644 index 000000000..1c5fbfe2e --- /dev/null +++ b/tests/cases/fsi/pipe_3d_partitioned/mesh/solid/mesh-surfaces/outside.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:29c18dc6c30b33a98016c3df4038433430307f05bc5fd3d2d39469b4a3d4bb63 +size 28813 diff --git a/tests/cases/fsi/pipe_3d_partitioned/mesh/solid/mesh-surfaces/start.vtp b/tests/cases/fsi/pipe_3d_partitioned/mesh/solid/mesh-surfaces/start.vtp new file mode 100644 index 000000000..37fd38c98 --- /dev/null +++ b/tests/cases/fsi/pipe_3d_partitioned/mesh/solid/mesh-surfaces/start.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9c679e42882a1a2db6382133d896879e84468058870863e670fc31212d4e29c2 +size 12335 diff --git a/tests/cases/fsi/pipe_3d_partitioned/result_fluid_001.vtu b/tests/cases/fsi/pipe_3d_partitioned/result_fluid_001.vtu new file mode 100644 index 000000000..126c374d6 --- /dev/null +++ b/tests/cases/fsi/pipe_3d_partitioned/result_fluid_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8cfb0cccdba4fc8c9878764ce5aaaf8552b119e9cf495658112a7eba254ec86c +size 92665 diff --git a/tests/cases/fsi/pipe_3d_partitioned/result_solid_001.vtu b/tests/cases/fsi/pipe_3d_partitioned/result_solid_001.vtu new file mode 100644 index 000000000..3b3d6a90a --- /dev/null +++ b/tests/cases/fsi/pipe_3d_partitioned/result_solid_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d861866a0635b039c78391389901eb4559b29747583b60ab8fb0f699c47d4146 +size 30258 diff --git a/tests/cases/fsi/pipe_3d_partitioned/solver.xml b/tests/cases/fsi/pipe_3d_partitioned/solver.xml new file mode 100644 index 000000000..094518d2d --- /dev/null +++ b/tests/cases/fsi/pipe_3d_partitioned/solver.xml @@ -0,0 +1,199 @@ + + + + + + + 0 + 3 + 1 + 1e-2 + 0.50 + STOP_SIM + true + result + 1 + 1 + 1 + 0 + 1 + 0 + 0 + + + + mesh/fluid/mesh-complete.mesh.vtu + + mesh/fluid/mesh-surfaces/start.vtp + + + mesh/fluid/mesh-surfaces/end.vtp + + + mesh/fluid/mesh-surfaces/interface.vtp + + 0 + + + + mesh/solid/mesh-complete.mesh.vtu + + mesh/solid/mesh-surfaces/start.vtp + + + mesh/solid/mesh-surfaces/end.vtp + + + mesh/solid/mesh-surfaces/interface.vtp + + + mesh/solid/mesh-surfaces/outside.vtp + + 1 + + + + lumen_wall + + + + + false + 1 + 100 + 1e-12 + + + fluid + 1.0 + + 0.04 + + 0.2 + + + + + fsils + + 1e-14 + 200 + 100 + + + + true + true + + + + Neu + 5.0e4 + + + + + + + + + false + 1 + 100 + 1e-12 + + + struct + + M94 + 1.0 + 1.0e7 + 0.3 + + + + + fsils + + 1e-14 + 500 + 100 + + + + true + true + + + + Dir + 0.0 + true + false + (0, 0, 1) + + + + Dir + 0.0 + true + false + (0, 0, 1) + + + + + + + false + 1 + 100 + 1e-12 + 0.3 + + + mesh + 0.0 + 1.0 + 0.3 + + + + + fsils + + 1e-12 + 400 + + + + true + + + + Dir + 0.0 + + + + Dir + 0.0 + + + + + + + 400 + 1e-6 + 0.01 + 1.0 + aitken + lumen_wall + wall_inner + + + diff --git a/tests/conftest.py b/tests/conftest.py index 38f058381..5df8a8dbe 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -108,20 +108,8 @@ def n_proc(request): return request.param -def run_by_name(folder, name, t_max, n_proc=1): - """ - Run a test case and return results - Args: - folder: location from which test will be executed - name: name of svMultiPhysics input file (.xml) - t_max: time step to compare - n_proc: number of processors - - Returns: - Simulation results - """ - - # remove old results folders if they exist +def _run_simulation(folder, name_inp, n_proc): + """Run the solver once, removing any previous output directory first.""" dir_path = os.path.join(folder, str(n_proc) + "-procs") if os.path.exists(dir_path): shutil.rmtree(dir_path) @@ -135,7 +123,7 @@ def run_by_name(folder, name, t_max, n_proc=1): "-np", str(n_proc), exe, - name, + name_inp, ] ) @@ -158,91 +146,48 @@ def run_by_name(folder, name, t_max, n_proc=1): "Exit code {}: {}\n".format(completed.returncode, completed.stderr) ) - # read results - fname = os.path.join( - folder, str(n_proc) + "-procs", "result_" + str(t_max).zfill(3) + ".vtu" - ) + +def _read_result(folder, n_proc, name_result, t_max): + """Read a single VTU result file from the n_proc output directory.""" + result_file = name_result if name_result else "result_" + str(t_max).zfill(3) + ".vtu" + fname = os.path.join(folder, str(n_proc) + "-procs", result_file) if not os.path.exists(fname): raise RuntimeError("No svMultiPhysics output: " + fname) return meshio.read(fname) -def run_with_reference( - base_folder, - test_folder, - fields, - n_proc=1, - t_max=1, - name_ref=None, - name_inp="solver.xml", -): - """ - Run a test case and compare it to a stored reference solution - Args: - folder: location from which test will be executed - fields: array fields to compare (e.g. ["Pressure", "Velocity"]) - n_proc: number of processors - t_max: time step to compare - name_inp: name of svMultiPhysics input file (.xml) - name_ref: name of refence file (.vtu) - """ - # default reference name - if not name_ref: - name_ref = "result_" + str(t_max).zfill(3) + ".vtu" - - # run simulation - folder = os.path.join("cases", base_folder, test_folder) - res = run_by_name(folder, name_inp, t_max, n_proc) - - # read reference - fname = os.path.join(folder, name_ref) - ref = meshio.read(fname) - - # check results +def _compare_fields(res, ref, fields): + """Compare fields between a result and reference mesh. Returns an error string.""" msg = "" for f in fields: - # extract field - if f not in res.point_data.keys(): + if f not in res.point_data: raise ValueError("Field " + f + " not in simulation result") a = res.point_data[f] - if f not in ref.point_data.keys(): + if f not in ref.point_data: raise ValueError("Field " + f + " not in reference result") b = ref.point_data[f] - # truncate last dimension if solution is 2D but reference is 3D if len(a.shape) == 2: if a.shape[1] == 2 and b.shape[1] == 3: assert not np.any(b[:, 2]) b = b[:, :2] - # pick tolerance for current field if f not in RTOL: raise ValueError("No tolerance defined for field " + f) rtol = RTOL[f] - # relative difference (as computed in np.isclose) - # note that we consider rtol as absolute zero (and as relative tolerance) a_fl = a.flatten() b_fl = b.flatten() rel_diff = np.abs(a_fl - b_fl) - rtol - rtol * np.abs(b_fl) - # throw error if not all results are within relative tolerance close = rel_diff <= 0.0 if not np.all(close): - # portion of individual results that are above the tolerance wrong = 1 - np.sum(close) / close.size - - # location of maximum relative difference i_max = rel_diff.argmax() - - # maximum relative difference max_rel = rel_diff[i_max] - - # maximum absolute difference at same location max_abs = np.abs(a_fl[i_max] - b_fl[i_max]) - # throw error message for pytest msg += "Test failed in field " + f + "." msg += " Results differ by more than rtol=" + str(rtol) msg += " in {:.1%}".format(wrong) @@ -250,6 +195,59 @@ def run_with_reference( msg += " Max. rel. difference is" msg += " {:.1e}".format(max_rel) msg += " (abs. {:.1e}".format(max_abs) + ")\n" - # check all fields first and then throw error if any failed + return msg + + +def run_by_name(folder, name, t_max, n_proc=1, name_result=None): + """Run a test case and return results (legacy single-file interface).""" + _run_simulation(folder, name, n_proc) + return _read_result(folder, n_proc, name_result, t_max) + + +def run_with_reference( + base_folder, + test_folder, + fields, + n_proc=1, + t_max=1, + name_ref=None, + name_inp="solver.xml", + name_result=None, + comparisons=None, +): + """ + Run a test case once and compare one or more output files to stored references. + + Args: + fields: fields to compare for the primary output file + n_proc: number of processors + t_max: timestep index used for default file naming + name_inp: solver input XML filename + name_ref: reference VTU filename (default: result_{t_max:03d}.vtu) + name_result: result VTU filename inside {n_proc}-procs/ (default: same as name_ref) + comparisons: list of dicts for multi-file comparison, each with keys: + "fields" — list of field names to compare + "name_ref" — reference VTU filename + "name_result" — result VTU filename (optional, defaults to name_ref) + When provided, fields/name_ref/name_result are ignored. + """ + folder = os.path.join("cases", base_folder, test_folder) + + # Build the comparison list from either the new or legacy parameters + if comparisons is None: + if not name_ref: + name_ref = "result_" + str(t_max).zfill(3) + ".vtu" + comparisons = [{"fields": fields, "name_ref": name_ref, "name_result": name_result}] + + # Run the simulation once + _run_simulation(folder, name_inp, n_proc) + + # Compare each requested output file against its reference + msg = "" + for comp in comparisons: + res = _read_result(folder, n_proc, comp.get("name_result"), t_max) + ref = meshio.read(os.path.join(folder, comp["name_ref"])) + msg += _compare_fields(res, ref, comp["fields"]) + if msg: raise AssertionError(msg) diff --git a/tests/test_fsi.py b/tests/test_fsi.py index ac11fcc04..76d41e85b 100644 --- a/tests/test_fsi.py +++ b/tests/test_fsi.py @@ -33,4 +33,29 @@ def test_pipe_3d_trilinos_ml(n_proc): def test_pipe_RCR_3d(n_proc): test_folder = "pipe_RCR_3d" t_max = 5 - run_with_reference(base_folder, test_folder, fields, n_proc, t_max) \ No newline at end of file + run_with_reference(base_folder, test_folder, fields, n_proc, t_max) + +def test_pipe_3d_partitioned(n_proc): + test_folder = "pipe_3d_partitioned" + t_max = 1 + # A single 1-proc result is the reference for all processor counts, compared + # at the global RTOL. Each sub-mesh (fluid/solid/mesh) is partitioned + # independently, so multi-proc runs only match the 1-proc run once the + # Dirichlet-Neumann coupling has converged tightly. Two settings in + # solver.xml make that possible at the original tolerances: + # - Time_step_size = 1e-2: a larger step weakens the added-mass coupling + # stiffness (which grows ~1/dt for partitioned FSI), so the Aitken + # iteration converges robustly and identically at 1, 3 and 4 procs + # instead of stagnating at 4 procs. + # - Coupling_tolerance = 1e-6: at this dt all proc counts reach the same + # coupling fixed point in ~10 iterations, leaving the partition-count + # round-off 5-10 orders below the global RTOL for every field. + run_with_reference(base_folder, test_folder, fields=[], n_proc=n_proc, t_max=t_max, + comparisons=[ + {"fields": ["Velocity", "Pressure"], + "name_ref": "result_fluid_001.vtu", + "name_result": "result_fluid_001.vtu"}, + {"fields": ["Displacement", "VonMises_stress"], + "name_ref": "result_solid_001.vtu", + "name_result": "result_solid_001.vtu"}, + ]) \ No newline at end of file diff --git a/tests/unitTests/integrator_tests/test_fsi_coupling.cpp b/tests/unitTests/integrator_tests/test_fsi_coupling.cpp new file mode 100644 index 000000000..3c1825ec9 --- /dev/null +++ b/tests/unitTests/integrator_tests/test_fsi_coupling.cpp @@ -0,0 +1,248 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +/** + * @brief Unit tests for FSI interface data exchange helpers. + * + * These tests avoid full mesh I/O so they can check the coupling data + * conventions directly: solid displacement extraction, residual-based fluid + * force extraction semantics, and solid residual force application. + */ + +#include "gtest/gtest.h" + +#include "fsi_coupling.h" + +#include +#include + +// --------------------------------------------------------------------------- +// Helpers +// --------------------------------------------------------------------------- +static faceType make_face(const std::vector& nodes) +{ + faceType face; + face.nNo = static_cast(nodes.size()); + face.gN.resize(face.nNo); + for (int a = 0; a < face.nNo; a++) { + face.gN(a) = nodes[a]; + } + return face; +} + +// =========================================================================== +// Tests +// =========================================================================== + +/// @brief Extract solid displacement at face nodes using the solid equation offset. +TEST(FSICoupling, ExtractSolidDisplacement) +{ + ComMod com_mod; + com_mod.nsd = 3; + + eqType solid_eq; + solid_eq.s = 2; + + auto solid_face = make_face({0, 2, 4}); + + SolutionStates solutions; + auto& Dn = solutions.current.get_displacement(); + Dn.resize(5, 5); + Dn = 0.0; + for (int a = 0; a < solid_face.nNo; a++) { + int Ac = solid_face.gN(a); + for (int i = 0; i < com_mod.nsd; i++) { + Dn(i + solid_eq.s, Ac) = 10.0 * (i + 1) + Ac; + } + } + + auto disp = fsi_coupling::extract_solid_displacement( + com_mod, solid_eq, solid_face, solutions); + + EXPECT_EQ(disp.nrows(), com_mod.nsd); + EXPECT_EQ(disp.ncols(), solid_face.nNo); + for (int a = 0; a < solid_face.nNo; a++) { + int Ac = solid_face.gN(a); + for (int i = 0; i < com_mod.nsd; i++) { + EXPECT_DOUBLE_EQ(disp(i, a), Dn(i + solid_eq.s, Ac)); + } + } +} + +/// @brief Partitioned solid time integration should match the FSI/fluid equation. +TEST(FSICoupling, CopyTimeIntegrationParameters) +{ + eqType source_eq; + source_eq.roInf = 0.5; + source_eq.am = 5.0 / 6.0; + source_eq.af = 2.0 / 3.0; + source_eq.beta = 49.0 / 144.0; + source_eq.gam = 2.0 / 3.0; + source_eq.dof = 4; + source_eq.sym = "FS"; + + eqType target_eq; + target_eq.roInf = 0.5; + target_eq.am = 1.0; + target_eq.af = 2.0 / 3.0; + target_eq.beta = 4.0 / 9.0; + target_eq.gam = 5.0 / 6.0; + target_eq.dof = 3; + target_eq.sym = "ST"; + + fsi_coupling::copy_time_integration_parameters(source_eq, target_eq); + + EXPECT_DOUBLE_EQ(target_eq.roInf, source_eq.roInf); + EXPECT_DOUBLE_EQ(target_eq.am, source_eq.am); + EXPECT_DOUBLE_EQ(target_eq.af, source_eq.af); + EXPECT_DOUBLE_EQ(target_eq.beta, source_eq.beta); + EXPECT_DOUBLE_EQ(target_eq.gam, source_eq.gam); + EXPECT_EQ(target_eq.dof, 3); + EXPECT_EQ(target_eq.sym, "ST"); +} + +/// @brief The residual-recovered force is the opposite of the fluid residual. +TEST(FSICoupling, ExtractFluidResidualForceUsesUnconstrainedRows) +{ + ComMod com_mod; + com_mod.nsd = 3; + com_mod.tnNo = 5; + + auto fluid_face = make_face({1, 3}); + + Array residual(4, com_mod.tnNo); + for (int a = 0; a < residual.ncols(); a++) { + for (int i = 0; i < residual.nrows(); i++) { + residual(i, a) = 100.0 * (i + 1) + a; + } + } + + auto force = fsi_coupling::extract_fluid_residual_force( + com_mod, fluid_face, residual); + + EXPECT_EQ(force.nrows(), com_mod.nsd); + EXPECT_EQ(force.ncols(), fluid_face.nNo); + for (int a = 0; a < fluid_face.nNo; a++) { + int Ac = fluid_face.gN(a); + for (int i = 0; i < com_mod.nsd; i++) { + EXPECT_DOUBLE_EQ(force(i, a), -residual(i, Ac)); + } + } + + // Simulate strong Dirichlet row clearing after extraction. The recovered force + // must retain the unconstrained residual values. + for (int a = 0; a < fluid_face.nNo; a++) { + int Ac = fluid_face.gN(a); + for (int i = 0; i < com_mod.nsd; i++) { + residual(i, Ac) = 0.0; + } + } + for (int a = 0; a < fluid_face.nNo; a++) { + int Ac = fluid_face.gN(a); + for (int i = 0; i < com_mod.nsd; i++) { + EXPECT_NE(force(i, a), residual(i, Ac)); + } + } +} + +/// @brief Residual recovery gives the expected total force for a pressure load. +TEST(FSICoupling, RecoveredFluidForceMatchesKnownPressureLoad) +{ + ComMod com_mod; + com_mod.nsd = 3; + com_mod.tnNo = 4; + + auto fluid_face = make_face({0, 1, 2, 3}); + + const double pressure = 2.5; + const double area = 4.0; + const double nodal_force = pressure * area / fluid_face.nNo; + + Array residual(4, com_mod.tnNo); + residual = 0.0; + for (int a = 0; a < fluid_face.nNo; a++) { + int Ac = fluid_face.gN(a); + residual(0, Ac) = -nodal_force; + residual(3, Ac) = 999.0; // pressure row is not part of the recovered force + } + + auto force = fsi_coupling::extract_fluid_residual_force( + com_mod, fluid_face, residual); + + double total_force[3] = {0.0, 0.0, 0.0}; + for (int a = 0; a < fluid_face.nNo; a++) { + for (int i = 0; i < com_mod.nsd; i++) { + total_force[i] += force(i, a); + } + } + + EXPECT_DOUBLE_EQ(total_force[0], pressure * area); + EXPECT_DOUBLE_EQ(total_force[1], 0.0); + EXPECT_DOUBLE_EQ(total_force[2], 0.0); +} + +/// @brief Applying an external force to the solid subtracts it from R. +TEST(FSICoupling, ApplyForceOnSolidSubtractsExternalForce) +{ + ComMod com_mod; + com_mod.nsd = 3; + com_mod.dof = 4; + com_mod.tnNo = 5; + com_mod.R.resize(com_mod.dof, com_mod.tnNo); + com_mod.R = 0.0; + + eqType solid_eq; + auto solid_face = make_face({0, 2, 4}); + + Array force(com_mod.nsd, solid_face.nNo); + for (int a = 0; a < solid_face.nNo; a++) { + for (int i = 0; i < com_mod.nsd; i++) { + force(i, a) = static_cast((i + 1) * (a + 1)); + } + } + + fsi_coupling::apply_traction_on_solid( + com_mod, solid_eq, solid_face, force); + + for (int a = 0; a < solid_face.nNo; a++) { + int Ac = solid_face.gN(a); + for (int i = 0; i < com_mod.nsd; i++) { + EXPECT_DOUBLE_EQ(com_mod.R(i, Ac), -force(i, a)); + } + } +} + +/// @brief Fluid ALE coordinates are staged from the mesh displacement increment. +TEST(FSICoupling, StagedFluidMeshCoordinatesUseThetaIncrement) +{ + const int nsd = 3; + const int mesh_s = 2; + const double theta = 2.0 / 3.0; + + Array x_ref(nsd, 2); + Array mesh_Do(mesh_s + nsd, 2); + Array mesh_Dn(mesh_s + nsd, 2); + + for (int a = 0; a < x_ref.ncols(); a++) { + for (int i = 0; i < nsd; i++) { + x_ref(i, a) = 10.0 * (i + 1) + a; + } + for (int i = 0; i < mesh_Do.nrows(); i++) { + mesh_Do(i, a) = 0.25 * (i + 1) + a; + mesh_Dn(i, a) = mesh_Do(i, a) + 0.1 * (i + 1) * (a + 1); + } + } + + auto x_stage = fsi_coupling::staged_fluid_mesh_coordinates( + x_ref, mesh_Dn, mesh_Do, mesh_s, nsd, theta); + + EXPECT_EQ(x_stage.nrows(), nsd); + EXPECT_EQ(x_stage.ncols(), x_ref.ncols()); + for (int a = 0; a < x_ref.ncols(); a++) { + for (int i = 0; i < nsd; i++) { + const double increment = mesh_Dn(i + mesh_s, a) + - mesh_Do(i + mesh_s, a); + EXPECT_DOUBLE_EQ(x_stage(i, a), x_ref(i, a) + theta * increment); + } + } +} diff --git a/tests/unitTests/integrator_tests/test_step_equation.cpp b/tests/unitTests/integrator_tests/test_step_equation.cpp new file mode 100644 index 000000000..07472b79e --- /dev/null +++ b/tests/unitTests/integrator_tests/test_step_equation.cpp @@ -0,0 +1,304 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +/** + * @brief Integration tests for Integrator::step_equation() + * + * Tests that step_equation() produces correct results by comparing against + * the monolithic step() method on real solver problems. + * + * These tests require MPI and access to test case data files. + * They are skipped if the test data directory is not available. + */ + +#include "gtest/gtest.h" + +#include "Integrator.h" +#include "Simulation.h" +#include "distribute.h" +#include "initialize.h" +#include "read_files.h" +#include "LinearAlgebra.h" +#include "set_bc.h" + +#include +#include +#include +#include + +// Path to test data, defined via CMake +#ifndef TEST_DATA_DIR +#define TEST_DATA_DIR "" +#endif + +// --------------------------------------------------------------------------- +// MPI environment: initializes MPI once before all tests, finalizes after +// --------------------------------------------------------------------------- +class MPIEnvironment : public ::testing::Environment { +public: + void SetUp() override { + int initialized = 0; + MPI_Initialized(&initialized); + if (!initialized) { + int argc = 0; + char** argv = nullptr; + MPI_Init(&argc, &argv); + } + } + void TearDown() override { + int finalized = 0; + MPI_Finalized(&finalized); + if (!finalized) { + MPI_Finalize(); + } + } +}; + +// Register the MPI environment (runs before any test) +static testing::Environment* const mpi_env = + testing::AddGlobalTestEnvironment(new MPIEnvironment); + +// --------------------------------------------------------------------------- +// Helper: set up a full simulation from an XML file +// --------------------------------------------------------------------------- +static Simulation* setup_simulation(const std::string& xml_path) +{ + // The solver reads mesh files relative to the XML directory, + // so we must chdir there before calling read_files. + std::string dir = xml_path.substr(0, xml_path.find_last_of('/')); + std::string file = xml_path.substr(xml_path.find_last_of('/') + 1); + char orig_dir[4096]; + getcwd(orig_dir, sizeof(orig_dir)); + chdir(dir.c_str()); + + auto simulation = new Simulation(); + read_files_ns::read_files(simulation, file); + distribute(simulation); + Vector init_time(3); + initialize(simulation, init_time); + for (int iEq = 0; iEq < simulation->com_mod.nEq; iEq++) { + add_eq_linear_algebra(simulation->com_mod, simulation->com_mod.eq[iEq]); + } + + chdir(orig_dir); + return simulation; +} + +static void teardown_simulation(Simulation* simulation) +{ + for (int iEq = 0; iEq < simulation->com_mod.nEq; iEq++) { + simulation->com_mod.eq[iEq].linear_algebra->finalize(); + } + delete simulation; +} + +// --------------------------------------------------------------------------- +// Helper: run one time step using step() +// --------------------------------------------------------------------------- +static void run_one_timestep_step(Simulation* simulation) +{ + auto& com_mod = simulation->com_mod; + auto& integrator = simulation->get_integrator(); + auto& solutions = integrator.get_solutions(); + + com_mod.cTS += 1; + com_mod.time += com_mod.dt; + com_mod.cEq = 0; + for (auto& eq : com_mod.eq) { + eq.itr = 0; + eq.ok = false; + } + + integrator.predictor(); + set_bc::set_bc_dir(com_mod, solutions); + integrator.step(); + + // Copy current -> old for next step + solutions.old.get_acceleration() = solutions.current.get_acceleration(); + solutions.old.get_velocity() = solutions.current.get_velocity(); + if (com_mod.dFlag) { + solutions.old.get_displacement() = solutions.current.get_displacement(); + } + com_mod.cplBC.xo = com_mod.cplBC.xn; +} + +// --------------------------------------------------------------------------- +// Helper: run one time step using step_equation() per equation +// --------------------------------------------------------------------------- +static void run_one_timestep_step_equation(Simulation* simulation, + int outer_iters = 1) +{ + auto& com_mod = simulation->com_mod; + auto& integrator = simulation->get_integrator(); + auto& solutions = integrator.get_solutions(); + + com_mod.cTS += 1; + com_mod.time += com_mod.dt; + com_mod.cEq = 0; + for (auto& eq : com_mod.eq) { + eq.itr = 0; + eq.ok = false; + } + + integrator.predictor(); + set_bc::set_bc_dir(com_mod, solutions); + + for (int outer = 0; outer < outer_iters; outer++) { + for (int iEq = 0; iEq < com_mod.nEq; iEq++) { + integrator.step_equation(iEq); + } + } + + // Copy current -> old for next step + solutions.old.get_acceleration() = solutions.current.get_acceleration(); + solutions.old.get_velocity() = solutions.current.get_velocity(); + if (com_mod.dFlag) { + solutions.old.get_displacement() = solutions.current.get_displacement(); + } + com_mod.cplBC.xo = com_mod.cplBC.xn; +} + +// --------------------------------------------------------------------------- +// Helper: compute relative difference between two solution arrays +// --------------------------------------------------------------------------- +static double rel_diff(const Array& a, const Array& b) +{ + double max_diff = 0.0; + double max_val = 0.0; + for (int j = 0; j < a.ncols(); j++) { + for (int i = 0; i < a.nrows(); i++) { + double d = std::abs(a(i,j) - b(i,j)); + if (d > max_diff) max_diff = d; + double v = std::abs(a(i,j)); + if (v > max_val) max_val = v; + } + } + return (max_val > 0) ? max_diff / max_val : max_diff; +} + +// --------------------------------------------------------------------------- +// Helper: check if test data directory exists +// --------------------------------------------------------------------------- +static bool test_data_available() +{ + std::string path = std::string(TEST_DATA_DIR); + if (path.empty()) return false; + struct stat st; + return (stat(path.c_str(), &st) == 0 && S_ISDIR(st.st_mode)); +} + +// =========================================================================== +// Tests +// =========================================================================== + +/// @brief For a single-equation problem, step_equation(0) must produce +/// bit-identical results to step(). +TEST(StepEquation, SingleEquationMatchesStep) +{ + if (!test_data_available()) GTEST_SKIP() << "Test data not available"; + + std::string xml = std::string(TEST_DATA_DIR) + "/fluid/newtonian/solver.xml"; + + // Run with step() + auto sim_a = setup_simulation(xml); + run_one_timestep_step(sim_a); + Array Yn_step = sim_a->get_integrator().get_solutions().current.get_velocity(); + teardown_simulation(sim_a); + + // Run with step_equation(0) + auto sim_b = setup_simulation(xml); + run_one_timestep_step_equation(sim_b, 1); + Array Yn_step_eq = sim_b->get_integrator().get_solutions().current.get_velocity(); + teardown_simulation(sim_b); + + double diff = rel_diff(Yn_step, Yn_step_eq); + EXPECT_LT(diff, 1e-12) << "step_equation(0) should match step() for single-equation problems"; +} + +/// @brief For coupled FSI, sequential step_equation converges to step() +/// result when outer coupling iterations are added. +TEST(StepEquation, CoupledFSIConvergesWithOuterIterations) +{ + if (!test_data_available()) GTEST_SKIP() << "Test data not available"; + + std::string xml = std::string(TEST_DATA_DIR) + "/fsi/pipe_3d/solver.xml"; + + // Run with step() as reference + auto sim_ref = setup_simulation(xml); + run_one_timestep_step(sim_ref); + Array Yn_ref = sim_ref->get_integrator().get_solutions().current.get_velocity(); + Array Dn_ref = sim_ref->get_integrator().get_solutions().current.get_displacement(); + teardown_simulation(sim_ref); + + // Run with 1 outer iteration (single pass) - should have coupling error + auto sim_1 = setup_simulation(xml); + run_one_timestep_step_equation(sim_1, 1); + Array Yn_1 = sim_1->get_integrator().get_solutions().current.get_velocity(); + Array Dn_1 = sim_1->get_integrator().get_solutions().current.get_displacement(); + teardown_simulation(sim_1); + + double diff_vel_1 = rel_diff(Yn_ref, Yn_1); + double diff_disp_1 = rel_diff(Dn_ref, Dn_1); + + // Run with 4 outer iterations - should converge to step() result + auto sim_4 = setup_simulation(xml); + run_one_timestep_step_equation(sim_4, 4); + Array Yn_4 = sim_4->get_integrator().get_solutions().current.get_velocity(); + Array Dn_4 = sim_4->get_integrator().get_solutions().current.get_displacement(); + teardown_simulation(sim_4); + + double diff_vel_4 = rel_diff(Yn_ref, Yn_4); + double diff_disp_4 = rel_diff(Dn_ref, Dn_4); + + // With 4 outer iterations, the coupling error should be orders of magnitude + // smaller than with 1 iteration + EXPECT_LT(diff_vel_4, diff_vel_1 * 1e-4) + << "4 outer iterations should reduce velocity coupling error by >4 orders"; + EXPECT_LT(diff_disp_4, diff_disp_1 * 1e-4) + << "4 outer iterations should reduce displacement coupling error by >4 orders"; + + // With 4 outer iterations, the result should match step() to near machine precision + EXPECT_LT(diff_vel_4, 1e-10) + << "Velocity should match step() after 4 outer iterations"; + EXPECT_LT(diff_disp_4, 1e-10) + << "Displacement should match step() after 4 outer iterations"; +} + +/// @brief The post_assembly callback fires on every Newton iteration. +TEST(StepEquation, PostAssemblyCallbackFires) +{ + if (!test_data_available()) GTEST_SKIP() << "Test data not available"; + + std::string xml = std::string(TEST_DATA_DIR) + "/fluid/newtonian/solver.xml"; + + auto sim = setup_simulation(xml); + auto& com_mod = sim->com_mod; + auto& integrator = sim->get_integrator(); + auto& solutions = integrator.get_solutions(); + + // Set up time step + com_mod.cTS += 1; + com_mod.time += com_mod.dt; + com_mod.cEq = 0; + for (auto& eq : com_mod.eq) { + eq.itr = 0; + eq.ok = false; + } + + integrator.predictor(); + set_bc::set_bc_dir(com_mod, solutions); + + // Count callback invocations + int callback_count = 0; + integrator.step_equation(0, [&callback_count]() { + callback_count++; + }); + + // Callback should fire once per Newton iteration (at least minItr times) + EXPECT_GE(callback_count, com_mod.eq[0].minItr) + << "Callback should fire at least minItr times"; + EXPECT_LE(callback_count, com_mod.eq[0].maxItr) + << "Callback should fire at most maxItr times"; + + teardown_simulation(sim); +}