diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index c546c2822..bb2c50fed 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -204,6 +204,7 @@ set(CSRCS stokes.h stokes.cpp sv_struct.h sv_struct.cpp svZeroD_interface.h svZeroD_interface.cpp + svOneD_interface.h svOneD_interface.cpp txt.h txt.cpp utils.h utils.cpp ustruct.h ustruct.cpp @@ -228,6 +229,7 @@ set(CSRCS SPLIT.c svZeroD_interface/LPNSolverInterface.h svZeroD_interface/LPNSolverInterface.cpp + svOneD_interface/OneDSolverInterface.h svOneD_interface/OneDSolverInterface.cpp BoundaryCondition.h BoundaryCondition.cpp RobinBoundaryCondition.h RobinBoundaryCondition.cpp diff --git a/Code/Source/solver/ComMod.cpp b/Code/Source/solver/ComMod.cpp index 0ac2d4365..93d8d8be7 100644 --- a/Code/Source/solver/ComMod.cpp +++ b/Code/Source/solver/ComMod.cpp @@ -190,10 +190,10 @@ rmshType::rmshType() ///////////////////////////////////////////////////////// -// s v Z e r o D S o l v e r I n t e r f a c e T y p e // +// s v Z e r o D S o l v e r I n t e r f a c e D a t a // ///////////////////////////////////////////////////////// -void svZeroDSolverInterfaceType::set_data(const svZeroDSolverInterfaceParameters& params) +void svZeroDSolverInterfaceData::set_data(const svZeroDSolverInterfaceParameters& params) { if (!params.defined()) { return; @@ -217,3 +217,13 @@ void svZeroDSolverInterfaceType::set_data(const svZeroDSolverInterfaceParameters has_data = true; } + +void svOneDSolverInterfaceData::set_data(const svOneDSolverInterfaceParameters& params) +{ + if (!params.defined()) { + return; + } + + solver_library = params.shared_library(); + has_data = true; +} diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 86db290f4..1a3a43893 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -204,6 +204,10 @@ class bcType // Coupled BC class CoupledBoundaryCondition coupled_bc; + + // svOneD: per-face 1D solver input file path (set when Time_dependence=Coupled + // and svOneDSolver_interface is active). + std::string oned_input_file; }; /// @brief Class storing data for B-Splines. @@ -771,15 +775,21 @@ class cplFaceType // RCR type BC rcrType RCR; + + // svOneD: path to the per-face 1D solver input file. + std::string oned_input_file; + + // Whether this face uses RCR (Windkessel) boundary condition. + bool isRCR = false; }; //---------------------------- -// svZeroDSolverInterfaceType +// svZeroDSolverInterfaceData //---------------------------- // This class stores information used to interface to // the svZeroDSolver. // -class svZeroDSolverInterfaceType +class svZeroDSolverInterfaceData { public: @@ -807,6 +817,32 @@ class svZeroDSolverInterfaceType void set_data(const svZeroDSolverInterfaceParameters& params); }; +/// \brief Stores information used to interface with svOneDSolver. +/// +/// This class stores the global svOneDSolver interface settings read from the +/// solver XML file. Per-face 1D input files are stored separately in +/// cplFaceType::oned_input_file. +class svOneDSolverInterfaceData +{ + public: + /// \brief Path to the svOneDSolver interface shared library. + /// + /// This may be provided either with the platform-specific extension + /// (.so/.dylib) or without it. + std::string solver_library; + + /// \brief True if the svOneDSolver interface settings were read from XML. + /// + /// This is set to true after the svOneDSolver_interface XML element has + /// been parsed successfully. + bool has_data = false; + + /// \brief Read svOneDSolver interface settings from parsed parameters. + /// + /// \param params Parsed svOneDSolver interface parameters. + void set_data(const svOneDSolverInterfaceParameters& params); +}; + /// @brief For coupled 0D-3D problems // class cplBCType @@ -822,6 +858,9 @@ class cplBCType // Whether to use svZeroD bool useSvZeroD = false; + // Whether to use svOneD (svOneDSolver) + bool useSvOneD = false; + // Whether to initialize RCR from flow data bool initRCR = false; @@ -851,7 +890,11 @@ class cplBCType std::string commuName; //std::string commuName = ".CPLBC_0D_3D.tmp"; - svZeroDSolverInterfaceType svzerod_solver_interface; + /// @brief Data structure used for coupling with svZeroD code + svZeroDSolverInterfaceData svzerod_solver_interface; + + /// @brief Data structure used for coupling with svOneD code + svOneDSolverInterfaceData svOneD_solver_interface; /// @brief The name of history file containing "X" std::string saveName; diff --git a/Code/Source/solver/CoupledBoundaryCondition.cpp b/Code/Source/solver/CoupledBoundaryCondition.cpp index 5f7602e3e..2d2ed353d 100644 --- a/Code/Source/solver/CoupledBoundaryCondition.cpp +++ b/Code/Source/solver/CoupledBoundaryCondition.cpp @@ -22,6 +22,17 @@ CoupledBoundaryCondition::CoupledBoundaryCondition(const CoupledBoundaryConditio , bc_type_(other.bc_type_) , block_name_(other.block_name_) , face_name_(other.face_name_) + , oned_input_file_(other.oned_input_file_) + , oned_ramp_steps_(other.oned_ramp_steps_) + , oned_ramp_ref_pressure_(other.oned_ramp_ref_pressure_) + , oned_relax_factor_(other.oned_relax_factor_) + , ramp_step_count_(other.ramp_step_count_) + , P_prev_sent_old_(other.P_prev_sent_old_) + , P_prev_sent_new_(other.P_prev_sent_new_) + , Q_prev_sent_(other.Q_prev_sent_) + , P_neu_prev_(other.P_neu_prev_) + , Q_input_prev_old_(other.Q_input_prev_old_) + , Q_input_prev_new_(other.Q_input_prev_new_) , Qo_(other.Qo_) , Qn_(other.Qn_) , Po_(other.Po_) @@ -53,6 +64,17 @@ CoupledBoundaryCondition& CoupledBoundaryCondition::operator=(const CoupledBound bc_type_ = other.bc_type_; block_name_ = other.block_name_; face_name_ = other.face_name_; + oned_input_file_ = other.oned_input_file_; + oned_ramp_steps_ = other.oned_ramp_steps_; + oned_ramp_ref_pressure_ = other.oned_ramp_ref_pressure_; + oned_relax_factor_ = other.oned_relax_factor_; + ramp_step_count_ = other.ramp_step_count_; + P_prev_sent_old_ = other.P_prev_sent_old_; + P_prev_sent_new_ = other.P_prev_sent_new_; + Q_prev_sent_ = other.Q_prev_sent_; + P_neu_prev_ = other.P_neu_prev_; + Q_input_prev_old_ = other.Q_input_prev_old_; + Q_input_prev_new_ = other.Q_input_prev_new_; Qo_ = other.Qo_; Qn_ = other.Qn_; Po_ = other.Po_; @@ -83,6 +105,17 @@ CoupledBoundaryCondition::CoupledBoundaryCondition(CoupledBoundaryCondition&& ot , bc_type_(other.bc_type_) , block_name_(std::move(other.block_name_)) , face_name_(std::move(other.face_name_)) + , oned_input_file_(std::move(other.oned_input_file_)) + , oned_ramp_steps_(other.oned_ramp_steps_) + , oned_ramp_ref_pressure_(other.oned_ramp_ref_pressure_) + , oned_relax_factor_(other.oned_relax_factor_) + , ramp_step_count_(other.ramp_step_count_) + , P_prev_sent_old_(other.P_prev_sent_old_) + , P_prev_sent_new_(other.P_prev_sent_new_) + , Q_prev_sent_(other.Q_prev_sent_) + , P_neu_prev_(other.P_neu_prev_) + , Q_input_prev_old_(other.Q_input_prev_old_) + , Q_input_prev_new_(other.Q_input_prev_new_) , Qo_(other.Qo_) , Qn_(other.Qn_) , Po_(other.Po_) @@ -129,6 +162,17 @@ CoupledBoundaryCondition& CoupledBoundaryCondition::operator=(CoupledBoundaryCon bc_type_ = other.bc_type_; block_name_ = std::move(other.block_name_); face_name_ = std::move(other.face_name_); + oned_input_file_ = std::move(other.oned_input_file_); + oned_ramp_steps_ = other.oned_ramp_steps_; + oned_ramp_ref_pressure_ = other.oned_ramp_ref_pressure_; + oned_relax_factor_ = other.oned_relax_factor_; + ramp_step_count_ = other.ramp_step_count_; + P_prev_sent_old_ = other.P_prev_sent_old_; + P_prev_sent_new_ = other.P_prev_sent_new_; + Q_prev_sent_ = other.Q_prev_sent_; + P_neu_prev_ = other.P_neu_prev_; + Q_input_prev_old_ = other.Q_input_prev_old_; + Q_input_prev_new_ = other.Q_input_prev_new_; Qo_ = other.Qo_; Qn_ = other.Qn_; Po_ = other.Po_; @@ -208,6 +252,16 @@ const std::string& CoupledBoundaryCondition::get_block_name() const return block_name_; } +const std::string& CoupledBoundaryCondition::get_oned_input_file() const +{ + return oned_input_file_; +} + +void CoupledBoundaryCondition::set_oned_input_file(const std::string& path) +{ + oned_input_file_ = path; +} + void CoupledBoundaryCondition::set_solution_ids(int flow_id, int pressure_id, double in_out_sign) { flow_sol_id_ = flow_id; @@ -401,6 +455,14 @@ void CoupledBoundaryCondition::distribute(const ComMod& com_mod, const CmMod& cm // Distribute block name cm.bcast(cm_mod, block_name_); + // Distribute 1D input file path + cm.bcast(cm_mod, oned_input_file_); + + // Distribute 1D ramp and relaxation parameters + cm.bcast(cm_mod, &oned_ramp_steps_); + cm.bcast(cm_mod, &oned_ramp_ref_pressure_); + cm.bcast(cm_mod, &oned_relax_factor_); + // Distribute face name cm.bcast(cm_mod, face_name_); @@ -986,7 +1048,7 @@ void CappingSurface::init_cap_face_quadrature(const ComMod& com_mod) try { if (nsd != cap_nsd_) { - throw CappingSurfaceBaseException("[CappingSurface::init_cap_face_quadrature] Cap surface requires nsd=3."); + throw CappingSurfaceGeometryException("[CappingSurface::init_cap_face_quadrature] Cap surface requires nsd=3."); } face_->nG = 1; @@ -1084,7 +1146,7 @@ std::pair> CappingSurface::compute_jacobian_and_normal(co Jac = utils::norm(n); if (utils::is_zero(Jac)) { - throw CappingSurfaceBaseException("[CappingSurface::compute_jacobian_and_normal] Zero Jacobian at Gauss point " + + throw CappingSurfaceGeometryException("[CappingSurface::compute_jacobian_and_normal] Zero Jacobian at Gauss point " + std::to_string(g)); } @@ -1311,3 +1373,23 @@ void CoupledBoundaryCondition::bcast_coupled_neumann_pressure(const CmMod& cm_mo cm.bcast(cm_mod, &pr); set_pressure(pr); } + +void CoupledBoundaryCondition::bcast_coupled_dir_flowrate(const CmMod& cm_mod, cmType& cm) +{ + if (cm.seq()) { + return; + } + using namespace consts; + if (get_bc_type() != BoundaryConditionType::bType_Dir) { + return; + } + double Qo = 0.0; + double Qn = 0.0; + if (cm.mas(cm_mod)) { + Qo = get_Qo(); + Qn = get_Qn(); + } + cm.bcast(cm_mod, &Qo); + cm.bcast(cm_mod, &Qn); + set_flowrates(Qo, Qn); +} diff --git a/Code/Source/solver/CoupledBoundaryCondition.h b/Code/Source/solver/CoupledBoundaryCondition.h index ae6d54a6d..858354a3f 100644 --- a/Code/Source/solver/CoupledBoundaryCondition.h +++ b/Code/Source/solver/CoupledBoundaryCondition.h @@ -98,6 +98,12 @@ class CappingSurfaceInvalidElementNodesException : public CappingSurfaceBaseExce std::to_string(eNoN) + " (expected " + std::to_string(expected) + ").") {} }; +/// @brief Geometry or Jacobian error during cap surface integration. +class CappingSurfaceGeometryException : public CappingSurfaceBaseException { +public: + explicit CappingSurfaceGeometryException(const std::string& detail) : CappingSurfaceBaseException(detail) {} +}; + /// @brief Cap face quadrature (shape functions on TRI3) setup failed. class CappingSurfaceQuadratureException : public CappingSurfaceBaseException { public: @@ -106,6 +112,18 @@ class CappingSurfaceQuadratureException : public CappingSurfaceBaseException { "functions: " + nested) {} }; +/// @brief Inconsistent cap connectivity or assembly indexing. +class CappingSurfaceAssemblyException : public CappingSurfaceBaseException { +public: + explicit CappingSurfaceAssemblyException(const std::string& detail) : CappingSurfaceBaseException(detail) {} +}; + +/// @brief Failure copying a \ref CappingSurface or its internal face. +class CappingSurfaceCopyException : public CappingSurfaceBaseException { +public: + explicit CappingSurfaceCopyException(std::string msg) : CappingSurfaceBaseException(std::move(msg)) {} +}; + /// @brief Capping surface geometry and integration for a coupled boundary. class CappingSurface { public: @@ -133,12 +151,11 @@ class CappingSurface { /// Surface velocity flux through the cap using \a st columns indexed by cap IEN / GlobalNodeID (master / serial). double integrate_velocity_flux(const CapGlobalMeshState& st, bool use_Yn_velocity, consts::MechanicalConfigurationType cfg); - + /// @brief Compute the cap contribution to the linear solver face (fills \ref valM_; safe under \c const *this). void compute_valM(consts::MechanicalConfigurationType cfg, const CapGlobalMeshState& st) const; /// @brief Get the cap face. - faceType* face() { return face_.get(); } const faceType* face() const { return face_.get(); } /// @brief Get the cap contribution. @@ -191,7 +208,27 @@ class CoupledBoundaryCondition { /// @brief svZeroD coupling data std::string block_name_; ///< Block name in svZeroDSolver configuration std::string face_name_; ///< Face name from the mesh - + /// @brief svOneD coupling data + std::string oned_input_file_; ///< Path to svOneDSolver input file (empty for svZeroD BCs) + /// @brief Pressure ramp parameters for 1D coupling initialization. + int oned_ramp_steps_ = 0; ///< Number of ramp steps (0 = disabled) + double oned_ramp_ref_pressure_ = 0.0; ///< Reference pressure at step used for ramping pressure + + /// @brief Under-relaxation factor for pressure or flowrate passed to the 1D solver. + /// Applied as: P_sent = omega * P_new + (1 - omega) * P_prev_sent. + /// Range: (0, 1]. Default 1.0 = no relaxation. + double oned_relax_factor_ = 1.0; + + /// @brief Ramp/relax runtime state (shared by both 0D and 1D coupling). + /// Updated only on committed ('L') time steps. + int ramp_step_count_ = 0; ///< Number of committed steps taken (used for ramp fraction). + double P_prev_sent_old_ = 0.0; ///< Under-relaxed pressure sent at t_old on last 'L' step (DIR input). + double P_prev_sent_new_ = 0.0; ///< Under-relaxed pressure sent at t_new on last 'L' step (DIR input). + double Q_prev_sent_ = 0.0; ///< Under-relaxed flow output on last 'L' step (DIR output). + double P_neu_prev_ = 0.0; ///< Under-relaxed pressure output on last 'L' step (NEU output). + double Q_input_prev_old_ = 0.0; ///< Under-relaxed flow sent at t_old on last 'L' step (NEU input). + double Q_input_prev_new_ = 0.0; ///< Under-relaxed flow sent at t_new on last 'L' step (NEU input). + /// @brief Flowrate data double Qo_ = 0.0; ///< Flowrate at old timestep (t_n) double Qn_ = 0.0; ///< Flowrate at new timestep (t_{n+1}) @@ -289,6 +326,82 @@ class CoupledBoundaryCondition { /// @return Block name const std::string& get_block_name() const; + /// @brief Get the svOneD input file path + /// @return Path to the 1D solver input file (empty for svZeroD BCs) + const std::string& get_oned_input_file() const; + + /// @brief Set the svOneD input file path + void set_oned_input_file(const std::string& path); + + bool is_svOneD_face() const { return !oned_input_file_.empty(); } + + /// @brief Get the pressure ramp step count (0 = disabled). + int get_oned_ramp_steps() const { return oned_ramp_steps_; } + + /// @brief Set the pressure ramp parameters for 1D coupling initialization. + /// @param steps Number of time steps over which to ramp (0 = disabled). + /// @param P_ref Reference pressure at step 0 (typically the 1D initial pressure). + void set_oned_ramp(int steps, double P_ref) { + oned_ramp_steps_ = steps; + oned_ramp_ref_pressure_ = P_ref; + } + + /// @brief Get the ramp reference pressure. + double get_oned_ramp_ref_pressure() const { return oned_ramp_ref_pressure_; } + + /// @brief Get the under-relaxation factor for DIR coupling pressure (1.0 = no relaxation). + double get_oned_relax_factor() const { return oned_relax_factor_; } + + /// @brief Set the under-relaxation factor for DIR coupling pressure. + /// @param omega Relaxation factor in (0, 1]. 1.0 disables relaxation. + void set_oned_relax_factor(double omega) { oned_relax_factor_ = omega; } + + // ========================================================================= + // Ramp/relax runtime state accessors (shared by 0D and 1D coupling) + // ========================================================================= + + /// @brief Get the number of committed time steps (used for ramp fraction). + int get_ramp_step_count() const { return ramp_step_count_; } + + /// @brief Increment the committed step counter (call once per 'L' step). + void increment_ramp_step_count() { ramp_step_count_++; } + + /// @brief Get the under-relaxed pressure sent at t_old on the last 'L' step (DIR input). + double get_P_prev_sent_old() const { return P_prev_sent_old_; } + + /// @brief Get the under-relaxed pressure sent at t_new on the last 'L' step (DIR input). + double get_P_prev_sent_new() const { return P_prev_sent_new_; } + + /// @brief Set the DIR input pressure history (call on 'L' steps). + void set_P_prev_sent(double old_val, double new_val) { + P_prev_sent_old_ = old_val; + P_prev_sent_new_ = new_val; + } + + /// @brief Get the under-relaxed flow output on the last 'L' step (DIR output). + double get_Q_prev_sent() const { return Q_prev_sent_; } + + /// @brief Set the DIR output flow history (call on 'L' steps). + void set_Q_prev_sent(double Q) { Q_prev_sent_ = Q; } + + /// @brief Get the under-relaxed pressure output on the last 'L' step (NEU output). + double get_P_neu_prev() const { return P_neu_prev_; } + + /// @brief Set the NEU output pressure history (call on 'L' steps). + void set_P_neu_prev(double P) { P_neu_prev_ = P; } + + /// @brief Get the under-relaxed flow sent at t_old on the last 'L' step (NEU input). + double get_Q_input_prev_old() const { return Q_input_prev_old_; } + + /// @brief Get the under-relaxed flow sent at t_new on the last 'L' step (NEU input). + double get_Q_input_prev_new() const { return Q_input_prev_new_; } + + /// @brief Set the NEU input flow history (call on 'L' steps). + void set_Q_input_prev(double old_val, double new_val) { + Q_input_prev_old_ = old_val; + Q_input_prev_new_ = new_val; + } + /// @brief Set the svZeroD solution IDs for flow and pressure /// @param flow_id Flow solution ID /// @param pressure_id Pressure solution ID @@ -414,6 +527,9 @@ class CoupledBoundaryCondition { /// @brief Master reads Neumann pressure, one scalar \c MPI_Bcast, all ranks set pressure (svZeroD sync). void bcast_coupled_neumann_pressure(const CmMod& cm_mod, cmType& cm); + /// @brief Master reads DIR flowrates (Qo, Qn), two scalar \c MPI_Bcast, all ranks set flowrates (svZeroD sync). + void bcast_coupled_dir_flowrate(const CmMod& cm_mod, cmType& cm); + }; #endif // COUPLED_BOUNDARY_CONDITION_H diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index e20dcf2f7..194dbc18d 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -466,8 +466,12 @@ void BoundaryConditionRCRParameters::set_values(tinyxml2::XMLElement* xml_elem) CouplingInterfaceParameters::CouplingInterfaceParameters() { bool required = false; - set_parameter("svZeroDSolver_block", "", !required, svzerod_solver_block); - set_parameter("Chamber_cap_surface", "", !required, chamber_cap_surface); + set_parameter("svZeroDSolver_block", "", !required, svzerod_solver_block); + set_parameter("Chamber_cap_surface", "", !required, chamber_cap_surface); + set_parameter("svOneDSolver_input_file","", !required, svoned_input_file); + set_parameter("Ramp_steps", 0, !required, coupling_ramp_steps); + set_parameter("Ramp_ref_pressure", 0.0, !required, coupling_ramp_ref_pressure); + set_parameter("Relax_factor", 1.0, !required, coupling_relax_factor); } void CouplingInterfaceParameters::set_values(tinyxml2::XMLElement* xml_elem) @@ -1175,6 +1179,37 @@ void svZeroDSolverInterfaceParameters::set_values(tinyxml2::XMLElement* xml_elem value_set = true; } +////////////////////////////////////////////////////////// +// svOneDSolverInterfaceParameters // +////////////////////////////////////////////////////////// + +const std::string svOneDSolverInterfaceParameters::xml_element_name_ = "svOneDSolver_interface"; + +svOneDSolverInterfaceParameters::svOneDSolverInterfaceParameters() +{ + bool required = true; + + set_parameter("Coupling_type", "", required, coupling_type); + set_parameter("Shared_library", "", required, shared_library); + // Unlike svZeroDSolver (which uses a single global Input_file for all faces), + // svOneDSolver requires each coupled face to supply its own input file. + // Per-face input files are specified via + // inside each element. +} + +void svOneDSolverInterfaceParameters::set_values(tinyxml2::XMLElement* xml_elem) +{ + std::string error_msg = "Unknown svOneDSolver_interface XML element '"; + + using std::placeholders::_1; + using std::placeholders::_2; + std::function ftpr = + std::bind(&svOneDSolverInterfaceParameters::set_parameter_value, *this, _1, _2); + xml_util_set_parameters(ftpr, xml_elem, error_msg); + + value_set = true; +} + ////////////////////////////////////////////////////////// // OutputParameters // ////////////////////////////////////////////////////////// @@ -1823,6 +1858,7 @@ void DomainParameters::print_parameters() fluid_viscosity.print_parameters(); solid_viscosity.print_parameters(); + } //------------ @@ -2415,6 +2451,9 @@ void EquationParameters::set_values(tinyxml2::XMLElement* eq_elem, DomainParamet } else if (name == svZeroDSolverInterfaceParameters::xml_element_name_) { svzerodsolver_interface_parameters.set_values(item); + } else if (name == svOneDSolverInterfaceParameters::xml_element_name_) { + svonedsolver_interface_parameters.set_values(item); + } else if (name == DomainParameters::xml_element_name_) { auto domain_params = new DomainParameters(); domain_params->set_values(item); diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index de383e80f..29811616f 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -679,6 +679,47 @@ class svZeroDSolverInterfaceParameters : public ParameterLists bool value_set = false; }; +//---------------------------------- +// svOneDSolverInterfaceParameters +//---------------------------------- +/// @brief Parameters for coupling to the svOneDSolver (1D blood-flow solver). +/// +/// XML element: \code {.xml} +/// +/// Implicit +/// /path/to/libsvoned_interface +/// +/// \endcode +/// +/// Notes +/// ----- +/// Coupling_type: "Explicit" | "Implicit" | "Semi-implicit" +/// Controls how the 3D Newton iteration couples to the 1D solver +/// (same semantics as the 0D coupling_type). +/// Shared_library: Path to the 1D interface shared library. The +/// extension (.so or .dylib) may be omitted; it will be appended +/// automatically based on the platform. +/// +/// Each coupled face specifies its own input file via +/// ... +/// inside the corresponding element. +// +class svOneDSolverInterfaceParameters : public ParameterLists +{ + public: + svOneDSolverInterfaceParameters(); + + static const std::string xml_element_name_; + + bool defined() const { return value_set; }; + void set_values(tinyxml2::XMLElement* xml_elem); + + Parameter coupling_type; + Parameter shared_library; + + bool value_set = false; +}; + /// @brief Body force over a mesh using the "Add_BF" command. /// /// \code {.xml} @@ -745,10 +786,12 @@ class BoundaryConditionRCRParameters : public ParameterLists /// @brief svZeroDSolver coupling options under Add_BC (with Time_dependence Coupled). /// /// \code {.xml} -/// -/// LV_IN -/// mesh/mesh-surfaces/endo_cap.vtp -/// +// +// OneDfilename.in +// 100 +// 0.0 +// 0.3 +// /// \endcode class CouplingInterfaceParameters : public ParameterLists { @@ -763,6 +806,24 @@ class CouplingInterfaceParameters : public ParameterLists Parameter svzerod_solver_block; Parameter chamber_cap_surface; + // Path to the svOneDSolver .in input file for this face (1D coupling). + Parameter svoned_input_file; + + // Ramp for 1D coupling initialization (both DIR and NEU coupling). + // Over the first Coupling_ramp_steps committed time steps the value passed + // to the 1D solver is linearly ramped: + // DIR: pressure P ramped from Coupling_ramp_ref_pressure to actual 3D P. + // NEU: output pressure P is ramped before being applied to the 3D domain. + // Set Coupling_ramp_steps = 0 (default) to disable. + Parameter coupling_ramp_steps; + Parameter coupling_ramp_ref_pressure; + + // Under-relaxation factor for the value passed to the 1D solver (both DIR and NEU coupling). + // DIR: P_sent = omega * P_target + (1 - omega) * P_prev_sent + // NEU: Q_sent = omega * Q_target + (1 - omega) * Q_prev_sent + // Range: (0, 1]. Default 1.0 = no relaxation. + Parameter coupling_relax_factor; + bool value_set = false; }; @@ -1528,6 +1589,8 @@ class EquationParameters : public ParameterLists svZeroDSolverInterfaceParameters svzerodsolver_interface_parameters; + svOneDSolverInterfaceParameters svonedsolver_interface_parameters; + DomainParameters* default_domain = nullptr; std::vector domains; diff --git a/Code/Source/solver/baf_ini.cpp b/Code/Source/solver/baf_ini.cpp index 2b3b1b071..63bf7aaa8 100644 --- a/Code/Source/solver/baf_ini.cpp +++ b/Code/Source/solver/baf_ini.cpp @@ -12,6 +12,7 @@ #include "set_bc.h" #include "utils.h" #include "svZeroD_interface.h" +#include "svOneD_interface.h" #include "fsils_api.hpp" #include "fils_struct.hpp" @@ -122,6 +123,10 @@ void baf_ini(Simulation* simulation, SolutionStates& solutions) com_mod.cplBC.fa[i].name = com_mod.msh[iM].fa[iFa].name; com_mod.cplBC.fa[i].y = 0.0; + if (!bc.oned_input_file.empty()) { + com_mod.cplBC.fa[i].oned_input_file = bc.oned_input_file; + } + if (utils::btest(bc.bType, iBC_Dir)) { com_mod.cplBC.fa[i].bGrp = CplBCType::cplBC_Dir; @@ -139,6 +144,9 @@ void baf_ini(Simulation* simulation, SolutionStates& solutions) com_mod.cplBC.fa[i].RCR.Rd = bc.RCR.Rd; com_mod.cplBC.fa[i].RCR.Pd = bc.RCR.Pd; com_mod.cplBC.fa[i].RCR.Xo = bc.RCR.Xo; + if (utils::btest(bc.bType, iBC_RCR)) { + com_mod.cplBC.fa[i].isRCR = true; + } } else { throw std::runtime_error("Not a compatible cplBC_type"); } @@ -162,6 +170,10 @@ void baf_ini(Simulation* simulation, SolutionStates& solutions) svZeroD::init_svZeroD(com_mod, cm_mod); } + if (com_mod.cplBC.useSvOneD) { + svOneD::init_svOneD(com_mod, cm_mod); + } + // Initialize cap integration for Coupled boundary conditions for (int iEq = 0; iEq < com_mod.nEq; iEq++) { auto& eq = com_mod.eq[iEq]; @@ -763,10 +775,21 @@ void fsi_ls_ini(ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, const faceTyp } fsils_bc_create(com_mod.lhs, lsPtr, lFa.nNo, nsd, BcType::BC_TYPE_Dir, gNodes, sVl); } + + } else if (btest(lBc.bType, iBC_Neu) || btest(lBc.bType, iBC_Coupled)) { + // For Coupled-DIR BCs: iBC_Dir was cleared in read_files but the face DOFs must still + // be excluded from the linear solve (Ax=b). Register as BC_TYPE_Dir so the + // preconditioner zeros out those rows/columns, preventing the solver from + // overwriting the velocity values set by set_bc_dir. + if (btest(lBc.bType, iBC_Coupled) && + lBc.coupled_bc.get_bc_type() == consts::BoundaryConditionType::bType_Dir) { + lsPtr = lsPtr + 1; + lBc.lsPtr = lsPtr; + sVl = 0.0; + fsils_bc_create(com_mod.lhs, lsPtr, lFa.nNo, nsd, BcType::BC_TYPE_Dir, gNodes, sVl); - } else if (btest(lBc.bType, iBC_Neu)) { - // Compute integral of normal vector over the face (needed for resistance BC/0D-coupling) - if (btest(lBc.bType, iBC_res)) { + } else if (btest(lBc.bType, iBC_res)) { + // Compute integral of normal vector over the face (needed for resistance BC/0D-coupling) sV = 0.0; for (int e = 0; e < lFa.nEl; e++) { if (lFa.eType == ElementType::NRB) { diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 64acce552..a2c7bd1c3 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -550,24 +550,30 @@ void distribute(Simulation* simulation) cm.bcast_enum(cm_mod, &cplBC.schm); cm.bcast(cm_mod, &cplBC.useGenBC); cm.bcast(cm_mod, &cplBC.useSvZeroD); + cm.bcast(cm_mod, &cplBC.useSvOneD); if (cplBC.useGenBC) { if (cm.slv(cm_mod)) { cplBC.nX = 0; cplBC.xo.resize(cplBC.nX); } + } - } else if (cplBC.useSvZeroD) { - if (cm.slv(cm_mod)) { - cplBC.nX = 0; - cplBC.xo.resize(cplBC.nX); - } + if (cplBC.useSvOneD) { + // Broadcast the svOneD solver interface data so that ALL ranks can call + // init_svOneD / calc_svOneD (which use MPI_Bcast collectives that require + // every rank to participate). Without this, slave processes have + // has_data = false and throw immediately inside init_svOneD. + cm.bcast(cm_mod, &cplBC.svOneD_solver_interface.has_data); + cm.bcast(cm_mod, cplBC.svOneD_solver_interface.solver_library); + } - } else { - // RCR (Windkessel): nX/xo sized in read_files from nFa; not genBC/svZeroD. + if (!cplBC.useGenBC) { + // Broadcast nX and xo: when RCR faces coexist with svZeroD/svOneD, nX > 0 + // and xo must be distributed to all slave processes so rcr_init works correctly. cm.bcast(cm_mod, &cplBC.nX); if (cplBC.xo.size() == 0) { - cplBC.xo.resize(cplBC.nX); + cplBC.xo.resize(cplBC.nX); } if (cplBC.nX != 0) { cm.bcast(cm_mod, cplBC.xo); diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 6c6723bce..fa542a90d 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -25,6 +25,7 @@ #include #include #include +#include "Core/Exception.h" namespace read_files_ns { @@ -128,12 +129,15 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, { using namespace consts; auto bc_type = bc_params->type.value(); + BoundaryConditionType coupled_bc_type = BoundaryConditionType::bType_Neu; if (std::set{"Dirichlet", "Dir"}.count(bc_type)) { - lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_Dir)); + lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_Dir)); + coupled_bc_type = BoundaryConditionType::bType_Dir; } else if (std::set{"Neumann", "Neu"}.count(bc_type)) { lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_Neu)); + coupled_bc_type = BoundaryConditionType::bType_Neu; if ((lEq.phys == EquationType::phys_fluid) || (lEq.phys == EquationType::phys_FSI)) { lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_bfs)); } @@ -245,38 +249,117 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, read_fourier_coeff_values_file(file_name, lBc); } - // Coupling to a 0D model: + // Coupling to a 0D/1D model: // - GenBC / cplBC: Time_dependence Coupled without (cplBC.fa). // - svZeroDSolver: Time_dependence Coupled + (CoupledBoundaryCondition). + // - svOneDSolver: Time_dependence Coupled + (CoupledBoundaryCondition). // } else if (ctmp == "Coupled") { auto& face_name = com_mod.msh[lBc.iM].fa[lBc.iFa].name; const bool svzd_iface = com_mod.cplBC.svzerod_solver_interface.has_data; + const bool svOneD_iface = com_mod.cplBC.svOneD_solver_interface.has_data; const bool ci_set = bc_params->coupling_interface.value_set; - const bool ci_has_block = ci_set && bc_params->coupling_interface.svzerod_solver_block.defined(); - - if (svzd_iface) { - // Coupled BC to svZeroDSolver - lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_Coupled)); - lBc.bType = utils::ibclr(lBc.bType, enum_int(BoundaryConditionType::bType_bfs)); + const bool ci_has_block = ci_set && bc_params->coupling_interface.svzerod_solver_block.defined(); + const bool ci_has_1d_file = ci_set && bc_params->coupling_interface.svoned_input_file.defined(); // Sanity check: must define - if (!ci_has_block) { + if (svzd_iface || svOneD_iface) { + // svZeroD and/or svOneD path: route each face individually based on its + // content. + + if (ci_has_block && ci_has_1d_file) { + svmp::raise( + SVMP_HERE, + std::string("[read_bc] on face '") + face_name + + "' defines both and . " + "Specify exactly one per face.", + svmp::StatusCode::InvalidArgument); + } + + if (ci_has_block) { + // 0D face: route to svZeroD. + if (!svzd_iface) { + svmp::raise( + SVMP_HERE, + std::string("[read_bc] Face '") + face_name + + "' specifies but no " + "is defined on the equation.", + svmp::StatusCode::InvalidArgument); + } + lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_Coupled)); + lBc.bType = utils::ibclr(lBc.bType, enum_int(BoundaryConditionType::bType_Dir)); + lBc.bType = utils::ibclr(lBc.bType, enum_int(BoundaryConditionType::bType_Neu)); + lBc.bType = utils::ibclr(lBc.bType, enum_int(BoundaryConditionType::bType_bfs)); + lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_cpl)); + // oned_input_file stays empty — marks this as a 0D face in the + // CoupledBoundaryCondition construction block below. + + } else if (ci_has_1d_file) { + // 1D face: route to svOneD. + if (!svOneD_iface) { + svmp::raise( + SVMP_HERE, + std::string("[read_bc] Face '") + face_name + + "' specifies but no " + "is defined on the equation.", + svmp::StatusCode::InvalidArgument); + } + lBc.oned_input_file = bc_params->coupling_interface.svoned_input_file.value(); + + // svOneD now uses the CoupledBoundaryCondition interface (same flags as svZeroD). + lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_Coupled)); + lBc.bType = utils::ibclr(lBc.bType, enum_int(BoundaryConditionType::bType_Dir)); + lBc.bType = utils::ibclr(lBc.bType, enum_int(BoundaryConditionType::bType_Neu)); + lBc.bType = utils::ibclr(lBc.bType, enum_int(BoundaryConditionType::bType_bfs)); + lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_cpl)); + // cplBCptr stays -1; coupled_bc is constructed in the block below. + + if (com_mod.cplBC.schm == CplBCType::cplBC_NA) { + svmp::raise( + SVMP_HERE, + std::string("[read_bc] A coupling method (e.g. svOneDSolver_interface coupling_type) " + "must be defined for Time_dependence Coupled on face '") + + face_name + "'.", + svmp::StatusCode::InvalidArgument); + } + + } else { + // is present but has neither field, or is absent entirely. if (ci_set) { - throw std::runtime_error(std::string("[read_bc] on face '") + face_name + - "' must define ."); + svmp::raise( + SVMP_HERE, + std::string("[read_bc] on face '") + face_name + + "' must define either (for 0D coupling) or " + " (for 1D coupling).", + svmp::StatusCode::InvalidArgument); } - throw std::runtime_error( - std::string("[read_bc] With , each svZeroD-coupled face needs " - " with (Time_dependence Coupled) on face '") + - face_name + "'."); - } + if (svzd_iface) { + svmp::raise( + SVMP_HERE, + std::string("[read_bc] With , each svZeroD-coupled face needs " + " with " + "(Time_dependence Coupled) on face '") + + face_name + "'.", + svmp::StatusCode::InvalidArgument); + } else { + svmp::raise( + SVMP_HERE, + std::string("[read_bc] With , each 1D-coupled face needs " + " with " + "(Time_dependence Coupled) on face '") + + face_name + "'.", + svmp::StatusCode::InvalidArgument); + } + } } else { - // Coupled BC to GenBC - if (ci_set) { - throw std::runtime_error( - "[read_bc] is only valid when is defined on the equation."); + // genBC / cplBC path: no svZeroD or svOneD interface defined. + if (bc_params->coupling_interface.value_set) { + svmp::raise( + SVMP_HERE, + "[read_bc] is only valid when or " + " is defined on the equation.", + svmp::StatusCode::InvalidArgument); } lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_cpl)); @@ -284,25 +367,36 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, lBc.cplBCptr = com_mod.cplBC.nFa - 1; if (com_mod.cplBC.schm == CplBCType::cplBC_NA) { - throw std::runtime_error( + svmp::raise( + SVMP_HERE, std::string("[read_bc] A coupling method (e.g. Couple_to_genBC) must be defined for Time_dependence " "Coupled on face '") + - face_name + "'."); + face_name + "'.", + svmp::StatusCode::InvalidArgument); } } } else if (ctmp == "Resistance") { lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_res)); if (!utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Neu))) { - throw std::runtime_error("[read_bc] Resistance is only defined for Neu BC."); + svmp::raise( + SVMP_HERE, + "[read_bc] Resistance is only defined for Neu BC.", + svmp::StatusCode::InvalidArgument); } if (utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Robin))) { - throw std::runtime_error("[read_bc] Resistance is not defined for Robin BC."); + svmp::raise( + SVMP_HERE, + "[read_bc] Resistance is not defined for Robin BC.", + svmp::StatusCode::InvalidArgument); } if (std::set{Equation_fluid,Equation_FSI,Equation_CMM}.count(lEq.phys) == 0) { - throw std::runtime_error("[read_bc] Resistance is only defined for fluid/CMM/SI equations."); + svmp::raise( + SVMP_HERE, + "[read_bc] Resistance is only defined for fluid/CMM/SI equations.", + svmp::StatusCode::InvalidArgument); } lBc.r = bc_params->value.value(); @@ -327,8 +421,9 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, lBc.RCR.Pd = bc_params->rcr.distal_pressure.value(); lBc.RCR.Xo = bc_params->rcr.initial_pressure.value(); - if (com_mod.cplBC.schm != CplBCType::cplBC_NA || com_mod.cplBC.xo.size() != 0) { - throw std::runtime_error("[read_bc] RCR cannot be used in conjunction with cplBC."); + if ((com_mod.cplBC.schm != CplBCType::cplBC_NA && !com_mod.cplBC.useSvOneD && !com_mod.cplBC.useSvZeroD) || + com_mod.cplBC.xo.size() != 0) { + throw std::runtime_error("[read_bc] An RCR boundary condition can only be used when coupled with the svOneDsolver."); } com_mod.cplBC.nFa = com_mod.cplBC.nFa + 1; lBc.cplBCptr = com_mod.cplBC.nFa - 1; @@ -403,71 +498,138 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, } - // Coupled BC (svZeroDSolver / CoupledBoundaryCondition) + // Coupled BC (svZeroDSolver or svOneDSolver via CoupledBoundaryCondition) + bool is_coupled_dir = false; if (utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Coupled))) { - // Get asssociated face name and block name + // Get associated face name const auto& face_name = com_mod.msh[lBc.iM].fa[lBc.iFa].name; - const std::string zd_block = bc_params->coupling_interface.svzerod_solver_block.value(); - - // Get follower pressure load flag if defined - bool cpl_flwP = false; - if (lEq.phys == Equation_struct || lEq.phys == Equation_ustruct) { - cpl_flwP = bc_params->follower_pressure_load.value(); - } - // Sanity check: CoupledBoundaryCondition is only supported for struct, ustruct, fluid, FSI, or CMM physics - const auto cpl_phys = lEq.phys; - if (cpl_phys != Equation_struct && cpl_phys != Equation_ustruct && cpl_phys != Equation_fluid && - cpl_phys != Equation_FSI && cpl_phys != Equation_CMM) { + // Sanity check: CoupledBoundaryCondition must be Dirichlet or Neumann + if (coupled_bc_type != BoundaryConditionType::bType_Dir && + coupled_bc_type != BoundaryConditionType::bType_Neu) { throw std::runtime_error( - std::string("[read_bc] CoupledBoundaryCondition (svZeroDSolver) is only supported for struct, ustruct, fluid, FSI, or CMM physics on face '") + - face_name + "'."); - } - - // Sanity check: svZeroDSolver coupling is currently implemented only for Neumann-type boundaries. - if (!utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Neu))) { - throw std::runtime_error( - std::string("[read_bc] CoupledBoundaryCondition (svZeroDSolver) currently requires boundary Neu on face '") + + std::string("[read_bc] CoupledBoundaryCondition requires boundary Dirichlet or Neumann on face '") + face_name + "'."); } - // Sanity check: Follower pressure load must be used for 0D coupling with struct/ustruct - if ((cpl_phys == Equation_struct || cpl_phys == Equation_ustruct) && !cpl_flwP) { - throw std::runtime_error( - std::string("[read_bc] Follower pressure load must be used for 0D coupling with struct/ustruct on face '") + - face_name + "'."); - } + const bool is_svOneD_face = !lBc.oned_input_file.empty(); - // Get cap face VTP file name if defined - std::string zerod_cap; - bool use_cap = false; - if (bc_params->coupling_interface.chamber_cap_surface.defined()) { - zerod_cap = bc_params->coupling_interface.chamber_cap_surface.value(); - use_cap = true; - } + if (is_svOneD_face) { + // ------------------------------------------------------------------ + // svOneD path: construct CoupledBoundaryCondition with empty + // block_name and store the 1D input file path. + // ------------------------------------------------------------------ - // Figure out the coupled BC type - BoundaryConditionType coupled_bc_type = BoundaryConditionType::bType_Neu; - if (utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Dir))) { - coupled_bc_type = BoundaryConditionType::bType_Dir; - } else if (utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Neu))) { - coupled_bc_type = BoundaryConditionType::bType_Neu; - } + // Sanity check: only fluid, FSI, or CMM physics supported + const auto cpl_phys = lEq.phys; + if (cpl_phys != Equation_fluid && cpl_phys != Equation_FSI && cpl_phys != Equation_CMM) { + throw std::runtime_error( + std::string("[read_bc] CoupledBoundaryCondition (svOneDSolver) is only supported for " + "fluid, FSI, or CMM physics on face '") + + face_name + "'."); + } - // Create the coupled boundary condition object - if (use_cap) { lBc.coupled_bc = CoupledBoundaryCondition(coupled_bc_type, com_mod.msh[lBc.iM].fa[lBc.iFa], - com_mod.msh[lBc.iM].fa[lBc.iFa].name, zd_block, zerod_cap, - lEq.phys, cpl_flwP); + com_mod.msh[lBc.iM].fa[lBc.iFa].name, + /*block_name=*/"", lEq.phys, /*follower_pressure_load=*/false); + lBc.coupled_bc.set_oned_input_file(lBc.oned_input_file); + + // Read optional pressure ramp parameters for 1D coupling initialization. + if (bc_params->coupling_interface.coupling_ramp_steps.defined() && + bc_params->coupling_interface.coupling_ramp_steps.value() > 0) { + int ramp_steps = bc_params->coupling_interface.coupling_ramp_steps.value(); + double ramp_P_ref = bc_params->coupling_interface.coupling_ramp_ref_pressure.defined() + ? bc_params->coupling_interface.coupling_ramp_ref_pressure.value() + : 0.0; + lBc.coupled_bc.set_oned_ramp(ramp_steps, ramp_P_ref); + } + + // Read optional under-relaxation factor for DIR coupling pressure. + if (bc_params->coupling_interface.coupling_relax_factor.defined()) { + double omega = bc_params->coupling_interface.coupling_relax_factor.value(); + lBc.coupled_bc.set_oned_relax_factor(omega); + } + + } else if (com_mod.cplBC.svzerod_solver_interface.has_data) { + + // Get block name + const std::string zd_block = bc_params->coupling_interface.svzerod_solver_block.value(); + + // Get follower pressure load flag if defined + bool cpl_flwP = false; + if (lEq.phys == Equation_struct || lEq.phys == Equation_ustruct) { + cpl_flwP = bc_params->follower_pressure_load.value(); + } + + // Sanity check: CoupledBoundaryCondition is only supported for struct, ustruct, fluid, FSI, or CMM physics + const auto cpl_phys = lEq.phys; + if (cpl_phys != Equation_struct && cpl_phys != Equation_ustruct && cpl_phys != Equation_fluid && + cpl_phys != Equation_FSI && cpl_phys != Equation_CMM) { + throw std::runtime_error( + std::string("[read_bc] CoupledBoundaryCondition (svZeroDSolver) is only supported for struct, ustruct, fluid, FSI, or CMM physics on face '") + + face_name + "'."); + } + + // Sanity check: Follower pressure load must be used for 0D coupling with struct/ustruct + if ((cpl_phys == Equation_struct || cpl_phys == Equation_ustruct) && !cpl_flwP) { + throw std::runtime_error( + std::string("[read_bc] Follower pressure load must be used for 0D coupling with struct/ustruct on face '") + + face_name + "'."); + } + + // Get cap face VTP file name if defined + std::string zd_cap; + bool use_cap = false; + if (bc_params->coupling_interface.chamber_cap_surface.defined()) { + zd_cap = bc_params->coupling_interface.chamber_cap_surface.value(); + use_cap = true; + } + + if (use_cap) { + lBc.coupled_bc = CoupledBoundaryCondition(coupled_bc_type, com_mod.msh[lBc.iM].fa[lBc.iFa], + com_mod.msh[lBc.iM].fa[lBc.iFa].name, zd_block, zd_cap, + lEq.phys, cpl_flwP); + } else { + lBc.coupled_bc = CoupledBoundaryCondition(coupled_bc_type, com_mod.msh[lBc.iM].fa[lBc.iFa], + com_mod.msh[lBc.iM].fa[lBc.iFa].name, zd_block, lEq.phys, + cpl_flwP); + } + + // Read optional pressure ramp parameters for 0D coupling initialization. + if (bc_params->coupling_interface.coupling_ramp_steps.defined() && + bc_params->coupling_interface.coupling_ramp_steps.value() > 0) { + int ramp_steps = bc_params->coupling_interface.coupling_ramp_steps.value(); + double ramp_P_ref = bc_params->coupling_interface.coupling_ramp_ref_pressure.defined() + ? bc_params->coupling_interface.coupling_ramp_ref_pressure.value() + : 0.0; + lBc.coupled_bc.set_oned_ramp(ramp_steps, ramp_P_ref); + } + + // Read optional under-relaxation factor for 0D coupling. + if (bc_params->coupling_interface.coupling_relax_factor.defined()) { + double omega = bc_params->coupling_interface.coupling_relax_factor.value(); + lBc.coupled_bc.set_oned_relax_factor(omega); + } + } else { - lBc.coupled_bc = CoupledBoundaryCondition(coupled_bc_type, com_mod.msh[lBc.iM].fa[lBc.iFa], - com_mod.msh[lBc.iM].fa[lBc.iFa].name, zd_block, lEq.phys, - cpl_flwP); + throw std::runtime_error( + std::string("[read_bc] bType_Coupled is set on face '") + face_name + + "' but neither svZeroDSolver_interface nor svOneDSolver_interface is defined."); } - } + // For DIR Coupled BCs (svZeroD or svOneD), the downstream solver always returns a + // volumetric flow rate Q [m³/s], not a velocity [m/s]. Without bType_flx, bc_ini + // sets gx(a) = 1, and set_bc_dir_l applies velocity = Q * nV which has the wrong + // dimensions. With bType_flx, bc_ini normalises gx(a) = 1/area, so the applied + // velocity = (Q/area) * nV [m/s] is physically correct. Enforce this automatically + // regardless of the user's setting. + if (coupled_bc_type == BoundaryConditionType::bType_Dir) { + lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_flx)); + is_coupled_dir = true; + } + } // To impose value or flux bool ltmp = bc_params->impose_flux.value(); @@ -475,13 +637,22 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_flx)); } - // To zero-out perimeter or not. Default is .true. for Dir/CMM + // To zero-out perimeter or not. Default is .true. for Dir/CMM and Coupled-DIR. + // For Coupled-DIR BCs, bType_Dir is cleared, but we still need bType_zp so that + // bc_ini zeros shared perimeter nodes before normalising gx. Without bType_zp, + // gx = 1/full_area and the wall BC later zeros perimeter nodes, reducing the + // effective flux to Q*(interior_area/full_area) instead of Q. // ltmp = false; - ltmp = utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Dir)); + ltmp = utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Dir)) || is_coupled_dir; + if (bc_params->zero_out_perimeter.defined()) { ltmp = bc_params->zero_out_perimeter.value(); } + + if (is_coupled_dir) { + ltmp = true; + } lBc.bType = utils::ibclr(lBc.bType, enum_int(BoundaryConditionType::bType_zp)); if (ltmp || utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_CMM))) { @@ -1425,13 +1596,38 @@ void read_eq(Simulation* simulation, EquationParameters* eq_params, eqType& lEq) cplBC.useGenBC = true; cplbc_type_str = eq_params->couple_to_genBC.type.value(); - } else if (eq_params->svzerodsolver_interface_parameters.defined()) { - cplBC.useSvZeroD = true; - cplbc_type_str = eq_params->svzerodsolver_interface_parameters.coupling_type.value(); - cplBC.svzerod_solver_interface.set_data(eq_params->svzerodsolver_interface_parameters); + } else { + // Allow svZeroDSolver_interface and svOneDSolver_interface to coexist + // (mixed coupling: some faces go to 0D, others to 1D). + // + // Processing order: svZeroD first, then svOneD. + // cplbc_type_str is set by whichever interface is processed first; + // the second interface must match that type. + if (eq_params->svzerodsolver_interface_parameters.defined()) { + cplBC.useSvZeroD = true; + cplbc_type_str = eq_params->svzerodsolver_interface_parameters.coupling_type.value(); + cplBC.svzerod_solver_interface.set_data(eq_params->svzerodsolver_interface_parameters); + } + + if (eq_params->svonedsolver_interface_parameters.defined()) { + const std::string svOneD_type = eq_params->svonedsolver_interface_parameters.coupling_type.value(); + // When both interfaces are defined, their Coupling_type must match. + if (!cplbc_type_str.empty() && cplbc_type_str != svOneD_type) { + throw std::runtime_error( + "[read_eq] svZeroDSolver_interface and svOneDSolver_interface must use the same " + "Coupling_type in a mixed-coupling simulation " + "(svZeroDSolver_interface has '" + cplbc_type_str + + "', svOneDSolver_interface has '" + svOneD_type + "')."); + } + cplbc_type_str = svOneD_type; + cplBC.useSvOneD = true; + cplBC.svOneD_solver_interface.set_data(eq_params->svonedsolver_interface_parameters); + } } - if (eq_params->couple_to_genBC.defined() || eq_params->svzerodsolver_interface_parameters.defined()) { + if (eq_params->couple_to_genBC.defined() || + eq_params->svzerodsolver_interface_parameters.defined() || + eq_params->svonedsolver_interface_parameters.defined()) { try { cplBC.schm = consts::cplbc_name_to_type.at(cplbc_type_str); } catch (const std::out_of_range& exception) { @@ -1448,7 +1644,7 @@ void read_eq(Simulation* simulation, EquationParameters* eq_params, eqType& lEq) cplBC.nX = 0; cplBC.xp.resize(cplBC.nX); - } else if (cplBC.useSvZeroD) { + } else { cplBC.nX = 0; } } @@ -1537,9 +1733,13 @@ void read_eq(Simulation* simulation, EquationParameters* eq_params, eqType& lEq) if (std::set{Equation_fluid,Equation_FSI,Equation_CMM}.count(lEq.phys) == 0) { throw std::runtime_error("RCR-type BC is allowed for fluid/CMM/FSI eq. only."); } - cplBC.schm = CplBCType::cplBC_SI; - if (lEq.useTLS) { - cplBC.schm = CplBCType::cplBC_E; + // Only set coupling scheme if not already configured by an external solver (e.g. svOneD, svZeroD). + // When svOneD/svZeroD and RCR coexist, the external solver owns the scheme; RCR uses the same scheme. + if (!cplBC.useSvOneD && !cplBC.useSvZeroD) { + cplBC.schm = CplBCType::cplBC_SI; + if (lEq.useTLS) { + cplBC.schm = CplBCType::cplBC_E; + } } cplBC.nX = cplBC.nFa; cplBC.nXp = cplBC.nFa + 1; diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 01d361ceb..704030bb0 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -18,6 +18,7 @@ #include "utils.h" #include #include "svZeroD_interface.h" +#include "svOneD_interface.h" namespace set_bc { @@ -103,11 +104,18 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& RCRflag = true; } } - - // Compute flowrates at 3D Neumann0D boundaries at timesteps n and n+1 for Coupled BCs + + // Compute flowrates/pressures at 3D coupled boundaries for Coupled BCs. + // DIR coupling (both svZeroD and svOneD): the downstream solver is driven + // by the 3D face average pressure, so compute_pressures() is required. + // NEU coupling: the downstream solver is driven by the 3D outflow Q. if (utils::btest(bc.bType, iBC_Coupled)) { - bc.coupled_bc.compute_flowrates(com_mod, cm_mod, solutions); - #ifdef debug_calc_der_cpl_bc + if (bc.coupled_bc.get_bc_type() == consts::BoundaryConditionType::bType_Dir) { + bc.coupled_bc.compute_pressures(com_mod, cm_mod, solutions); + } else { + bc.coupled_bc.compute_flowrates(com_mod, cm_mod, solutions); + } + #ifdef debug_calc_der_cpl_bc dmsg << "iBC_Coupled "; dmsg << "coupled_bc.Qo: " << bc.coupled_bc.get_Qo(); dmsg << "coupled_bc.Qn: " << bc.coupled_bc.get_Qn(); @@ -171,10 +179,21 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& // Call genBC or cplBC to get updated pressures or flowrates. if (cplBC.useGenBC) { set_bc::genBC_Integ_X(com_mod, cm_mod, "D"); - } else if (cplBC.useSvZeroD) { - svZeroD::calc_svZeroD(com_mod, cm_mod, 'D'); } else { - set_bc::cplBC_Integ_X(com_mod, cm_mod, RCRflag); + // In mixed-coupling simulations both useSvZeroD and useSvOneD can be true. + // Call each active solver independently. + if (cplBC.useSvZeroD) { + svZeroD::calc_svZeroD(com_mod, cm_mod, 'D'); + } + if (cplBC.useSvOneD) { + svOneD::calc_svOneD(com_mod, cm_mod, 'D'); + } + if (!cplBC.useSvZeroD && !cplBC.useSvOneD) { + set_bc::cplBC_Integ_X(com_mod, cm_mod, RCRflag); + } else if (RCRflag) { + // Also integrate any RCR faces that coexist with svZeroD/svOneD faces. + set_bc::cplBC_Integ_X(com_mod, cm_mod, true); + } } // Compute the epsilon parameter (diff) for the finite difference calculation @@ -198,7 +217,7 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& diff = diff*relTol; } - // Store the original pressures and flowrates + // Store the original pressures and flowrates for cplBC std::vector orgY(cplBC.fa.size()); std::vector orgQ(cplBC.fa.size()); @@ -221,6 +240,16 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& set_bc::genBC_Integ_X(com_mod, cm_mod, "D"); } else if (cplBC.useSvZeroD) { svZeroD::calc_svZeroD(com_mod, cm_mod, 'D'); + // Also integrate any RCR faces that coexist with svZeroD faces. + if (RCRflag) { + set_bc::cplBC_Integ_X(com_mod, cm_mod, true); + } + } else if (cplBC.useSvOneD) { + svOneD::calc_svOneD(com_mod, cm_mod, 'D'); + // Also integrate any RCR faces that coexist with svOneD faces. + if (RCRflag) { + set_bc::cplBC_Integ_X(com_mod, cm_mod, true); + } } else { set_bc::cplBC_Integ_X(com_mod, cm_mod, RCRflag); } @@ -262,7 +291,11 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& // Perturb flowrate and compute new pressure bc.coupled_bc.perturb_flowrate(diff); - svZeroD::calc_svZeroD(com_mod, cm_mod, 'D'); + if (bc.coupled_bc.is_svOneD_face()) { + svOneD::calc_svOneD(com_mod, cm_mod, 'D'); + } else { + svZeroD::calc_svZeroD(com_mod, cm_mod, 'D'); + } // Finite difference: dP/dQ bc.r = (bc.coupled_bc.get_pressure() - orig_state.pressure) / diff; @@ -281,6 +314,7 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& } /// @brief RCR (Windkessel) integration for the non-genBC / non-svZeroD branch. +/// The legacy external Fortran Couple_to_cplBC file coupling has been removed. void cplBC_Integ_X(ComMod& com_mod, const CmMod& cm_mod, const bool RCRflag) { using namespace consts; @@ -476,26 +510,40 @@ void RCR_Integ_X(ComMod& com_mod, const CmMod& cm_mod, int istat) double tt = fmax(time - dt, 0.0); double dtt = dt / static_cast(nTS); - int nX = cplBC.nFa; + + // Collect indices of RCR faces only. When svOneD and RCR are mixed, some + // cplBC.fa[] entries belong to the 1D solver and must be skipped here to + // avoid accessing uninitialised RCR parameters (Rp/C/Rd/Pd = 0) and + // causing division-by-zero inside the RK4 loop. + std::vector rcrIdx; + rcrIdx.reserve(cplBC.nFa); + for (int i = 0; i < cplBC.nFa; i++) { + if (cplBC.fa[i].isRCR) { + rcrIdx.push_back(i); + } + } + int nX = static_cast(rcrIdx.size()); Vector Rp(nX), C(nX), Rd(nX), Pd(nX); Vector X(nX), Xrk(nX); Array frk(nX,4), Qrk(nX,4); - for (int i = 0; i < nX; i++) { - Rp(i) = cplBC.fa[i].RCR.Rp; - C(i) = cplBC.fa[i].RCR.C; - Rd(i) = cplBC.fa[i].RCR.Rd; - Pd(i) = cplBC.fa[i].RCR.Pd; + for (int k = 0; k < nX; k++) { + int i = rcrIdx[k]; + Rp(k) = cplBC.fa[i].RCR.Rp; + C(k) = cplBC.fa[i].RCR.C; + Rd(k) = cplBC.fa[i].RCR.Rd; + Pd(k) = cplBC.fa[i].RCR.Pd; + X(k) = cplBC.xo[i]; } - X = cplBC.xo; for (int n = 0; n < nTS; n++) { for (int i = 0; i < 4; i++) { double r = static_cast(i) / 3.0; r = (static_cast(n) + r) / static_cast(nTS); - for (int j = 0; j < Qrk.nrows(); j++) { - Qrk(j,i) = cplBC.fa[j].Qo + (cplBC.fa[j].Qn - cplBC.fa[j].Qo) * r; + for (int j = 0; j < nX; j++) { + int fi = rcrIdx[j]; + Qrk(j,i) = cplBC.fa[fi].Qo + (cplBC.fa[fi].Qn - cplBC.fa[fi].Qo) * r; } } @@ -512,21 +560,21 @@ void RCR_Integ_X(ComMod& com_mod, const CmMod& cm_mod, int istat) trk = tt + dtt / 3.0; Xrk = X + dtt * frk.col(0) / 3.0; - for (int j = 0; j < Qrk.nrows(); j++) { + for (int j = 0; j < nX; j++) { frk(j,1) = (Qrk(j,1) - (Xrk(j)-Pd(j)) / Rd(j)) / C(j); } // RK-4 3rd pass trk = tt + 2.0 * dtt / 3.0; Xrk = X - dtt * frk.col(0) / 3.0 + dtt * frk.col(1); - for (int j = 0; j < Qrk.nrows(); j++) { + for (int j = 0; j < nX; j++) { frk(j,2) = (Qrk(j,2) - (Xrk(j) - Pd(j)) / Rd(j)) / C(j); } // RK-4 4th pass trk = tt + dtt; Xrk = X + dtt * frk.col(0) - dtt * frk.col(1) + dtt * frk.col(2); - for (int j = 0; j < Qrk.nrows(); j++) { + for (int j = 0; j < nX; j++) { frk(j,3) = (Qrk(j,3) - (Xrk(j) - Pd(j)) / Rd(j)) / C(j); } @@ -534,8 +582,8 @@ void RCR_Integ_X(ComMod& com_mod, const CmMod& cm_mod, int istat) X = X + r*(frk.col(0) + 3.0*(frk.col(1) + frk.col(2)) + frk.col(3)); tt = tt + dtt; - for (int i = 0; i < nX; i++) { - if (isnan(X(i))) { + for (int k = 0; k < nX; k++) { + if (isnan(X(k))) { throw std::runtime_error("ERROR: NaN detected in RCR integration"); istat = -1; return; @@ -543,13 +591,14 @@ void RCR_Integ_X(ComMod& com_mod, const CmMod& cm_mod, int istat) } } - cplBC.xn = X; - cplBC.xp(0) = tt; - - for (int i = 0; i < nX; i++) { - cplBC.xp(i+1) = Qrk(i,3); //cplBC.fa(i).Qn - cplBC.fa[i].y = X(i) + (cplBC.fa[i].Qn * Rp(i)); + // Write results back using the original (sparse) face indices. + for (int k = 0; k < nX; k++) { + int i = rcrIdx[k]; + cplBC.xn[i] = X(k); + cplBC.xp(i+1) = Qrk(k,3); + cplBC.fa[i].y = X(k) + (cplBC.fa[i].Qn * Rp(k)); } + cplBC.xp(0) = tt; } @@ -743,10 +792,16 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) } } - - // Compute flowrates at 3D Neumann0D boundaries at timesteps n and n+1 for Coupled BCs + // Compute flowrates/pressures at 3D coupled boundaries for Coupled BCs. + // DIR coupling (both svZeroD and svOneD): the downstream solver is driven + // by the 3D face average pressure, so compute_pressures() is required. + // NEU coupling: the downstream solver is driven by the 3D outflow Q. if (utils::btest(bc.bType, iBC_Coupled)) { - bc.coupled_bc.compute_flowrates(com_mod, cm_mod, solutions); + if (bc.coupled_bc.get_bc_type() == consts::BoundaryConditionType::bType_Dir) { + bc.coupled_bc.compute_pressures(com_mod, cm_mod, solutions); + } else { + bc.coupled_bc.compute_flowrates(com_mod, cm_mod, solutions); + } } if (ptr != -1) { @@ -791,10 +846,21 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) // Updates pressure or flowrates stored in cplBC.fa[i].y if (cplBC.useGenBC) { set_bc::genBC_Integ_X(com_mod, cm_mod, "D"); - } else if (cplBC.useSvZeroD){ - svZeroD::calc_svZeroD(com_mod, cm_mod, 'D'); } else { - set_bc::cplBC_Integ_X(com_mod, cm_mod, RCRflag); + // In mixed-coupling simulations both useSvZeroD and useSvOneD can be true. + // Call each active solver independently. + if (cplBC.useSvZeroD) { + svZeroD::calc_svZeroD(com_mod, cm_mod, 'D'); + } + if (cplBC.useSvOneD) { + svOneD::calc_svOneD(com_mod, cm_mod, 'D'); + } + if (!cplBC.useSvZeroD && !cplBC.useSvOneD) { + set_bc::cplBC_Integ_X(com_mod, cm_mod, RCRflag); + } else if (RCRflag) { + // Also integrate any RCR faces that coexist with svZeroD/svOneD faces. + set_bc::cplBC_Integ_X(com_mod, cm_mod, true); + } } } @@ -802,9 +868,16 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) auto& bc = eq.bc[iBc]; int iFa = bc.iFa; - // For Coupled BC, get pressure from CoupledBoundaryCondition (set by svZeroD interface) + // For Coupled BC, get the coupling value from CoupledBoundaryCondition: + // NEU: 0D/1D solver returns pressure P → apply as Neumann traction + // DIR: 0D/1D solver returns flow rate Q → apply as Dirichlet velocity (bc.g = Q, then + // set_bc_dir_l multiplies by gx(a)*nV to produce the nodal velocity profile) if (utils::btest(bc.bType, iBC_Coupled)) { - bc.g = bc.coupled_bc.get_pressure(); + if (bc.coupled_bc.get_bc_type() == BoundaryConditionType::bType_Dir) { + bc.g = bc.coupled_bc.get_Qn(); + } else { + bc.g = bc.coupled_bc.get_pressure(); + } } // For other coupled BCs (Dir, Neu), get from cplBC.fa else { @@ -881,7 +954,12 @@ void set_bc_dir(ComMod& com_mod, SolutionStates& solutions) } } // END bType_CMM - if (!utils::btest(bc.bType, iBC_Dir)) { + // Allow Coupled BCs whose internal coupling type is DIR (svZeroD/svOneD DIR + // coupling): iBC_Dir is cleared for these in read_files but the velocity + // profile must still be applied via set_bc_dir_l. + bool isCoupledDir = utils::btest(bc.bType, iBC_Coupled) && + (bc.coupled_bc.get_bc_type() == BoundaryConditionType::bType_Dir); + if (!utils::btest(bc.bType, iBC_Dir) && !isCoupledDir) { continue; } @@ -1390,7 +1468,11 @@ void set_bc_neu(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& solu if (utils::btest(bc.bType, iBC_Ris0D)) {continue;} - if (utils::btest(bc.bType, iBC_Neu) || utils::btest(bc.bType, iBC_Coupled)) { + // Coupled BCs with DIR type must be handled by set_bc_dir (velocity profile), + // not here. Only NEU Coupled BCs get a Neumann pressure traction. + bool isCoupledDir = utils::btest(bc.bType, iBC_Coupled) && + (bc.coupled_bc.get_bc_type() == BoundaryConditionType::bType_Dir); + if ((utils::btest(bc.bType, iBC_Neu) || utils::btest(bc.bType, iBC_Coupled)) && !isCoupledDir) { #ifdef debug_set_bc_neu dmsg << "iM: " << iM+1; dmsg << "iFa: " << iFa+1; @@ -1459,7 +1541,38 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const } } - } else if (utils::btest(lBc.bType,iBC_res)) { + } else if (utils::btest(lBc.bType,iBC_Coupled)) { + // New-style Coupled NEU BC (svOneD/svZeroD): apply the actual 1D/0D + // pressure from the most recent 'D' solver call. bc.g is updated + // every Newton iteration by set_bc_cpl / calc_der_cpl_bc. + + //h(0) = lBc.g; + + double Q_3D = all_fun::integ(com_mod, cm_mod, lFa, Yn, eq.s, solutions, + eq.s+nsd-1, false, + consts::MechanicalConfigurationType::reference); + + + h(0) = lBc.g; + //h(0) = lBc.g - lBc.r * std::abs(Q_3D); + + // Backflow kinetic energy correction: when backflow is detected + // (Q < 0), subtract the face-averaged dynamic pressure to further + // reduce the applied traction and damp the incoming flow. + if (Q_3D < 0.0) { + int iM = lFa.iM; + int cDmn_local = all_fun::domain(com_mod, com_mod.msh[iM], cEq, lFa.gE(0)); + double rho = eq.dmn[cDmn_local].prop.at( + consts::PhysicalProperyType::fluid_density); + double beta = eq.dmn[cDmn_local].prop.at( + consts::PhysicalProperyType::backflow_stab); + double A = lFa.area; + if (A > 0.0) { + double u_n = Q_3D / A; // face-averaged normal velocity (< 0) + h(0) -= 0.5 * beta * rho * u_n * u_n; + } + } + } else if (utils::btest(lBc.bType,iBC_res)) { h(0) = lBc.r * all_fun::integ(com_mod, cm_mod, lFa, Yn, eq.s, solutions, eq.s+nsd-1, false, consts::MechanicalConfigurationType::reference); } else if (utils::btest(lBc.bType,iBC_std)) { diff --git a/Code/Source/solver/svOneD_interface.cpp b/Code/Source/solver/svOneD_interface.cpp new file mode 100644 index 000000000..c0004759c --- /dev/null +++ b/Code/Source/solver/svOneD_interface.cpp @@ -0,0 +1,460 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +// 3D-1D coupling subroutines. +// +// These routines interface the 3D finite-element solver (svMultiPhysics) with +// the 1D blood-flow solver (svOneDSolver) via a dynamically loaded shared +// library (libsvOneDSolver_interface.so/.dylib). +// +// Coupling overview +// ----------------- +// NEU coupling (1D inlet driven by 3D outflow): +// 3D → 1D : flow rate Q (as params[3..4]) +// 1D → 3D : pressure P (as cpl_value) +// BC applied on 3D face : Neumann (pressure traction) +// +// DIR coupling (1D outlet driven by 3D pressure): +// 3D → 1D : pressure P (as params[3..4]) +// 1D → 3D : flow rate Q (as cpl_value) +// BC applied on 3D face : Dirichlet (velocity profile) +// +// Parallelism model +// ----------------- +// Unlike the 0D solver (which is solved once on the master rank), each +// 1D model is INDEPENDENT and has its own input file. Multiple 1D models +// are therefore read, initialized, and solved in parallel: +// +// Initialization (init_svOneD): +// Phase 1 – parallel init: +// - Collect all svOneD-coupled faces into a list indexed 0..N-1. +// - Assign face (model) k to MPI rank k % nProcs. +// - Each rank reads and initializes ONLY its owned model(s) with no +// MPI synchronization, so all ranks work simultaneously. +// Phase 2 – batch metadata exchange: +// - After every rank has finished initializing its own model(s), +// share system_size and coupled_dof via MPI_Bcast so that all +// ranks know the sizes needed for subsequent result broadcasts. +// +// Time-stepping (calc_svOneD): +// Phase 1 – parallel solve: +// - Each rank runs run_simulation() for its owned model(s) with no MPI +// calls, so model k on rank A and model k+1 on rank B truly run +// concurrently. +// Phase 2 – batch result exchange: +// - After every rank has finished solving, results are shared via +// MPI_Bcast so that ALL ranks know the result, and each BC's +// coupled_bc.set_pressure() is updated accordingly. +// +// params array passed to run_1d_simulation_step_1d_: +// params[0] = 2.0 (number of time points) +// params[1] = t_old (time at start of step) +// params[2] = t_new (time at end of step) +// params[3] = BC_val_old (Q or P at t_old) +// params[4] = BC_val_new (Q or P at t_new) + +#include "svOneD_interface.h" + +#include +#include +#include +#include +#include +#include +#include + +#include "ComMod.h" +#include "consts.h" +#include "utils.h" +#include "svOneD_interface/OneDSolverInterface.h" +#include "Core/Exception.h" + +#include "mpi.h" + +namespace svOneD { + +// --------------------------------------------------------------------------- +// Per-model state. Each entry corresponds to one svOneD-coupled face (one 1D +// model). Indexed by the sequential order in which coupled faces were found +// in eq[0].bc[]. +// --------------------------------------------------------------------------- + +struct OneDModelState { + // Interface object (null on ranks that do not own this model). + OneDSolverInterface* interface = nullptr; + + // Problem identifier returned by the 1D library. + int problem_id = 0; + + // Total DOF count (nodes × 2). + int system_size = 0; + + // "NEU" or "DIR". + std::string coupling_type; + + // Index in the solution vector that corresponds to the coupled BC DOF. + int coupled_dof = 0; + + // Authoritative solution vector (only valid on the owning rank). + std::vector solution; + + // Owning MPI rank for this model. + int owner_rank = 0; + + // Index into eq[0].bc[] for the BC this model services. + int iBc = -1; + + // Pressure ramp for 1D coupling initialization. + // DIR: over the first ramp_steps committed time steps the pressure sent to + // the 1D solver is linearly interpolated from ramp_ref_pressure to the + // actual 3D pressure value. + // NEU: over the first ramp_steps committed time steps the flow rate sent to + // the 1D solver is linearly interpolated from 0 to the actual 3D Q. + // Zero means no ramping. + int ramp_steps = 0; + double ramp_ref_pressure = 0.0; + int step_count = 0; ///< Number of committed (BCFlag=='L') steps taken. + + // Under-relaxation (omega in (0, 1]). Default 1.0 = no relaxation. + // + // DIR coupling: + // Input (P sent to 1D): P_sent = omega * P_target + (1-omega) * P_prev_sent + // Output (Q from 1D) : Q_relax = omega * Q_raw + (1-omega) * Q_prev_sent + // + // NEU coupling: + // Output (P from 1D) : P_applied = omega * P_target + (1-omega) * P_neu_prev + // (P_target already includes the ramp from ramp_ref_pressure to P_raw) + double relax_factor = 1.0; + double P_prev_sent_old = 0.0; ///< Under-relaxed pressure sent at params[3] (t_old) on last 'L' step (DIR). + double P_prev_sent_new = 0.0; ///< Under-relaxed pressure sent at params[4] (t_new) on last 'L' step (DIR). + double Q_prev_sent = 0.0; ///< Under-relaxed flow rate output on last 'L' step (DIR only). + double P_neu_prev = 0.0; ///< Under-relaxed pressure applied to 3D on last 'L' step (NEU only). + double Q_prev_sent_old = 0.0; ///< Under-relaxed Q sent at params[3] (t_old) on last 'L' step (NEU). + double Q_prev_sent_new = 0.0; ///< Under-relaxed Q sent at params[4] (t_new) on last 'L' step (NEU). +}; + +// --------------------------------------------------------------------------- +// Module-level state. +// --------------------------------------------------------------------------- + +// One entry per svOneD-coupled face, filled during init_svOneD(). +static std::vector oned_models; + +// Shared library handle (one per process, loaded once). +static OneDSolverInterface* shared_lib_instance = nullptr; + +// Simulation time (advanced only on 'L' steps). +static double svOneDTime = 0.0; + +// --------------------------------------------------------------------------- +// Helper: resolve the shared-library path (.so / .dylib / as-is). +// --------------------------------------------------------------------------- +static std::string resolve_lib_path(const std::string& lib_base) +{ + if (std::filesystem::is_regular_file(lib_base + ".so")) { + return lib_base + ".so"; + } + + if (std::filesystem::is_regular_file(lib_base + ".dylib")) { + return lib_base + ".dylib"; + } + return lib_base; // already has extension, or will fail at dlopen time +} + +// --------------------------------------------------------------------------- +// init_svOneD +// --------------------------------------------------------------------------- +void init_svOneD(ComMod& com_mod, const CmMod& cm_mod) +{ + using namespace consts; + + auto& cplBC = com_mod.cplBC; + auto& solver_if = cplBC.svOneD_solver_interface; + auto& cm = com_mod.cm; + const int nProcs = cm.nProcs; + const int myRank = cm.taskId; + + if (!solver_if.has_data) { + throw svmp::CoreException( + "[svOneD::init_svOneD] svOneD solver interface data is missing. Please check the XML input file.", + svmp::StatusCode::InvalidState, + __FILE__, + __LINE__, + __func__); + } + + // Initialize the 1D simulation clock from the 3D solver's current time so + // that restarts and non-zero start times are handled correctly. + svOneDTime = com_mod.time; + + // ----- Collect the list of svOneD-coupled faces ----- + // Iterate over eq[0]'s BCs and pick those with iBC_Coupled and a non-empty + // oned_input_file (stored in coupled_bc). + { + const int iEq = 0; + const auto& eq = com_mod.eq[iEq]; + for (int iBc = 0; iBc < eq.nBc; iBc++) { + const auto& bc = eq.bc[iBc]; + if (!utils::btest(bc.bType, iBC_Coupled)) continue; + if (bc.coupled_bc.get_oned_input_file().empty()) continue; + + OneDModelState st; + st.iBc = iBc; + st.coupling_type = (bc.coupled_bc.get_bc_type() == BoundaryConditionType::bType_Neu) ? "NEU" : "DIR"; + st.ramp_steps = bc.coupled_bc.get_oned_ramp_steps(); + st.ramp_ref_pressure = bc.coupled_bc.get_oned_ramp_ref_pressure(); + st.relax_factor = bc.coupled_bc.get_oned_relax_factor(); + oned_models.push_back(std::move(st)); + } + } + + if (oned_models.empty()) { + throw svmp::CoreException( + "[svOneD::init_svOneD] No svOneD-coupled faces with input files found. Please check the XML input file.", + svmp::StatusCode::InvalidState, + __FILE__, + __LINE__, + __func__); + } + + // ----- Guard: require at least one MPI rank per 1D model ----- + // Each rank owns exactly one model (owner_rank = k % nProcs). If nProcs < N + // a single rank would own multiple models and call shared_lib_instance->initialize() + // more than once, corrupting the static problem-ID state inside the shared library. + const int nTotalModels = static_cast(oned_models.size()); + if (nProcs < nTotalModels) { + throw svmp::CoreException( + "[svOneD::init_svOneD] Number of MPI processes (" + std::to_string(nProcs) + + ") is less than the number of svOneD-coupled faces (" + + std::to_string(nTotalModels) + + "). Please run with at least " + std::to_string(nTotalModels) + + " MPI processes.", + svmp::StatusCode::InvalidState, + __FILE__, + __LINE__, + __func__); + } + + // ----- Load shared library (once per process) ----- + const std::string lib_path = resolve_lib_path(solver_if.solver_library); + shared_lib_instance = new OneDSolverInterface(); + shared_lib_instance->load_library(lib_path); + + // ----- Assign ranks and initialize owned models (Phase 1: parallel) ----- + // No MPI calls in this loop. All ranks proceed simultaneously, each + // reading and initializing only the model(s) it owns. Rank k owns model k + // (assigned via k % nProcs), so for N models and N ranks every rank handles + // exactly one model with no inter-rank synchronization. + const int iEq = 0; + auto& eq = com_mod.eq[iEq]; + for (int k = 0; k < nTotalModels; k++) { + auto& st = oned_models[k]; + st.owner_rank = k % nProcs; + + if (myRank != st.owner_rank) continue; + + // This rank owns model k: read the input file and initialize. + const std::string& input_file = eq.bc[st.iBc].coupled_bc.get_oned_input_file(); + int problem_id = 0; + int system_size = 0; + + shared_lib_instance->initialize(input_file, problem_id, system_size, + st.coupling_type); + st.problem_id = problem_id; + st.system_size = system_size; + st.interface = shared_lib_instance; + + shared_lib_instance->set_external_step_size(problem_id, com_mod.dt); + shared_lib_instance->extract_coupled_dof(problem_id, st.coupled_dof, + st.coupling_type); + + st.solution.resize(system_size, 0.0); + shared_lib_instance->return_solution(problem_id, st.solution.data(), system_size); + + // Initial coupled value = 0; first calc_svOneD call sets the real value. + eq.bc[st.iBc].coupled_bc.set_pressure(0.0); + } + + // ----- Broadcast metadata for all models (Phase 2: batch exchange) ----- + // All initialization is complete. Now share system_size and coupled_dof + // from each owner so that every rank knows the sizes needed for consistent + // result broadcasts in calc_svOneD. + for (int k = 0; k < nTotalModels; k++) { + auto& st = oned_models[k]; + MPI_Bcast(&st.system_size, 1, MPI_INT, st.owner_rank, cm.com()); + MPI_Bcast(&st.coupled_dof, 1, MPI_INT, st.owner_rank, cm.com()); + } + + // Run one 'D' step to populate the initial resistance term bc.r. + if (cplBC.schm != CplBCType::cplBC_E) { + calc_svOneD(com_mod, cm_mod, 'D'); + } +} + +// --------------------------------------------------------------------------- +// calc_svOneD +// --------------------------------------------------------------------------- +void calc_svOneD(ComMod& com_mod, const CmMod& cm_mod, char BCFlag) +{ + using namespace consts; + + auto& cm = com_mod.cm; + const int myRank = cm.taskId; + + const double t_old = svOneDTime; + const double t_new = svOneDTime + com_mod.dt; + const int nTotalModels = static_cast(oned_models.size()); + + const int iEq = 0; + auto& eq = com_mod.eq[iEq]; + + // ----- Phase 1: each rank runs its own models without blocking ----- + // All ranks proceed through this loop simultaneously, each executing only + // the models it owns. No MPI call here, so model k on rank A and model k+1 + // on rank B truly run at the same time. + std::vector cpl_values(nTotalModels, 0.0); + + for (int k = 0; k < nTotalModels; k++) { + auto& st = oned_models[k]; + auto& bc = eq.bc[st.iBc]; + + if (myRank != st.owner_rank) continue; + + // Build params = [2, t_old, t_new, BC_val_old, BC_val_new] + double params[5]; + params[0] = 2.0; + params[1] = t_old; + params[2] = t_new; + + if (bc.coupled_bc.get_bc_type() == BoundaryConditionType::bType_Neu) { + // NEU coupling: apply under-relaxation to the 3D flow rate Q sent to + // the 1D solver to damp timestep-to-timestep oscillations in the input. + // No ramping is applied here; ramping is applied to the *output* pressure + // P that the 1D solver returns (Phase 2 below). + const double omega = st.relax_factor; + params[3] = omega * bc.coupled_bc.get_Qo() + (1.0 - omega) * st.Q_prev_sent_new; + params[4] = omega * bc.coupled_bc.get_Qn() + (1.0 - omega) * st.Q_prev_sent_new; + } else { + double raw_P_old = bc.coupled_bc.get_Po(); + double raw_P_new = bc.coupled_bc.get_Pn(); + + // Step 1: apply pressure ramp (scales amplitude from ramp_ref_pressure + // to actual 3D pressure over the first ramp_steps steps). + double P_target_old, P_target_new; + if (st.ramp_steps > 0) { + double ramp_factor = std::min(1.0, static_cast(st.step_count) / st.ramp_steps); + double P_ref = st.ramp_ref_pressure; + P_target_old = P_ref + ramp_factor * (raw_P_old - P_ref); + P_target_new = P_ref + ramp_factor * (raw_P_new - P_ref); + } else { + P_target_old = raw_P_old; + P_target_new = raw_P_new; + } + + // Step 2: apply under-relaxation (damps timestep-to-timestep oscillations). + // P_sent = omega * P_target + (1 - omega) * P_prev_sent + const double omega = st.relax_factor; + params[3] = omega * P_target_old + (1.0 - omega) * st.P_prev_sent_new; + params[4] = omega * P_target_new + (1.0 - omega) * st.P_prev_sent_new; + } + + // Working copy of solution so that 'D' steps don't corrupt the + // committed state. + std::vector work_sol = st.solution; + st.interface->update_solution(st.problem_id, work_sol.data(), st.system_size); + + int save_incr = com_mod.saveIncr; + int error_code = 0; + + try { + st.interface->run_simulation(st.problem_id, t_old, save_incr, + st.coupling_type, params, + work_sol.data(), cpl_values[k], BCFlag, error_code); + } catch (const std::exception& e) { + svmp::raise( + SVMP_HERE, + "[svOneD::calc_svOneD] 1D solver step failed for svOneD input file '" + + bc.coupled_bc.get_oned_input_file() + "' with coupling type '" + + st.coupling_type + "' at time " + std::to_string(t_new) + ". " + + "Original error: " + e.what()); + } + + // Commit the updated solution only on the final iteration. + if (BCFlag == 'L') { + st.solution = work_sol; + // Update the under-relaxation history with the values actually sent. + // NEU: Q_prev_sent updated here; P history updated in Phase 2. + // DIR: P_prev_sent tracks the pressure value sent to the 1D solver. + if (st.coupling_type == "NEU") { + st.Q_prev_sent_old = params[3]; + st.Q_prev_sent_new = params[4]; + } else { + st.P_prev_sent_old = params[3]; + st.P_prev_sent_new = params[4]; + } + } + } + + // ----- Phase 2: broadcast all results and update coupled BCs ----- + // After every rank has finished solving its own models, gather the results. + // Each MPI_Bcast here is a cheap scalar transfer; the expensive 1D solver + // work has already been done concurrently in Phase 1. + for (int k = 0; k < nTotalModels; k++) { + auto& st = oned_models[k]; + MPI_Bcast(&cpl_values[k], 1, MPI_DOUBLE, st.owner_rank, cm.com()); + auto& cpl_bc = eq.bc[st.iBc].coupled_bc; + if (cpl_bc.get_bc_type() == BoundaryConditionType::bType_Dir) { + // 1D solver returns flow Q for DIR coupling; store it as flowrate so that + // set_bc can read get_Qn() and build the nodal velocity profile. + // Negate the sign: the 1D solver returns Q > 0 for inflow, but the 3D + // code applies velocity as Q * gx * outward_normal, so Q must be negative + // for an inlet face (outward normal points away from the domain). + // This matches the svZeroD convention: in_out = -1 for DIR (outlet of 0D + // = inlet of 3D), giving QCoupled = -1 * lpn_state_y[flow_id]. + double Q_raw = -cpl_values[k]; + // Apply under-relaxation to the Q output to damp timestep-to-timestep oscillations. + const double omega = st.relax_factor; + double Qn_relaxed = omega * Q_raw + (1.0 - omega) * st.Q_prev_sent; + double Qo_prev = cpl_bc.get_Qn(); + cpl_bc.set_flowrates(Qo_prev, Qn_relaxed); + if (BCFlag == 'L') { + st.Q_prev_sent = Qn_relaxed; + } + } else { + // 1D solver returns pressure P for NEU coupling. + // + // Step 1: apply pressure ramp (scales output P from ramp_ref_pressure + // to the actual 1D pressure over the first ramp_steps committed + // steps). This prevents a large sudden pressure jump from being + // imposed on the 3D domain at startup, which is the primary cause + // of oscillations in Neumann coupling. + double P_raw = cpl_values[k]; + double P_target; + if (st.ramp_steps > 0) { + double ramp_factor = std::min(1.0, static_cast(st.step_count) / st.ramp_steps); + double P_ref = st.ramp_ref_pressure; + P_target = P_ref + ramp_factor * (P_raw - P_ref); + } else { + P_target = P_raw; + } + // Step 2: apply under-relaxation to damp timestep-to-timestep oscillations. + // P_applied = omega * P_target + (1 - omega) * P_prev_applied + const double omega = st.relax_factor; + double P_relaxed = omega * P_target + (1.0 - omega) * st.P_neu_prev; + cpl_bc.set_pressure(P_relaxed); + if (BCFlag == 'L') { + st.P_neu_prev = P_relaxed; + } + } + } + + // Advance the simulation clock after the final iteration. + if (BCFlag == 'L') { + svOneDTime += com_mod.dt; + for (auto& st : oned_models) { + st.step_count++; + } + } +} + +} // namespace svOneD diff --git a/Code/Source/solver/svOneD_interface.h b/Code/Source/solver/svOneD_interface.h new file mode 100644 index 000000000..6b8009cd0 --- /dev/null +++ b/Code/Source/solver/svOneD_interface.h @@ -0,0 +1,25 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef SV1D_SUBROUTINES_H +#define SV1D_SUBROUTINES_H + +#include "Simulation.h" +#include "consts.h" +#include "svOneD_interface/OneDSolverInterface.h" + +namespace svOneD { + +/// @brief Initialize the 1D solver and populate the initial cplBC state. +/// Called once from baf_ini() after the BC data structures are set up. +void init_svOneD(ComMod& com_mod, const CmMod& cm_mod); + +/// @brief Advance the 1D solver by one time step and update the coupled BC value. +/// +/// @param BCFlag 'D' - derivative / perturbation step (state is NOT committed). +/// 'L' - last Newton iteration (state IS committed, time advances). +void calc_svOneD(ComMod& com_mod, const CmMod& cm_mod, char BCFlag); + +} // namespace svOneD + +#endif // SV1D_SUBROUTINES_H diff --git a/Code/Source/solver/svOneD_interface/OneDSolverInterface.cpp b/Code/Source/solver/svOneD_interface/OneDSolverInterface.cpp new file mode 100644 index 000000000..24f193263 --- /dev/null +++ b/Code/Source/solver/svOneD_interface/OneDSolverInterface.cpp @@ -0,0 +1,176 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +#include "OneDSolverInterface.h" + +#include +#include +#include +#include +#include "Core/Exception.h" + +OneDSolverInterface::~OneDSolverInterface() +{ + if (library_handle_) { + dlclose(library_handle_); + library_handle_ = nullptr; + } +} + +void OneDSolverInterface::load_library(const std::string& interface_lib) +{ + library_handle_ = dlopen(interface_lib.c_str(), RTLD_LAZY); + if (!library_handle_) { + throw svmp::CoreException( + std::string("[OneDSolverInterface] Could not load shared library '") + + interface_lib + "': " + dlerror(), + svmp::StatusCode::DependencyError, + __FILE__, + __LINE__, + __func__); + } + + // Clear any existing error. + dlerror(); + + auto load_sym = [&](const char* name) -> void* { + void* sym = dlsym(library_handle_, name); + const char* err = dlerror(); + if (err) { + throw svmp::CoreException( + std::string("[OneDSolverInterface] Could not find the 1D interface function '") + + name + "' in the shared library '" + interface_lib + "': " + err, + svmp::StatusCode::DependencyError, + __FILE__, + __LINE__, + __func__); + } + return sym; + }; + + *(void**)(&initialize_1d_) = load_sym("initialize_1d"); + *(void**)(&set_external_step_size_1d_) = load_sym("set_external_step_size_1d"); + *(void**)(&return_1d_solution_) = load_sym("return_1d_solution"); + *(void**)(&update_1d_solution_) = load_sym("update_1d_solution"); + *(void**)(&run_1d_simulation_step_1d_) = load_sym("run_1d_simulation_step_1d"); + *(void**)(&extract_coupled_dof_) = load_sym("extract_coupled_dof"); +} + +void OneDSolverInterface::initialize(const std::string& input_file, + int& problem_id, + int& system_size, + const std::string& coupling_type) +{ + if (!initialize_1d_) { + throw svmp::CoreException( + "[OneDSolverInterface] Cannot initialize the 1D solver because the " + "1D interface library has not been loaded. Call load_library() with " + "the svOneDSolver interface shared library before initialize(), and " + "check that the library path is correct.", + svmp::StatusCode::DependencyError, + __FILE__, + __LINE__, + __func__); + } + initialize_1d_(input_file.c_str(), problem_id, system_size, + coupling_type.c_str()); + problem_id_ = problem_id; + system_size_ = system_size; +} + +void OneDSolverInterface::set_external_step_size(int problem_id, double dt) +{ + if (!set_external_step_size_1d_) { + throw svmp::CoreException( + "[OneDSolverInterface] set_external_step_size_1d not loaded", + svmp::StatusCode::DependencyError, + __FILE__, + __LINE__, + __func__); + } + set_external_step_size_1d_(problem_id, dt); +} + +void OneDSolverInterface::return_solution(int problem_id, double* solution, int size) +{ + if (!return_1d_solution_) { + throw svmp::CoreException( + "[OneDSolverInterface] return_1d_solution not loaded", + svmp::StatusCode::DependencyError, + __FILE__, + __LINE__, + __func__); + } + return_1d_solution_(problem_id, solution, size); +} + +void OneDSolverInterface::update_solution(int problem_id, double* solution, int size) +{ + if (!update_1d_solution_) { + throw svmp::CoreException( + "[OneDSolverInterface] update_1d_solution not loaded", + svmp::StatusCode::DependencyError, + __FILE__, + __LINE__, + __func__); + } + update_1d_solution_(problem_id, solution, size); +} + +void OneDSolverInterface::run_simulation(int problem_id, double current_time, + int save_incr, + const std::string& coupling_type, + double* params, double* solution, + double& cpl_value, char last_flag, + int& error_code) +{ + if (!run_1d_simulation_step_1d_) { + throw svmp::CoreException( + "[OneDSolverInterface] Cannot run the 1D solver because the " + "run_1d_simulation_step_1d interface function is not available. " + "Check that the svOneDSolver interface shared library was loaded correctly.", + svmp::StatusCode::DependencyError, + __FILE__, + __LINE__, + __func__); + } + + std::vector ctype_buf(coupling_type.begin(), coupling_type.end()); + ctype_buf.push_back('\0'); + + char flag_buf[2] = { last_flag, '\0' }; + + error_code = 0; + + run_1d_simulation_step_1d_(problem_id, current_time, save_incr, + ctype_buf.data(), params, solution, + cpl_value, flag_buf, error_code); + + if (error_code != 0) { + throw svmp::CoreException( + "[OneDSolverInterface] svOneDSolver failed while advancing the 1D model. " + "The 1D solver returned error code " + std::to_string(error_code) + ".", + svmp::StatusCode::DependencyError, + __FILE__, + __LINE__, + __func__); + } +} + +void OneDSolverInterface::extract_coupled_dof(int problem_id, int& coupled_dof, + const std::string& coupling_type) +{ + if (!extract_coupled_dof_) { + throw svmp::CoreException( + "[OneDSolverInterface] extract_coupled_dof not loaded", + svmp::StatusCode::DependencyError, + __FILE__, + __LINE__, + __func__); + } + // Copy into a mutable buffer; the shared-library function signature uses + // char* (not const char*) so we must pass a writable copy. + std::vector buf(coupling_type.begin(), coupling_type.end()); + buf.push_back('\0'); + extract_coupled_dof_(problem_id, coupled_dof, buf.data()); +} diff --git a/Code/Source/solver/svOneD_interface/OneDSolverInterface.h b/Code/Source/solver/svOneD_interface/OneDSolverInterface.h new file mode 100644 index 000000000..a13dd30d0 --- /dev/null +++ b/Code/Source/solver/svOneD_interface/OneDSolverInterface.h @@ -0,0 +1,117 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef ONEDSOLVER_INTERFACE_H +#define ONEDSOLVER_INTERFACE_H + +#include +#include +#include + +/// @brief Wrapper class for dynamically loading and calling the 1D solver shared library. +/// +/// This class loads the svOneDSolver shared library and +/// provides a C++ wrapper around its exported C interface. +/// +/// The interface is used to couple 3D Navier-Stokes simulations in +/// svMultiPhysics with reduced-order 1D blood flow models in +/// svOneDSolver. +/// +/// Supported coupling modes: +/// - NEU coupling: 3D flow rate -> 1D model, 1D pressure -> 3D solver +/// - DIR coupling: 3D pressure -> 1D model, 1D flow rate -> 3D solver +/// +/// Each initialized 1D model is identified by a problem_id assigned +/// by svOneDSolver and used in subsequent library calls. +/// +/// Shared library functions: +/// - initialize_1d(input_file, problem_id, system_size, coupling_type) +/// - set_external_step_size_1d(problem_id, dt) +/// - return_1d_solution(problem_id, solution, size) +/// - update_1d_solution(problem_id, solution, size) +/// - run_1d_simulation_step_1d(problem_id, time, save_flag, coupling_type, +/// params, solution, cpl_value, error_code) +/// - extract_coupled_dof(problem_id, coupled_dof, coupling_type) +// +class OneDSolverInterface { + public: + OneDSolverInterface() = default; + ~OneDSolverInterface(); + + /// @brief Load the 1D solver shared library from the given path. + void load_library(const std::string& interface_lib); + + /// @brief Initialize the 1D solver from an input file. + /// + /// @param[in] input_file Path to the 1D solver .in file. + /// @param[out] problem_id Problem identifier assigned by the solver. + /// @param[out] system_size Total number of DOFs (nodes * 2: flow + area). + /// @param[in] coupling_type "NEU" or "DIR" coupling direction. + void initialize(const std::string& input_file, int& problem_id, + int& system_size, const std::string& coupling_type); + + /// @brief Synchronize the 1D solver's internal time step with the 3D solver. + void set_external_step_size(int problem_id, double dt); + + /// @brief Copy the current 1D solution into the caller-provided buffer. + void return_solution(int problem_id, double* solution, int size); + + /// @brief Push a solution vector into the 1D solver as the current state. + void update_solution(int problem_id, double* solution, int size); + + /// @brief Advance the 1D solver by one time step. + /// + /// @param[in] problem_id Problem identifier. + /// @param[in] save_incr VTK output interval (Increment_in_saving_VTK_files from solver.xml); + /// the 1D library decides internally whether to write output. + /// @param[in] coupling_type "NEU" or "DIR". + /// @param[in] params Array [N, t1, t2, ..., val1, val2, ...] where N=2. + /// @param[in,out] solution Solution vector updated after the step. + /// @param[out] cpl_value The BC value returned by the 1D solver + /// (pressure for NEU, flow for DIR). + /// @param[in] last_flag 'L' for the final (committed) iteration, 'D' for + /// derivative / predictor steps. + /// @param[out] error_code Non-zero on failure. + void run_simulation(int problem_id, double current_time, int save_incr, + const std::string& coupling_type, double* params, + double* solution, double& cpl_value, char last_flag, + int& error_code); + + /// @brief Retrieve the index within the solution vector that corresponds to + /// the coupled boundary DOF. + void extract_coupled_dof(int problem_id, int& coupled_dof, + const std::string& coupling_type); + + + private: + /// @brief Problem identifier assigned by the 1D solver during initialization. + int problem_id_ = 0; + /// @brief Total number of DOFs (nodes * 2: flow + area) assigned by the 1D solver during initialization. + int system_size_ = 0; + + void* library_handle_ = nullptr; + + // Function pointers to shared-library symbols. + + /// @brief Initialize the 1D solver + void (*initialize_1d_)(const char*, int&, int&, const char*) = nullptr; + + /// @brief Set the 1D solver's internal time step size + void (*set_external_step_size_1d_)(int, double) = nullptr; + + /// @brief Return the current 1D solution vector + void (*return_1d_solution_)(int, double*, int) = nullptr; + + /// @brief Update the 1D solver's current solution vector + void (*update_1d_solution_)(int, double*, int) = nullptr; + + /// @brief Advance the 1D solver by one 3D time step + void (*run_1d_simulation_step_1d_)(int, double, int, const char*, double*, + double*, double&, char*, int&) = nullptr; + + /// @brief Retrieve the index within the solution vector that corresponds to + /// the coupled boundary DOF. + void (*extract_coupled_dof_)(int, int&, char*) = nullptr; +}; + +#endif // ONEDSOLVER_INTERFACE_H diff --git a/Code/Source/solver/svZeroD_interface.cpp b/Code/Source/solver/svZeroD_interface.cpp index 470f5fff8..2a3c8da8f 100644 --- a/Code/Source/solver/svZeroD_interface.cpp +++ b/Code/Source/solver/svZeroD_interface.cpp @@ -27,7 +27,13 @@ static void build_svzero_coupled_bc_idxs(ComMod& com_mod) cpl.svZeroD_coupled_bc_idxs.clear(); for (int iEq = 0; iEq < com_mod.nEq; iEq++) { for (int iBc = 0; iBc < com_mod.eq[iEq].nBc; iBc++) { - if (utils::btest(com_mod.eq[iEq].bc[iBc].bType, consts::iBC_Coupled)) { + const auto& bc = com_mod.eq[iEq].bc[iBc]; + // A Coupled face belongs to svZeroD only when it is not a 1D face. + // The is_svOneD_face() helper encapsulates the routing invariant: + // 1D face → oned_input_file_ non-empty, block_name_ empty + // 0D face → block_name_ non-empty, oned_input_file_ empty + if (utils::btest(bc.bType, consts::iBC_Coupled) && + !bc.coupled_bc.is_svOneD_face()) { cpl.svZeroD_coupled_bc_idxs.emplace_back(iEq, iBc); } } @@ -386,6 +392,14 @@ void calc_svZeroD(ComMod& com_mod, const CmMod& cm_mod, char BCFlag) double params[2]; double times[2]; int error_code; + + // Temporary arrays to record relaxed values for history update on 'L' steps. + std::vector P_sent_old_arr(numCoupledSrfs, 0.0); + std::vector P_sent_new_arr(numCoupledSrfs, 0.0); + std::vector Q_input_sent_old_arr(numCoupledSrfs, 0.0); + std::vector Q_input_sent_new_arr(numCoupledSrfs, 0.0); + std::vector Q_relaxed_arr(numCoupledSrfs, 0.0); + std::vector P_relaxed_arr(numCoupledSrfs, 0.0); get_coupled_QP(com_mod, QCoupled, QnCoupled, PCoupled, PnCoupled); @@ -416,11 +430,53 @@ void calc_svZeroD(ComMod& com_mod, const CmMod& cm_mod, char BCFlag) double sign = bc->coupled_bc.get_in_out_sign(); if (is_dirichlet) { - params[0] = PCoupled[i]; - params[1] = PnCoupled[i]; + double raw_P_old = PCoupled[i]; + double raw_P_new = PnCoupled[i]; + + // Step 1: apply pressure ramp (scales amplitude from ramp_ref_pressure + // to actual 3D pressure over the first ramp_steps steps). + int ramp_steps = bc->coupled_bc.get_oned_ramp_steps(); + double P_target_old, P_target_new; + if (ramp_steps > 0) { + double ramp_factor = std::min(1.0, static_cast(bc->coupled_bc.get_ramp_step_count()) / ramp_steps); + double P_ref = bc->coupled_bc.get_oned_ramp_ref_pressure(); + P_target_old = P_ref + ramp_factor * (raw_P_old - P_ref); + P_target_new = P_ref + ramp_factor * (raw_P_new - P_ref); + } else { + P_target_old = raw_P_old; + P_target_new = raw_P_new; + } + + // Step 2: apply under-relaxation (damps timestep-to-timestep oscillations). + // P_sent = omega * P_target + (1 - omega) * P_prev_sent + const double omega = bc->coupled_bc.get_oned_relax_factor(); + params[0] = omega * P_target_old + (1.0 - omega) * bc->coupled_bc.get_P_prev_sent_old(); + params[1] = omega * P_target_new + (1.0 - omega) * bc->coupled_bc.get_P_prev_sent_new(); + P_sent_old_arr[i] = params[0]; + P_sent_new_arr[i] = params[1]; } else { - params[0] = sign * QCoupled[i]; - params[1] = sign * QnCoupled[i]; + double raw_Q_old = sign * QCoupled[i]; + double raw_Q_new = sign * QnCoupled[i]; + + // Step 1: apply flow ramp (scales Q from ramp_ref over the first ramp_steps steps). + int ramp_steps = bc->coupled_bc.get_oned_ramp_steps(); + double Q_target_old, Q_target_new; + if (ramp_steps > 0) { + double ramp_factor = std::min(1.0, static_cast(bc->coupled_bc.get_ramp_step_count()) / ramp_steps); + double Q_ref = 0.0; // ref value (0.0 by default) + Q_target_old = Q_ref + ramp_factor * (raw_Q_old - Q_ref); + Q_target_new = Q_ref + ramp_factor * (raw_Q_new - Q_ref); + } else { + Q_target_old = raw_Q_old; + Q_target_new = raw_Q_new; + } + + // Step 2: apply under-relaxation. + const double omega = bc->coupled_bc.get_oned_relax_factor(); + params[0] = omega * Q_target_old + (1.0 - omega) * bc->coupled_bc.get_Q_input_prev_old(); + params[1] = omega * Q_target_new + (1.0 - omega) * bc->coupled_bc.get_Q_input_prev_new(); + Q_input_sent_old_arr[i] = params[0]; + Q_input_sent_new_arr[i] = params[1]; } update_svZeroD_block_params(svzd_blk_names[i], times, params); } @@ -447,12 +503,20 @@ void calc_svZeroD(ComMod& com_mod, const CmMod& cm_mod, char BCFlag) } if (bc->coupled_bc.get_bc_type() == consts::BoundaryConditionType::bType_Neu) { - PCoupled[i] = lpn_state_y[pressure_id]; - bc->coupled_bc.set_pressure(PCoupled[i]); + double P_raw = lpn_state_y[pressure_id]; + // Apply under-relaxation to the pressure output. + const double omega = bc->coupled_bc.get_oned_relax_factor(); + double P_relaxed = omega * P_raw + (1.0 - omega) * bc->coupled_bc.get_P_neu_prev(); + bc->coupled_bc.set_pressure(P_relaxed); + P_relaxed_arr[i] = P_relaxed; } else if (bc->coupled_bc.get_bc_type() == consts::BoundaryConditionType::bType_Dir) { - QCoupled[i] = in_out * lpn_state_y[flow_id]; + double Q_raw = in_out * lpn_state_y[flow_id]; + // Apply under-relaxation to the flow output to damp timestep-to-timestep oscillations. + const double omega = bc->coupled_bc.get_oned_relax_factor(); + double Q_relaxed = omega * Q_raw + (1.0 - omega) * bc->coupled_bc.get_Q_prev_sent(); double Qo_prev = bc->coupled_bc.get_Qn(); - bc->coupled_bc.set_flowrates(Qo_prev, QCoupled[i]); + bc->coupled_bc.set_flowrates(Qo_prev, Q_relaxed); + Q_relaxed_arr[i] = Q_relaxed; } else { throw std::runtime_error("ERROR: [calc_svZeroD] Invalid Coupled BC type."); } @@ -463,6 +527,20 @@ void calc_svZeroD(ComMod& com_mod, const CmMod& cm_mod, char BCFlag) interface->return_ydot(last_state_ydot); std::copy(lpn_state_y.begin(), lpn_state_y.end(), last_state_y.begin()); + // Update ramp/relax history for all coupled BCs. + for (int i = 0; i < numCoupledSrfs; ++i) { + bcType* bc_hist = nullptr; + if (!nth_coupled_bc(com_mod, i, &bc_hist)) continue; + if (bc_hist->coupled_bc.get_bc_type() == consts::BoundaryConditionType::bType_Dir) { + bc_hist->coupled_bc.set_P_prev_sent(P_sent_old_arr[i], P_sent_new_arr[i]); + bc_hist->coupled_bc.set_Q_prev_sent(Q_relaxed_arr[i]); + } else { + bc_hist->coupled_bc.set_Q_input_prev(Q_input_sent_old_arr[i], Q_input_sent_new_arr[i]); + bc_hist->coupled_bc.set_P_neu_prev(P_relaxed_arr[i]); + } + bc_hist->coupled_bc.increment_ramp_step_count(); + } + if (writeSvZeroD == 1) { // Write the state vector to a file int arg = 1; @@ -486,6 +564,9 @@ void calc_svZeroD(ComMod& com_mod, const CmMod& cm_mod, char BCFlag) if (utils::btest(bc.bType, iBC_Coupled) && bc.coupled_bc.get_bc_type() == BoundaryConditionType::bType_Neu) { bc.coupled_bc.bcast_coupled_neumann_pressure(cm_mod, cm); + } else if (utils::btest(bc.bType, iBC_Coupled) && + bc.coupled_bc.get_bc_type() == BoundaryConditionType::bType_Dir) { + bc.coupled_bc.bcast_coupled_dir_flowrate(cm_mod, cm); } } } diff --git a/Code/Source/solver/txt.cpp b/Code/Source/solver/txt.cpp index 61e3060a0..89f03c907 100644 --- a/Code/Source/solver/txt.cpp +++ b/Code/Source/solver/txt.cpp @@ -13,6 +13,7 @@ #include "utils.h" #include #include "svZeroD_interface.h" +#include "svOneD_interface.h" namespace txt_ns { @@ -174,16 +175,52 @@ void txt(Simulation* simulation, const bool init_write, const SolutionStates& so if (!init_write) { if (cplBC.useGenBC) { set_bc::genBC_Integ_X(com_mod, cm_mod, "L"); - } else if (cplBC.useSvZeroD) { - svZeroD::calc_svZeroD(com_mod, cm_mod, 'L'); } else { - for (auto& bc : com_mod.eq[0].bc) { - if (utils::btest(bc.bType, iBC_RCR)) { - ltmp = true; - break; + // In mixed-coupling simulations both useSvZeroD and useSvOneD can be true. + // Call each active solver independently. + if (cplBC.useSvZeroD) { + svZeroD::calc_svZeroD(com_mod, cm_mod, 'L'); + } + + if (cplBC.useSvOneD) { + // Update NEU coupling flowrates from the final converged velocity + // field (Yn) before committing the 1D solution. During the Newton + // loop, compute_flowrates() is called at the *start* of each + // iteration using a partially-converged Yn. After the last picc + // correction Yn is fully converged, but Qn_ is never refreshed. + // Re-computing here ensures the 1D solver receives the same + // flowrate that write_boundary_integral_data() reports in + // B_NS_Velocity_flux.txt. + for (auto& bc : com_mod.eq[0].bc) { + if (utils::btest(bc.bType, iBC_Coupled) && + bc.coupled_bc.is_svOneD_face() && + bc.coupled_bc.get_bc_type() == BoundaryConditionType::bType_Neu) { + bc.coupled_bc.compute_flowrates(com_mod, cm_mod, solutions); + } + } + svOneD::calc_svOneD(com_mod, cm_mod, 'L'); + } + + // Also integrate any RCR faces coexisting with svZeroD/svOneD faces. + if (cplBC.useSvZeroD || cplBC.useSvOneD) { + for (auto& bc : com_mod.eq[0].bc) { + if (utils::btest(bc.bType, iBC_RCR)) { + ltmp = true; + break; + } + } + if (ltmp) { + set_bc::cplBC_Integ_X(com_mod, cm_mod, true); + } + } else { + for (auto& bc : com_mod.eq[0].bc) { + if (utils::btest(bc.bType, iBC_RCR)) { + ltmp = true; + break; + } } + set_bc::cplBC_Integ_X(com_mod, cm_mod, ltmp); } - set_bc::cplBC_Integ_X(com_mod, cm_mod, ltmp); } } } diff --git a/tests/cases/fluid/Tpipe_svOneD_2faces/inlet_1d.in b/tests/cases/fluid/Tpipe_svOneD_2faces/inlet_1d.in new file mode 100644 index 000000000..2bcf832c0 --- /dev/null +++ b/tests/cases/fluid/Tpipe_svOneD_2faces/inlet_1d.in @@ -0,0 +1,125 @@ +# ================================ +# inlet_1d MODEL - UNITS IN CGS +# ================================ + +# ========== +# MODEL CARD +# ========== +# - Name of the model (string) + +MODEL inlet_1d + + + +### DO NOT CHANGE THIS SECTION - generated automatically +# +# ========== +# NODE CARD +# ========== +# - Node Name (double) +# - Node X Coordinate (double) +# - Node Y Coordinate (double) +# - Node Z Coordinate (double) + +NODE 0 -6.320068093537426e-18 5.27257554949756e-07 -10.0 +NODE 1 2.3869191998230702e-18 5.27257554949756e-07 0.0 + + +### DO NOT CHANGE THIS SECTION - generated automatically +# +# ========== +# JOINT CARD +# ========== +# - Joint Name (string) +# - Joint Node (double) +# - Joint Inlet Name (string) +# - Joint Outlet Name (string) + + + +### DO NOT CHANGE THIS SECTION - generated automatically +# +# ================================ +# JOINTINLET AND JOINTOUTLET CARDS +# ================================ +# - Inlet/Outlet Name (string) +# - Total Number of segments (int) +# - List of segments (list of int) + +# ============ +# SEGMENT CARD +# ============ +# - Segment Name (string) +# - Segment ID (int) +# - Segment Length (double) +# - Total Finite Elements in Segment (int) +# - Segment Inlet Node (int) +# - Segment Outlet Node (int) +# - Segment Inlet Area (double) +# - Segment Outlet Area (double) +# - Segment Inflow Value (double) +# - Segment Material (string) +# - Type of Loss (string - 'NONE','STENOSIS','BRANCH_THROUGH_DIVIDING','BRANCH_SIDE_DIVIDING','BRANCH_THROUGH_CONVERGING', +# 'BRANCH_SIDE_CONVERGING','BIFURCATION_BRANCH') +# - Branch Angle (double) +# - Upstream Segment ID (int) +# - Branch Segment ID (int) +# - Boundary Condition Type (string - 'NOBOUND','PRESSURE','AREA','FLOW','RESISTANCE','RESISTANCE_TIME','PRESSURE_WAVE', +# 'WAVE','RCR','CORONARY','IMPEDANCE','PULMONARY') +# - Data Table Name (string) + +SEGMENT branch0_seg0 0 10.0 20 0 1 3.130298366662493 3.1302983666630864 0.0 MAT1 NONE 0.0 0 0 COUPLED NONE + + + +DATATABLE INFLOW LIST +0.0 100.0 +10.0 100.0 +ENDDATATABLE + + +# ================== +# SOLVEROPTIONS CARD +# ================== +# - Solver Time Step (double), +# - Steps Between Saves (int), +# - Max Number of Steps (int) +# - Number of quadrature points for finite elements (int), +# - Name of Datatable for inlet conditions (string) +# - Type of boundary condition (string - 'NOBOUND','PRESSURE','AREA','FLOW','RESISTANCE','RESISTANCE_TIME','PRESSURE_WAVE', +# 'WAVE','RCR','CORONARY','IMPEDANCE','PULMONARY') +# - Convergence tolerance (double), +# - Formulation Type (int - 0 Advective, 1 Conservative), +# - Stabilization (int - 0 No stabilization, 1 With stabilization) + +SOLVEROPTIONS 0.005 20 100 2 INFLOW FLOW 1.0e-5 1 1 + +COUPLING ON DIR 1 + +# ============= +# MATERIAL CARD +# ============= +# - Material Name (string) +# - Material Type (string - 'LINEAR','OLUFSEN') +# - Material Density (double) +# - Material Viscosity (double) +# - Material PRef (double) +# - Material Exponent (double) +# - Material Parameter 1 (double) +# - Material Parameter 2 (double) +# - Material Parameter 3 (double) + +MATERIAL MAT1 LINEAR 1.06 0.04 0.0 1.0 10000000.0 0.0 0.0 + +# ============ +# OUTPUT CARD +# ============ +# +# 1. Output file format. The following output types are supported: +# TEXT. The output of every segment is written in separate text files for the flow rate, pressure, area and Reynolds number. The rows contain output values at varying locations along the segment while columns contains results at various time instants. +# VTK. The results for all time steps are plotted to a 3D-like model using the XML VTK file format. +# 2. VTK export option. Two options are available for VTK file outputs: +# 0 - Multiple files (default). A separate file is written for each saved increment. A pvd file is also provided which contains the time information of the sequence. This is the best option to create animations. +# 1 - The results for all time steps are plotted to a single XML VTK file. + +OUTPUT TEXT diff --git a/tests/cases/fluid/Tpipe_svOneD_2faces/inlet_1d_00050.vtp b/tests/cases/fluid/Tpipe_svOneD_2faces/inlet_1d_00050.vtp new file mode 100644 index 000000000..6114f8b75 --- /dev/null +++ b/tests/cases/fluid/Tpipe_svOneD_2faces/inlet_1d_00050.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ad8a9596cd4f77a927c217f02e2f540fe82e7a20927f1bb1457ff20d11fda578 +size 66000 diff --git a/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-complete.exterior.vtp b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-complete.exterior.vtp new file mode 100644 index 000000000..324c6ec2b --- /dev/null +++ b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-complete.exterior.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b6d3ed0de2bbb08a7dc17388e1de4d7504096fbe5abe8d0bf88cab534353efda +size 119781 diff --git a/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-complete.mesh.vtu b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-complete.mesh.vtu new file mode 100644 index 000000000..58e35a91a --- /dev/null +++ b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1ba04aa43f6fcc96feb1c349840de4acac0329b6d3f388bf9ec17bc9ff8018ff +size 642824 diff --git a/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-surfaces/cap_main_3d.vtp b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-surfaces/cap_main_3d.vtp new file mode 100644 index 000000000..b98c01a7f --- /dev/null +++ b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-surfaces/cap_main_3d.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:732c605b92af8941d9cce53caf89928fb9de1876498539ad6be1ae6b4c814d69 +size 4910 diff --git a/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-surfaces/cap_main_3d_2.vtp b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-surfaces/cap_main_3d_2.vtp new file mode 100644 index 000000000..4efbbd59b --- /dev/null +++ b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-surfaces/cap_main_3d_2.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5724183e2c702cb93b14c8d35b10cccfbada94c2a54556cf35e3524faac76f35 +size 4828 diff --git a/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-surfaces/cap_side_3d.vtp b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-surfaces/cap_side_3d.vtp new file mode 100644 index 000000000..ee0e3e1cb --- /dev/null +++ b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-surfaces/cap_side_3d.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5bfc9a71ddd1b378205bd551f932db3da5f9945590aff9f80d28c294dd18c0f3 +size 3511 diff --git a/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-surfaces/wall_main_3d.vtp b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-surfaces/wall_main_3d.vtp new file mode 100644 index 000000000..5e0c767f5 --- /dev/null +++ b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-surfaces/wall_main_3d.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:be46a91f418bd88cf8711de48e158b16937c9487278f6d4eafad395d30f68c35 +size 95467 diff --git a/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-surfaces/wall_side_3d.vtp b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-surfaces/wall_side_3d.vtp new file mode 100644 index 000000000..7a2f94350 --- /dev/null +++ b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/mesh-surfaces/wall_side_3d.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:56f42230ac103ace09256b6c49fa5213c1b7945d155565358816f50a17800163 +size 20913 diff --git a/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/walls_combined.vtp b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/walls_combined.vtp new file mode 100644 index 000000000..08ef317b4 --- /dev/null +++ b/tests/cases/fluid/Tpipe_svOneD_2faces/mesh-complete/walls_combined.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4c049bf671ca4bbc51054e32707305a879cf50fbdf5514d5de5aedadcda34c52 +size 113058 diff --git a/tests/cases/fluid/Tpipe_svOneD_2faces/outlet_1d.in b/tests/cases/fluid/Tpipe_svOneD_2faces/outlet_1d.in new file mode 100644 index 000000000..1e6e03e71 --- /dev/null +++ b/tests/cases/fluid/Tpipe_svOneD_2faces/outlet_1d.in @@ -0,0 +1,137 @@ +# ================================ +# main_1d MODEL - UNITS IN CGS +# ================================ + +# ========== +# MODEL CARD +# ========== +# - Name of the model (string) + +MODEL outlet_1d + + + +### DO NOT CHANGE THIS SECTION - generated automatically +# +# ========== +# NODE CARD +# ========== +# - Node Name (double) +# - Node X Coordinate (double) +# - Node Y Coordinate (double) +# - Node Z Coordinate (double) + +NODE 0 -6.320068093537426e-18 5.27257554949756e-07 22.520139694213867 +NODE 1 -6.72370548127077e-10 2.772360858216416e-05 27.36005401611328 +NODE 2 -0.0007925734389573336 1.7465894416091032e-05 28.852317810058594 +NODE 3 2.3869191998230702e-18 5.27257554949756e-07 30.0 + + +### DO NOT CHANGE THIS SECTION - generated automatically +# +# ========== +# JOINT CARD +# ========== +# - Joint Name (string) +# - Joint Node (double) +# - Joint Inlet Name (string) +# - Joint Outlet Name (string) + + + +### DO NOT CHANGE THIS SECTION - generated automatically +# +# ================================ +# JOINTINLET AND JOINTOUTLET CARDS +# ================================ +# - Inlet/Outlet Name (string) +# - Total Number of segments (int) +# - List of segments (list of int) + +JOINT J0 1 IN0 OUT0 +JOINTINLET IN0 1 0 +JOINTOUTLET OUT0 1 1 + +JOINT J1 2 IN1 OUT1 +JOINTINLET IN1 1 1 +JOINTOUTLET OUT1 1 2 + +# ============ +# SEGMENT CARD +# ============ +# - Segment Name (string) +# - Segment ID (int) +# - Segment Length (double) +# - Total Finite Elements in Segment (int) +# - Segment Inlet Node (int) +# - Segment Outlet Node (int) +# - Segment Inlet Area (double) +# - Segment Outlet Area (double) +# - Segment Inflow Value (double) +# - Segment Material (string) +# - Type of Loss (string - 'NONE','STENOSIS','BRANCH_THROUGH_DIVIDING','BRANCH_SIDE_DIVIDING','BRANCH_THROUGH_CONVERGING', +# 'BRANCH_SIDE_CONVERGING','BIFURCATION_BRANCH') +# - Branch Angle (double) +# - Upstream Segment ID (int) +# - Branch Segment ID (int) +# - Boundary Condition Type (string - 'NOBOUND','PRESSURE','AREA','FLOW','RESISTANCE','RESISTANCE_TIME','PRESSURE_WAVE', +# 'WAVE','RCR','CORONARY','IMPEDANCE','PULMONARY') +# - Data Table Name (string) + +SEGMENT branch0_seg0 0 4.8399151774568026 5 0 1 3.130298369053221 3.1302983615097206 0.0 MAT1 NONE 0.0 0 0 NOBOUND NONE +SEGMENT branch0_seg1 1 1.49226421600263 5 1 2 3.130298345279223 3.130298345279223 0.0 MAT1 NONE 0.0 0 0 NOBOUND NONE +SEGMENT branch0_seg2 2 1.1476824649695037 5 2 3 3.13029837283059 3.1302983690532207 0.0 MAT1 NONE 0.0 0 0 RESISTANCE RESISTANCE_0 + + +DATATABLE RESISTANCE_0 LIST +0.0 100 +0.0 0.0 +ENDDATATABLE + + + +# ================== +# SOLVEROPTIONS CARD +# ================== +# - Solver Time Step (double), +# - Steps Between Saves (int), +# - Max Number of Steps (int) +# - Number of quadrature points for finite elements (int), +# - Name of Datatable for inlet conditions (string) +# - Type of boundary condition (string - 'NOBOUND','PRESSURE','AREA','FLOW','RESISTANCE','RESISTANCE_TIME','PRESSURE_WAVE', +# 'WAVE','RCR','CORONARY','IMPEDANCE','PULMONARY') +# - Convergence tolerance (double), +# - Formulation Type (int - 0 Advective, 1 Conservative), +# - Stabilization (int - 0 No stabilization, 1 With stabilization) + +SOLVEROPTIONS 0.005 20 100 2 NONE COUPLED 1.0e-5 1 1 + +COUPLING ON NEU 1 + +# ============= +# MATERIAL CARD +# ============= +# - Material Name (string) +# - Material Type (string - 'LINEAR','OLUFSEN') +# - Material Density (double) +# - Material Viscosity (double) +# - Material PRef (double) +# - Material Exponent (double) +# - Material Parameter 1 (double) +# - Material Parameter 2 (double) +# - Material Parameter 3 (double) + +MATERIAL MAT1 LINEAR 1.06 0.04 0.0 1.0 10000000.0 0.0 0.0 + +# ============ +# OUTPUT CARD +# ============ +# +# 1. Output file format. The following output types are supported: +# TEXT. The output of every segment is written in separate text files for the flow rate, pressure, area and Reynolds number. The rows contain output values at varying locations along the segment while columns contains results at various time instants. +# VTK. The results for all time steps are plotted to a 3D-like model using the XML VTK file format. +# 2. VTK export option. Two options are available for VTK file outputs: +# 0 - Multiple files (default). A separate file is written for each saved increment. A pvd file is also provided which contains the time information of the sequence. This is the best option to create animations. +# 1 - The results for all time steps are plotted to a single XML VTK file. + +OUTPUT VTK 0 diff --git a/tests/cases/fluid/Tpipe_svOneD_2faces/outlet_1d_00050.vtp b/tests/cases/fluid/Tpipe_svOneD_2faces/outlet_1d_00050.vtp new file mode 100644 index 000000000..f7573ec53 --- /dev/null +++ b/tests/cases/fluid/Tpipe_svOneD_2faces/outlet_1d_00050.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0528a3be307a2e54b02f17b763c3ddd98064b2336ea1bd79a80e542e343923b3 +size 57578 diff --git a/tests/cases/fluid/Tpipe_svOneD_2faces/result_050.vtu b/tests/cases/fluid/Tpipe_svOneD_2faces/result_050.vtu new file mode 100644 index 000000000..a15ed0974 --- /dev/null +++ b/tests/cases/fluid/Tpipe_svOneD_2faces/result_050.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6fb6056a48016488c81f4005c78f08ef506a8af453341e247b1cef9de9f1a838 +size 1768310 diff --git a/tests/cases/fluid/Tpipe_svOneD_2faces/solver.xml b/tests/cases/fluid/Tpipe_svOneD_2faces/solver.xml new file mode 100644 index 000000000..30afe711b --- /dev/null +++ b/tests/cases/fluid/Tpipe_svOneD_2faces/solver.xml @@ -0,0 +1,112 @@ + + + false + 3 + 500 + 0.005 + 0.5 + 10000 + 1 + true + result + 50 + 0.5 + + + mesh-complete/mesh-complete.mesh.vtu + + mesh-complete/mesh-surfaces/cap_main_3d.vtp + + + mesh-complete/mesh-surfaces/cap_main_3d_2.vtp + + + mesh-complete/mesh-surfaces/cap_side_3d.vtp + + + mesh-complete/mesh-surfaces/wall_main_3d.vtp + + + mesh-complete/mesh-surfaces/wall_side_3d.vtp + + + + true + 3 + 20 + 1e-5 + 0.2 + 1.06 + + 0.04 + + + + fsils + + 15 + 1e-4 + 200 + 10 + 1e-3 + 300 + 1e-3 + + + + semi-implicit + /Users/taeoukkim/Desktop/Simvascular_solvers/1Dsolver/svOneDSolver/build/lib/libsvoned_interface.dylib + + + + + true + true + true + true + true + true + + + true + true + true + + + Dir + Coupled + + inlet_1d.in + 100 + 0.0 + 0.3 + + Parabolic + + + Neu + Coupled + + outlet_1d.in + 100 + 0.0 + 0.3 + + + + Neu + Resistance + 100.0 + + + Dirichlet + Steady + 0 + + + Dirichlet + Steady + 0 + + + diff --git a/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-complete.exterior.vtp b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-complete.exterior.vtp new file mode 100644 index 000000000..324c6ec2b --- /dev/null +++ b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-complete.exterior.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b6d3ed0de2bbb08a7dc17388e1de4d7504096fbe5abe8d0bf88cab534353efda +size 119781 diff --git a/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-complete.mesh.vtu b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-complete.mesh.vtu new file mode 100644 index 000000000..58e35a91a --- /dev/null +++ b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1ba04aa43f6fcc96feb1c349840de4acac0329b6d3f388bf9ec17bc9ff8018ff +size 642824 diff --git a/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-surfaces/cap_main_3d.vtp b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-surfaces/cap_main_3d.vtp new file mode 100644 index 000000000..b98c01a7f --- /dev/null +++ b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-surfaces/cap_main_3d.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:732c605b92af8941d9cce53caf89928fb9de1876498539ad6be1ae6b4c814d69 +size 4910 diff --git a/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-surfaces/cap_main_3d_2.vtp b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-surfaces/cap_main_3d_2.vtp new file mode 100644 index 000000000..4efbbd59b --- /dev/null +++ b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-surfaces/cap_main_3d_2.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5724183e2c702cb93b14c8d35b10cccfbada94c2a54556cf35e3524faac76f35 +size 4828 diff --git a/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-surfaces/cap_side_3d.vtp b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-surfaces/cap_side_3d.vtp new file mode 100644 index 000000000..ee0e3e1cb --- /dev/null +++ b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-surfaces/cap_side_3d.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5bfc9a71ddd1b378205bd551f932db3da5f9945590aff9f80d28c294dd18c0f3 +size 3511 diff --git a/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-surfaces/wall_main_3d.vtp b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-surfaces/wall_main_3d.vtp new file mode 100644 index 000000000..5e0c767f5 --- /dev/null +++ b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-surfaces/wall_main_3d.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:be46a91f418bd88cf8711de48e158b16937c9487278f6d4eafad395d30f68c35 +size 95467 diff --git a/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-surfaces/wall_side_3d.vtp b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-surfaces/wall_side_3d.vtp new file mode 100644 index 000000000..7a2f94350 --- /dev/null +++ b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/mesh-surfaces/wall_side_3d.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:56f42230ac103ace09256b6c49fa5213c1b7945d155565358816f50a17800163 +size 20913 diff --git a/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/walls_combined.vtp b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/walls_combined.vtp new file mode 100644 index 000000000..08ef317b4 --- /dev/null +++ b/tests/cases/fluid/Tpipe_svZeroD_Dir/mesh-complete/walls_combined.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4c049bf671ca4bbc51054e32707305a879cf50fbdf5514d5de5aedadcda34c52 +size 113058 diff --git a/tests/cases/fluid/Tpipe_svZeroD_Dir/result_050.vtu b/tests/cases/fluid/Tpipe_svZeroD_Dir/result_050.vtu new file mode 100644 index 000000000..156caee26 --- /dev/null +++ b/tests/cases/fluid/Tpipe_svZeroD_Dir/result_050.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:649dc557eb87737a69bbd9cbd1fdbd9ac95f1a3c7b0d37d409938b0ff2489be2 +size 1758470 diff --git a/tests/cases/fluid/Tpipe_svZeroD_Dir/solver.xml b/tests/cases/fluid/Tpipe_svZeroD_Dir/solver.xml new file mode 100644 index 000000000..c4499cf8f --- /dev/null +++ b/tests/cases/fluid/Tpipe_svZeroD_Dir/solver.xml @@ -0,0 +1,110 @@ + + + false + 3 + 300 + 0.001 + 0.5 + 10000 + 1 + true + result + 50 + 0.5 + + + mesh-complete/mesh-complete.mesh.vtu + + mesh-complete/mesh-surfaces/cap_main_3d.vtp + + + mesh-complete/mesh-surfaces/cap_main_3d_2.vtp + + + mesh-complete/mesh-surfaces/cap_side_3d.vtp + + + mesh-complete/mesh-surfaces/wall_main_3d.vtp + + + mesh-complete/mesh-surfaces/wall_side_3d.vtp + + + + true + 3 + 20 + 1e-5 + 0.2 + 1.06 + + 0.04 + + + + fsils + + 15 + 1e-4 + 200 + 10 + 1e-3 + 300 + 1e-3 + + + + semi-implicit + solver_0d.json + /Users/taeoukkim/Desktop/Simvascular_solvers/0Dsolver/svZeroDSolver/build/src/interface/libsvzero_interface.dylib + 0.0 + 0.0 + + + + + true + true + true + true + true + true + + + true + true + true + + + Dir + Coupled + + cp1 + 100 + 0.0 + 0.2 + + Parabolic + + + Neu + Resistance + 100.0 + + + Neu + Resistance + 100.0 + + + Dirichlet + Steady + 0 + + + Dirichlet + Steady + 0 + + + diff --git a/tests/cases/fluid/Tpipe_svZeroD_Dir/solver_0d.json b/tests/cases/fluid/Tpipe_svZeroD_Dir/solver_0d.json new file mode 100644 index 000000000..5a3f74346 --- /dev/null +++ b/tests/cases/fluid/Tpipe_svZeroD_Dir/solver_0d.json @@ -0,0 +1,52 @@ +{ + "boundary_conditions": [ + { + "bc_name": "INFLOW", + "bc_type": "FLOW", + "bc_values": { + "Q": [ + 100.0, + 100.0 + ], + "t": [ + 0.0, + 1.0 + ] + } + } + ], + "external_solver_coupling_blocks": [ + { + "name": "cp1", + "type": "PRESSURE", + "location": "outlet", + "connected_block": "vessel1", + "periodic": false, + "values": { + "t": [0.0, 1.0], + "P": [1.0, 1.0] + } + } + ], + "junctions": [], + "simulation_parameters": { + "coupled_simulation": true, + "number_of_time_pts": 10, + "output_all_cycles": true, + "steady_initial": false + }, + "vessels": [ + { + "boundary_conditions": { + "inlet": "INFLOW" + }, + "vessel_id": 0, + "vessel_length": 1.0, + "vessel_name": "vessel1", + "zero_d_element_type": "BloodVessel", + "zero_d_element_values": { + "R_poiseuille": 100.0 + } + } + ] +} \ No newline at end of file