diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index c546c2822..05a3f1e9c 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_subroutines.h svOneD_subroutines.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..7d543dd15 100644 --- a/Code/Source/solver/ComMod.cpp +++ b/Code/Source/solver/ComMod.cpp @@ -217,3 +217,13 @@ void svZeroDSolverInterfaceType::set_data(const svZeroDSolverInterfaceParameters has_data = true; } + +void svOneDSolverInterfaceType::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 8b41e5c6a..81ca89a53 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,6 +775,12 @@ 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; }; //---------------------------- @@ -807,6 +817,26 @@ class svZeroDSolverInterfaceType void set_data(const svZeroDSolverInterfaceParameters& params); }; +//---------------------------- +// svOneDSolverInterfaceType +//---------------------------- +// This class stores information used to interface to the svOneDSolver. +// +class svOneDSolverInterfaceType +{ + public: + // Path to the 1D solver shared library (without .so/.dylib extension, + // or with extension if the full path is provided). + std::string solver_library; + + // If the data has been set for the interface. Set to true after + // the svOneDSolver_interface XML element has been parsed. + // Note: the per-face input files are stored in cplFaceType::oned_input_file. + bool has_data = false; + + void set_data(const svOneDSolverInterfaceParameters& params); +}; + /// @brief For coupled 0D-3D problems // class cplBCType @@ -822,6 +852,9 @@ class cplBCType // Whether to use svZeroD bool useSvZeroD = false; + // Whether to use svOneD (svOneDSolver) + bool useSv1D = false; + // Whether to initialize RCR from flow data bool initRCR = false; @@ -853,6 +886,8 @@ class cplBCType svZeroDSolverInterfaceType svzerod_solver_interface; + svOneDSolverInterfaceType sv1d_solver_interface; + /// @brief The name of history file containing "X" std::string saveName; //std::string(LEN=stdL) :: saveName = "LPN.dat"; @@ -1633,6 +1668,9 @@ class ComMod { bool usePrecomp = false; //----- int members -----// + /// @brief scalar initial pressure is explicitly provided in the input file + bool have_initial_pressure = false; + /// @brief Current domain int cDmn = 0; @@ -1715,6 +1753,9 @@ class ComMod { /// @brief Time double time = 0.0; + /// @brief initial scalar pressure value + double initial_pressure = 0.0; + //----- string members -----// diff --git a/Code/Source/solver/CoupledBoundaryCondition.cpp b/Code/Source/solver/CoupledBoundaryCondition.cpp index 5f7602e3e..56f917afa 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_); @@ -875,8 +937,27 @@ CappingSurface::CappingSurface(const CappingSurface& other) , valM_(other.valM_) , normals_(other.normals_) { - if (other.face_) { - face_ = std::make_unique(*other.face_); + if (other.face_ != nullptr) { + try { + face_ = std::make_unique(*other.face_); + } catch (const std::exception& e) { + throw CappingSurfaceCopyException("[CappingSurface::copy constructor] Failed to copy face_: " + + std::string(e.what())); + } + if (face_ != nullptr) { + if (face_->nNo > 0 && face_->gN.size() != face_->nNo) { + throw CappingSurfaceCopyException("[CappingSurface::copy constructor] Invalid face_: gN.size()=" + + std::to_string(face_->gN.size()) + " != nNo=" + + std::to_string(face_->nNo)); + } + if (face_->nEl > 0 && face_->IEN.ncols() != face_->nEl) { + throw CappingSurfaceCopyException("[CappingSurface::copy constructor] Invalid face_: IEN.ncols()=" + + std::to_string(face_->IEN.ncols()) + " != nEl=" + + std::to_string(face_->nEl)); + } + } + } else { + face_.reset(); } } @@ -886,8 +967,13 @@ CappingSurface& CappingSurface::operator=(const CappingSurface& other) global_node_ids_ = other.global_node_ids_; valM_ = other.valM_; normals_ = other.normals_; - if (other.face_) { - face_ = std::make_unique(*other.face_); + if (other.face_ != nullptr) { + try { + face_ = std::make_unique(*other.face_); + } catch (const std::exception& e) { + throw CappingSurfaceCopyException("[CappingSurface::operator=] Failed to copy face_: " + + std::string(e.what())); + } } else { face_.reset(); } @@ -986,7 +1072,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; @@ -1034,6 +1120,10 @@ Array CappingSurface::update_element_position_global(int e, consts::Mech { using namespace consts; + if (mesh_x.nrows() < cap_nsd_ || mesh_Do.nrows() < cap_nsd_ || mesh_Dn.nrows() < cap_nsd_) { + throw CappingSurfaceGeometryException( + "[CappingSurface::update_element_position_global] Mesh arrays must have at least 3 rows."); + } Array xl(cap_nsd_, face_->eNoN); for (int a = 0; a < face_->eNoN; a++) { @@ -1063,6 +1153,13 @@ Array CappingSurface::update_element_position_global(int e, consts::Mech /// @return The Jacobian and normal vector. std::pair> CappingSurface::compute_jacobian_and_normal(const Array& xl, int e, int g) const { + if (xl.nrows() != cap_nsd_ || xl.ncols() != face_->eNoN) { + throw CappingSurfaceGeometryException("[CappingSurface::compute_jacobian_and_normal] xl has wrong dimensions: " + + std::to_string(xl.nrows()) + "x" + std::to_string(xl.ncols()) + + " (expected " + std::to_string(cap_nsd_) + "x" + + std::to_string(face_->eNoN) + ")."); + } + // Get the shape function derivatives for the Gauss point. Array Nx_g = face_->Nx.rslice(g); Array xXi(cap_nsd_, cap_insd_); @@ -1084,7 +1181,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)); } @@ -1092,6 +1189,11 @@ std::pair> CappingSurface::compute_jacobian_and_normal(co // Check if the initial normals are provided and if they are valid. if (normals_.ncols() > 0 && normals_.nrows() == cap_nsd_) { + if (e < 0 || e >= normals_.ncols()) { + throw CappingSurfaceGeometryException("[CappingSurface::compute_jacobian_and_normal] Element index e=" + + std::to_string(e) + " is out of bounds for normals_ (ncols=" + + std::to_string(normals_.ncols()) + ")."); + } Vector n0(cap_nsd_); for (int i = 0; i < cap_nsd_; i++) { @@ -1175,7 +1277,19 @@ void CappingSurface::compute_valM(consts::MechanicalConfigurationType cfg, const auto [Jac, n] = compute_jacobian_and_normal(xl, e, g); for (int a = 0; a < face_->eNoN; a++) { int gnNo_idx = face_->IEN(a, e); - int cap_a = gnNo_to_cap_local.at(gnNo_idx); + auto it = gnNo_to_cap_local.find(gnNo_idx); + if (it == gnNo_to_cap_local.end()) { + throw CappingSurfaceAssemblyException("[CappingSurface::compute_valM] IEN entry (element " + + std::to_string(e) + ", node " + std::to_string(a) + + ") contains invalid gnNo index " + std::to_string(gnNo_idx) + + " not found in cap face nodes."); + } + int cap_a = it->second; + if (cap_a < 0 || cap_a >= cap_nNo) { + throw CappingSurfaceAssemblyException("[CappingSurface::compute_valM] Invalid cap face-local index cap_a=" + + std::to_string(cap_a) + " (cap_nNo=" + std::to_string(cap_nNo) + + ")"); + } for (int i = 0; i < cap_nsd_; i++) { valM_(i, cap_a) += face_->N(a, g) * face_->w(g) * Jac * n(i); } @@ -1311,3 +1425,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..764877096 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). + + /// @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 (DIR coupling only). + int oned_ramp_steps_ = 0; ///< Number of ramp steps (0 = disabled) + double oned_ramp_ref_pressure_ = 0.0; ///< Reference pressure at step 0 + + /// @brief Under-relaxation factor for pressure passed to the 1D solver (DIR coupling only). + /// 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_sv1d_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 38ec0fe5d..7e910e711 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); @@ -2519,6 +2558,7 @@ GeneralSimulationParameters::GeneralSimulationParameters() set_parameter("Include_xml", "", !required, include_xml); set_parameter("Increment_in_saving_restart_files", 0, !required, increment_in_saving_restart_files); set_parameter("Increment_in_saving_VTK_files", 0, !required, increment_in_saving_vtk_files); + set_parameter("Initial_pressure_scalar", 0.0, !required, initial_pressure); set_parameter("Name_prefix_of_saved_VTK_files", "", !required, name_prefix_of_saved_vtk_files); set_parameter("Number_of_initialization_time_steps", 0, !required, number_of_initialization_time_steps, {0,int_inf}); diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index b726ba47f..dfcc8db3a 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -679,7 +679,47 @@ class svZeroDSolverInterfaceParameters : public ParameterLists bool value_set = false; }; -/// @brief Body force over a mesh using the "Add_BF" command. +//---------------------------------- +// 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; +}; + /// /// \code {.xml} /// @@ -763,6 +803,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: flow rate Q ramped from 0 to actual 3D Q (ref pressure not used). + // 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 +1586,8 @@ class EquationParameters : public ParameterLists svZeroDSolverInterfaceParameters svzerodsolver_interface_parameters; + svOneDSolverInterfaceParameters svonedsolver_interface_parameters; + DomainParameters* default_domain = nullptr; std::vector domains; @@ -1595,6 +1655,8 @@ class GeneralSimulationParameters : public ParameterLists Parameter spectral_radius_of_infinite_time_step; Parameter time_step_size; + Parameter initial_pressure; + Parameter include_xml; Parameter increment_in_saving_restart_files; diff --git a/Code/Source/solver/Simulation.cpp b/Code/Source/solver/Simulation.cpp index 7e6d56523..8cef36530 100644 --- a/Code/Source/solver/Simulation.cpp +++ b/Code/Source/solver/Simulation.cpp @@ -60,6 +60,12 @@ void Simulation::set_module_parameters() com_mod.nITs = general.number_of_initialization_time_steps.value(); com_mod.startTS = general.starting_time_step.value(); com_mod.dt = general.time_step_size.value(); + com_mod.have_initial_pressure = general.initial_pressure.defined(); // Set the flag to indicate whether initial pressure is provided in the input file + if (com_mod.have_initial_pressure) { + com_mod.initial_pressure = general.initial_pressure.value(); + } else { + com_mod.initial_pressure = 0.0; // Default pressure value if not provided + } com_mod.stopTrigName = general.searched_file_name_to_trigger_stop.value(); com_mod.ichckIEN = general.check_ien_order.value(); diff --git a/Code/Source/solver/baf_ini.cpp b/Code/Source/solver/baf_ini.cpp index 2b3b1b071..06d113bbb 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_subroutines.h" #include "fsils_api.hpp" #include "fils_struct.hpp" @@ -122,6 +123,11 @@ 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; + // Copy per-face 1D input file if present (svOneD coupling). + 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,19 +145,26 @@ 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"); } } } + // Create temporary SolutionStates from initialized 3D state. + SolutionStates init_solutions; + init_solutions.old.get_acceleration() = Ao; + init_solutions.old.get_displacement() = Do; + init_solutions.old.get_velocity() = Yo; + init_solutions.current.get_acceleration() = Ao; + init_solutions.current.get_displacement() = Do; + init_solutions.current.get_velocity() = Yo; + if (!com_mod.stFileFlag) { - // Create temporary SolutionStates for set_bc calls - SolutionStates temp_solutions; - temp_solutions.old.get_acceleration() = Ao; - temp_solutions.old.get_displacement() = Do; - temp_solutions.old.get_velocity() = Yo; - set_bc::rcr_init(com_mod, cm_mod, temp_solutions); + set_bc::rcr_init(com_mod, cm_mod, init_solutions); } if (com_mod.cplBC.useGenBC) { @@ -159,7 +172,11 @@ void baf_ini(Simulation* simulation, SolutionStates& solutions) } if (com_mod.cplBC.useSvZeroD) { - svZeroD::init_svZeroD(com_mod, cm_mod); + svZeroD::init_svZeroD(com_mod, cm_mod, init_solutions); + } + + if (com_mod.cplBC.useSv1D) { + svOneD::init_svOneD(com_mod, cm_mod, init_solutions); } // Initialize cap integration for Coupled boundary conditions @@ -763,10 +780,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)) { 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 234b27c74..32e160e20 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -331,6 +331,9 @@ void distribute(Simulation* simulation) cm.bcast(cm_mod, &com_mod.nEq); cm.bcast(cm_mod, &com_mod.dt); cm.bcast(cm_mod, &com_mod.precompDt); + cm.bcast(cm_mod, &com_mod.have_initial_pressure); + cm.bcast(cm_mod, &com_mod.initial_pressure); + cm.bcast(cm_mod, &com_mod.zeroAve); cm.bcast(cm_mod, &com_mod.cmmInit); @@ -550,24 +553,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.useSv1D); 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.useSv1D) { + // Broadcast the sv1D 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.sv1d_solver_interface.has_data); + cm.bcast(cm_mod, cplBC.sv1d_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/initialize.cpp b/Code/Source/solver/initialize.cpp index de31ec97e..17438db2a 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -937,15 +937,22 @@ void zero_init(Simulation* simulation, SolutionStates& solutions) } } - if (com_mod.Pinit.size() != 0) { + if (com_mod.Pinit.size() != 0) {// read initial pressure from vtk file. nodewise pressure value #ifdef debug_zero_init dmsg << "Initialize Yo to provided P solution"; #endif for (int a = 0; a < com_mod.tnNo; a++) { - for (int i = 0; i < nsd; i++) { - Yo(nsd,a) = com_mod.Pinit(a); - } + Yo(nsd,a) = com_mod.Pinit(a); + } + } else if (com_mod.have_initial_pressure) {// read initial pressure from solver.xml. scalar initial pressure value + #ifdef debug_zero_init + dmsg << "Initialize Yo to scalar initial pressure from solver.xml"; + #endif + for (int a = 0; a < com_mod.tnNo; a++) { + Yo(nsd,a) = com_mod.initial_pressure; } + std::string msg = "Initialize Yo(nsd,a) to initial_pressure value from solver.xml: " + std::to_string(com_mod.initial_pressure); + std::cout << msg << std::endl; } if (com_mod.Dinit.size() != 0) { diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 65b9aaaf1..e2fc91a7c 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -128,12 +128,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)); } @@ -236,38 +239,100 @@ 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). // } 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 sv1d_iface = com_mod.cplBC.sv1d_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 (ci_set) { - throw std::runtime_error(std::string("[read_bc] on face '") + face_name + - "' must define ."); - } + if (svzd_iface || sv1d_iface) { + // svZeroD and/or svOneD path: route each face individually based on its + // content. + + if (ci_has_block && ci_has_1d_file) { throw std::runtime_error( - std::string("[read_bc] With , each svZeroD-coupled face needs " - " with (Time_dependence Coupled) on face '") + - face_name + "'."); + std::string("[read_bc] on face '") + face_name + + "' defines both and . " + "Specify exactly one per face."); } + if (ci_has_block) { + // 0D face: route to svZeroD. + if (!svzd_iface) { + throw std::runtime_error( + std::string("[read_bc] Face '") + face_name + + "' specifies but no " + "is defined on the equation."); + } + 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 (!sv1d_iface) { + throw std::runtime_error( + std::string("[read_bc] Face '") + face_name + + "' specifies but no " + "is defined on the equation."); + } + 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) { + throw std::runtime_error( + std::string("[read_bc] A coupling method (e.g. svOneDSolver_interface coupling_type) " + "must be defined for Time_dependence Coupled on face '") + + face_name + "'."); + } + + } 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 either (for 0D coupling) or " + " (for 1D coupling)."); + } + + if (svzd_iface) { + throw std::runtime_error( + std::string("[read_bc] With , each svZeroD-coupled face needs " + " with " + "(Time_dependence Coupled) on face '") + + face_name + "'."); + } else { + throw std::runtime_error( + std::string("[read_bc] With , each 1D-coupled face needs " + " with " + "(Time_dependence Coupled) on face '") + + face_name + "'."); + } + } } else { - // Coupled BC to GenBC - if (ci_set) { + // genBC / cplBC path: no svZeroD or svOneD interface defined. + if (bc_params->coupling_interface.value_set) { throw std::runtime_error( - "[read_bc] is only valid when is defined on the equation."); + "[read_bc] is only valid when or " + " is defined on the equation."); } lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_cpl)); @@ -318,8 +383,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.useSv1D && !com_mod.cplBC.useSvZeroD) || + com_mod.cplBC.xo.size() != 0) { + throw std::runtime_error("[read_bc] RCR cannot be used in conjunction with cplBC (except alongside svOneD or svZeroD)."); } com_mod.cplBC.nFa = com_mod.cplBC.nFa + 1; lBc.cplBCptr = com_mod.cplBC.nFa - 1; @@ -394,67 +460,140 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, } - // Coupled BC (svZeroDSolver / CoupledBoundaryCondition) + // Coupled BC (svZeroDSolver or svOneDSolver via CoupledBoundaryCondition) 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) { - 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))) { + // 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) 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_sv1d_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_sv1d_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) { + // ------------------------------------------------------------------ + // svZeroD path (existing). + // ------------------------------------------------------------------ + + // Sanity check: must be defined + if (!bc_params->coupling_interface.value_set) { + throw std::runtime_error( + std::string("[read_bc] CoupledBoundaryCondition requires for face '") + face_name + + "'."); + } + + // Sanity check: must be defined + if (!bc_params->coupling_interface.svzerod_solver_block.defined()) { + throw std::runtime_error( + std::string("[read_bc] must define for face '") + face_name + + "'."); + } + + // 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."); } } @@ -466,10 +605,31 @@ 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 + // 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 (utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Coupled)) && + coupled_bc_type == BoundaryConditionType::bType_Dir) { + lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_flx)); + } + + // 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)); + // Also enable zero-perimeter for Coupled-DIR (svZeroD or svOneD DIR coupling). + if (!ltmp && + utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Coupled)) && + coupled_bc_type == BoundaryConditionType::bType_Dir) { + ltmp = true; + } if (bc_params->zero_out_perimeter.defined()) { ltmp = bc_params->zero_out_perimeter.value(); } @@ -1421,13 +1581,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 sv1d_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 != sv1d_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 '" + sv1d_type + "')."); + } + cplbc_type_str = sv1d_type; + cplBC.useSv1D = true; + cplBC.sv1d_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) { @@ -1444,7 +1629,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; } } @@ -1533,9 +1718,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.useSv1D && !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..35673ca27 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_subroutines.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 useSv1D can be true. + // Call each active solver independently. + if (cplBC.useSvZeroD) { + svZeroD::calc_svZeroD(com_mod, cm_mod, 'D'); + } + if (cplBC.useSv1D) { + svOneD::calc_svOneD(com_mod, cm_mod, 'D'); + } + if (!cplBC.useSvZeroD && !cplBC.useSv1D) { + 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.useSv1D) { + 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_sv1d_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 useSv1D can be true. + // Call each active solver independently. + if (cplBC.useSvZeroD) { + svZeroD::calc_svZeroD(com_mod, cm_mod, 'D'); + } + if (cplBC.useSv1D) { + svOneD::calc_svOneD(com_mod, cm_mod, 'D'); + } + if (!cplBC.useSvZeroD && !cplBC.useSv1D) { + 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/OneDSolverInterface.cpp b/Code/Source/solver/svOneD_interface/OneDSolverInterface.cpp new file mode 100644 index 000000000..da6e098c2 --- /dev/null +++ b/Code/Source/solver/svOneD_interface/OneDSolverInterface.cpp @@ -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 + +#include "OneDSolverInterface.h" + +#include +#include +#include +#include + +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 std::runtime_error(std::string("[OneDSolverInterface] Could not load shared library '") + + interface_lib + "': " + dlerror()); + } + + // 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 std::runtime_error(std::string("[OneDSolverInterface] Could not load symbol '") + + name + "': " + err); + } + 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 std::runtime_error("[OneDSolverInterface] initialize_1d not loaded"); + } + 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 std::runtime_error("[OneDSolverInterface] set_external_step_size_1d not loaded"); + } + set_external_step_size_1d_(problem_id, dt); +} + +void OneDSolverInterface::return_solution(int problem_id, double* solution, int size) +{ + if (!return_1d_solution_) { + throw std::runtime_error("[OneDSolverInterface] return_1d_solution not loaded"); + } + return_1d_solution_(problem_id, solution, size); +} + +void OneDSolverInterface::update_solution(int problem_id, double* solution, int size) +{ + if (!update_1d_solution_) { + throw std::runtime_error("[OneDSolverInterface] update_1d_solution not loaded"); + } + update_1d_solution_(problem_id, solution, size); +} + +void OneDSolverInterface::run_step(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 std::runtime_error("[OneDSolverInterface] run_1d_simulation_step_1d not loaded"); + } + // Copy coupling_type into a mutable buffer (shared-library uses char*). + std::vector ctype_buf(coupling_type.begin(), coupling_type.end()); + ctype_buf.push_back('\0'); + // Copy last_flag into a mutable single-character buffer. + char flag_buf[2] = { last_flag, '\0' }; + run_1d_simulation_step_1d_(problem_id, current_time, save_incr, + ctype_buf.data(), params, solution, + cpl_value, flag_buf, error_code); +} + +void OneDSolverInterface::extract_coupled_dof(int problem_id, int& coupled_dof, + const std::string& coupling_type) +{ + if (!extract_coupled_dof_) { + throw std::runtime_error("[OneDSolverInterface] extract_coupled_dof not loaded"); + } + // 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..29fcb3dd3 --- /dev/null +++ b/Code/Source/solver/svOneD_interface/OneDSolverInterface.h @@ -0,0 +1,90 @@ +// 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 uses dlopen/dlsym to load the libsvoned_interface shared library at runtime +/// and provides C++ method wrappers for each of its exported C functions. +/// +/// 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 input_file Path to the 1D solver .in file. + /// @param problem_id Output: problem identifier assigned by the solver. + /// @param system_size Output: total number of DOFs (nodes * 2: flow + area). + /// @param 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 problem_id Problem identifier. + /// @param current_time Current simulation time (start of the step). + /// @param save_incr VTK output interval (Increment_in_saving_VTK_files from solver.xml); + /// the 1D library decides internally whether to write output. + /// @param coupling_type "NEU" or "DIR". + /// @param params Array [N, t1, t2, ..., val1, val2, ...] where N=2. + /// @param solution In/out: solution vector updated after the step. + /// @param cpl_value Output: the BC value returned by the 1D solver + /// (pressure for NEU, flow for DIR). + /// @param last_flag 'L' for the final (committed) iteration, 'D' for + /// derivative / predictor steps. + /// @param error_code Output: non-zero on failure. + void run_step(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); + + // Public data members set after initialize(). + int problem_id_ = 0; + int system_size_ = 0; + + private: + void* library_handle_ = nullptr; + + // Function pointers to shared-library symbols. + void (*initialize_1d_)(const char*, int&, int&, const char*) = nullptr; + void (*set_external_step_size_1d_)(int, double) = nullptr; + void (*return_1d_solution_)(int, double*, int) = nullptr; + void (*update_1d_solution_)(int, double*, int) = nullptr; + void (*run_1d_simulation_step_1d_)(int, double, int, const char*, double*, + double*, double&, char*, int&) = nullptr; + void (*extract_coupled_dof_)(int, int&, char*) = nullptr; +}; + +#endif // ONEDSOLVER_INTERFACE_H diff --git a/Code/Source/solver/svOneD_subroutines.cpp b/Code/Source/solver/svOneD_subroutines.cpp new file mode 100644 index 000000000..eed2e7932 --- /dev/null +++ b/Code/Source/solver/svOneD_subroutines.cpp @@ -0,0 +1,451 @@ +// 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_step() 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_subroutines.h" + +#include +#include +#include +#include +#include +#include + +#include "ComMod.h" +#include "all_fun.h" +#include "consts.h" +#include "utils.h" +#include "svOneD_interface/OneDSolverInterface.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::ifstream(lib_base + ".so").good()) return lib_base + ".so"; + if (std::ifstream(lib_base + ".dylib").good()) 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, const SolutionStates& solutions) +{ + using namespace consts; + + auto& cplBC = com_mod.cplBC; + auto& solver_if = cplBC.sv1d_solver_interface; + auto& cm = com_mod.cm; + const int nProcs = cm.nProcs; + const int myRank = cm.taskId; + + if (!solver_if.has_data) { + throw std::runtime_error("[svOneD::init_svOneD] svOneD solver interface data is missing."); + } + + // 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(); + + const auto& face = com_mod.msh[bc.iM].fa[bc.iFa]; + const double area = face.area; + if (area > 0.0) { + const auto& Yo = solutions.old.get_velocity(); + const double P_init = all_fun::integ(com_mod, cm_mod, face, Yo, com_mod.nsd, + solutions, std::nullopt, false, + MechanicalConfigurationType::reference) / area; + st.P_prev_sent_old = P_init; + st.P_prev_sent_new = P_init; + st.P_neu_prev = P_init; + } + + oned_models.push_back(std::move(st)); + } + } + + if (oned_models.empty()) { + throw std::runtime_error("[svOneD::init_svOneD] No svOneD-coupled faces with input files found."); + } + + // ----- 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 std::runtime_error( + "[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."); + } + + // ----- 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); + + // Keep the BC pressure state consistent with the seeded pressure history. + eq.bc[st.iBc].coupled_bc.set_pressure(st.P_neu_prev); + } + + // ----- 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; + + st.interface->run_step(st.problem_id, t_old, save_incr, + st.coupling_type, params, + work_sol.data(), cpl_values[k], BCFlag, error_code); + + if (error_code != 0) { + throw std::runtime_error( + "[svOneD::calc_svOneD] 1D solver step for face '" + + bc.coupled_bc.get_oned_input_file() + "' failed with error code " + + std::to_string(error_code)); + } + + // 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_subroutines.h b/Code/Source/solver/svOneD_subroutines.h new file mode 100644 index 000000000..a18bc8596 --- /dev/null +++ b/Code/Source/solver/svOneD_subroutines.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, const SolutionStates& solutions); + +/// @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/svZeroD_interface.cpp b/Code/Source/solver/svZeroD_interface.cpp index 470f5fff8..d3c5c6c39 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_sv1d_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_sv1d_face()) { cpl.svZeroD_coupled_bc_idxs.emplace_back(iEq, iBc); } } @@ -220,7 +226,7 @@ void print_svZeroD(int* nSrfs, const std::vector& surfID, double Q[], doubl // init_svZeroD //-------------- // -void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod) +void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& solutions) { using namespace consts; @@ -252,6 +258,11 @@ void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod) int ids[2]; double in_out; + const bool have_init_flow = solver_interface.have_initial_flows; + const double init_flow = solver_interface.initial_flows; + const bool have_init_press = solver_interface.have_initial_pressures; + const double init_press = solver_interface.initial_pressures; + if (cm.mas(cm_mod)) { nsrflistCoupled.clear(); svzd_blk_names.clear(); @@ -265,6 +276,7 @@ void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod) throw std::runtime_error("ERROR: [init_svZeroD] svZeroDSolver interface data is missing."); } + #ifdef debug_init_svZeroD dmsg << "#### Use XML data #### "; #endif @@ -272,10 +284,6 @@ void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod) svzerod_library = solver_interface.solver_library; svzerod_file = solver_interface.configuration_file; - const bool have_init_flow = solver_interface.have_initial_flows; - const double init_flow = solver_interface.initial_flows; - const bool have_init_press = solver_interface.have_initial_pressures; - const double init_press = solver_interface.initial_pressures; for (int s = 0; s < numCoupledSrfs; ++s) { bcType* bc = nullptr; @@ -324,6 +332,7 @@ void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod) if (!nth_coupled_bc(com_mod, s, &bc)) { throw std::runtime_error("ERROR: [init_svZeroD] Internal error resolving Coupled BC for initialization."); } + if (bc->coupled_bc.get_bc_type() == consts::BoundaryConditionType::bType_Neu) { if (have_init_flow) { lpn_state_y[sol_IDs[2 * s]] = init_flow; @@ -353,6 +362,26 @@ void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod) } } + for (int s = 0; s < numCoupledSrfs; ++s) { + bcType* bc = nullptr; + if (!nth_coupled_bc(com_mod, s, &bc)) { + throw std::runtime_error("ERROR: [init_svZeroD] Internal error resolving Coupled BC for initialization history."); + } + + // For cap coupling in parallel, flow integration requires all ranks. + if (cm.seq() || !bc->coupled_bc.has_cap()) { + bc->coupled_bc.compute_flowrates(com_mod, cm_mod, solutions); + } + bc->coupled_bc.compute_pressures(com_mod, cm_mod, solutions); + const double Q_init = bc->coupled_bc.get_Qo(); + const double P_init = bc->coupled_bc.get_Po(); + bc->coupled_bc.set_Q_prev_sent(Q_init); + bc->coupled_bc.set_Q_input_prev(Q_init, Q_init); + bc->coupled_bc.set_P_prev_sent(P_init, P_init); + bc->coupled_bc.set_P_neu_prev(P_init); + bc->coupled_bc.set_pressure(P_init); + } + // Broadcast initial values to follower processes if (!cm.seq()) { // Coupled BCs - broadcast Neumann pressures (one scalar bcast per BC). @@ -386,6 +415,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 +453,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 = bc->coupled_bc.get_oned_ramp_ref_pressure(); // 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 +526,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 +550,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 +587,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/svZeroD_interface.h b/Code/Source/solver/svZeroD_interface.h index 1de624f85..d807891ee 100644 --- a/Code/Source/solver/svZeroD_interface.h +++ b/Code/Source/solver/svZeroD_interface.h @@ -17,7 +17,7 @@ void get_coupled_QP(ComMod& com_mod, double QCoupled[], double QnCoupled[], doub void print_svZeroD(int* nSrfs, const std::vector& surfID, double Q[], double P[]); -void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod); +void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& solutions); void calc_svZeroD(ComMod& com_mod, const CmMod& cm_mod, char BCFlag); diff --git a/Code/Source/solver/txt.cpp b/Code/Source/solver/txt.cpp index 61e3060a0..bbe6b033d 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_subroutines.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; + } else { + // In mixed-coupling simulations both useSvZeroD and useSv1D can be true. + // Call each active solver independently. + if (cplBC.useSvZeroD) { + svZeroD::calc_svZeroD(com_mod, cm_mod, 'L'); + } + + if (cplBC.useSv1D) { + // 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_sv1d_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.useSv1D) { + 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/pipe_svOneD_2faces/solver.xml b/tests/cases/fluid/pipe_svOneD_2faces/solver.xml new file mode 100644 index 000000000..bc15d6bad --- /dev/null +++ b/tests/cases/fluid/pipe_svOneD_2faces/solver.xml @@ -0,0 +1,142 @@ + + + + + + + false + 3 + 100 + 0.001 + 0.50 + STOP_SIM + + 1 + result + 10 + 1 + + 100 + 0 + + 1 + 0 + 0 + + + + + + mesh-complete/mesh-complete.mesh.vtu + + + mesh-complete/mesh-surfaces/lumen_inlet.vtp + + + + mesh-complete/mesh-surfaces/lumen_outlet1.vtp + + + + mesh-complete/mesh-surfaces/lumen_outlet2.vtp + + + + mesh-complete/mesh-surfaces/lumen_wall.vtp + + + + + + 1 + 3 + 10 + 1e-3 + 0.2 + + 1.06 + + 0.04 + + + + true + true + true + true + + + + + fsils + + 10 + 3 + 500 + 1e-3 + 1e-3 + 1e-3 + 50 + + + + + semi-implicit + /path/to/svOneDSolver/build/lib/libsvOneDSolver_interface.so + + + + + Dir + Unsteady + lumen_inlet.flw + true + true + + + + + Neu + Coupled + + 1dmodel1.in + + + + + + Neu + Coupled + + 1dmodel2.in + + + + + + Dir + Steady + 0.0 + + + + +