diff --git a/_codeql_detected_source_root b/_codeql_detected_source_root new file mode 120000 index 000000000..945c9b46d --- /dev/null +++ b/_codeql_detected_source_root @@ -0,0 +1 @@ +. \ No newline at end of file diff --git a/docs/pages/activation_functions.md b/docs/pages/activation_functions.md new file mode 100644 index 000000000..5ddb0c91f --- /dev/null +++ b/docs/pages/activation_functions.md @@ -0,0 +1,162 @@ +# Activation Functions for Chamber Models + +## Overview + +Chamber models in svZeroDSolver support different activation functions to model the time-varying elastance of cardiac chambers. The activation function determines how the chamber contracts and relaxes over the cardiac cycle. + +## Available Activation Functions + +### 1. Half Cosine Activation + +A half cosine wave activation function to model contraction: + +``` +A(t) = -0.5 * cos(2π * t_contract / t_twitch) + 0.5 if t_contract ≤ t_twitch +A(t) = 0 otherwise +``` + +where `t_contract = max(0, t - t_active)` + +**Parameters:** +- `t_active`: Time when activation begins within the cardiac cycle +- `t_twitch`: Duration of the contraction twitch + +**Example usage in JSON:** +```json +{ + "type": "ChamberElastanceInductor", + "name": "ventricle", + "values": { + "Emax": 1.057, + "Emin": 0.091, + "Vrd": 26.1, + "Vrs": 18.0, + "Impedance": 0.000351787, + "activation_function_type": "half_cosine", + "activation_function_values": { + "t_active": 0.2, + "t_twitch": 0.3 + } + } +} +``` + +### 2. Piecewise Cosine Activation + +A piecewise activation function that models separate contraction and relaxation phases: + +``` +φ(t) = 0.5 * [1 - cos(π * (t - t_C) / T_C)] during contraction +φ(t) = 0.5 * [1 + cos(π * (t - t_R) / T_R)] during relaxation +φ(t) = 0 otherwise +``` + +**Parameters:** +- `contract_start`: Time when contraction starts +- `relax_start`: Time when relaxation starts +- `contract_duration`: Duration of contraction phase +- `relax_duration`: Duration of relaxation phase + +**Example usage in JSON:** +```json +{ + "type": "LinearElastanceChamber", + "name": "left_atrium", + "values": { + "Emax": 199.95, + "Epass": 260.59, + "Vrest": 26.24, + "activation_function_type": "piecewise_cosine", + "activation_function_values": { + "contract_start": 0.025, + "relax_start": 0.08625, + "contract_duration": 0.06125, + "relax_duration": 0.18375 + } + } +} +``` + +### 3. Two Hill Activation + +A more flexible and physiologically realistic activation function based on the two-hill model. This allows for more precise control over the activation waveform shape. + +The activation is computed as: +``` +A(t) = C * [g₁(t) / (1 + g₁(t))] * [1 / (1 + g₂(t))] +``` + +where: +``` +g₁(t) = (t_shifted / τ₁)^m₁ +g₂(t) = (t_shifted / τ₂)^m₂ +t_shifted = (t - t_shift) mod T_cardiac +``` + +C is a normalization constant ensuring maximum activation is 1. + +**Parameters:** +- `t_shift`: Time shift parameter to align activation with cardiac cycle +- `tau_1`: Time constant for first hill (controls rise) +- `tau_2`: Time constant for second hill (controls fall) +- `m1`: Exponent for first hill (controls steepness of rise) +- `m2`: Exponent for second hill (controls steepness of fall) + +**Example usage in JSON:** +```json +{ + "type": "ChamberElastanceInductor", + "name": "ventricle", + "values": { + "Emax": 1.057, + "Emin": 0.091, + "Vrd": 26.1, + "Vrs": 18.0, + "Impedance": 0.000351787, + "activation_function_type": "two_hill", + "activation_function_values": { + "t_shift": 0.15, + "tau_1": 0.25, + "tau_2": 0.35, + "m1": 1.5, + "m2": 10.0 + } + } +} +``` + +## Switching Between Activation Functions + +**Both `ChamberElastanceInductor` and `LinearElastanceChamber` require the `activation_function_type` parameter to be explicitly specified as a string:** + +- `activation_function_type: "half_cosine"` - Half Cosine activation +- `activation_function_type: "piecewise_cosine"` - Piecewise Cosine activation +- `activation_function_type: "two_hill"` - Two Hill activation + +All activation function-specific parameters must be grouped under `activation_function_values` as a nested dictionary. + +### Important Notes + +- The `activation_function_type` parameter is **required** for both chamber types +- Only string values are accepted (e.g., "half_cosine", not numeric values) +- Parameters must be specified within the `activation_function_values` nested dictionary + +## Elastance Calculation + +For all activation functions, the time-varying elastance E(t) is computed as: + +**For ChamberElastanceInductor:** +``` +E(t) = (Emax - Emin) * A(t) + Emin +Vrest(t) = (1 - A(t)) * (Vrd - Vrs) + Vrs +``` + +**For LinearElastanceChamber:** +``` +E(t) = Epass + Emax * φ(t) +``` + +## References + +The two-hill activation function is described in: +- Kaiser, A. D., et al. (2022). "A design-based model of the aortic valve for fluid-structure interaction." Biomechanics and Modeling in Mechanobiology. https://link.springer.com/article/10.1007/s10439-022-03047-3 diff --git a/src/model/ActivationFunction.h b/src/model/ActivationFunction.h new file mode 100644 index 000000000..ff1a3950b --- /dev/null +++ b/src/model/ActivationFunction.h @@ -0,0 +1,383 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +/** + * @file ActivationFunction.h + * @brief Activation function classes for cardiac chamber models + */ + +#ifndef SVZERODSOLVER_MODEL_ACTIVATIONFUNCTION_HPP_ +#define SVZERODSOLVER_MODEL_ACTIVATIONFUNCTION_HPP_ + +#include +#include +#include +#include +#include + +/** + * @brief Enum for activation function types + */ +enum class ActivationType { + HALF_COSINE = 0, + PIECEWISE_COSINE = 1, + TWO_HILL = 2 +}; + +/** + * @brief Parameters for creating activation functions + */ +struct ActivationFunctionParams { + ActivationType type; + double cardiac_period; + + // Half cosine parameters + double t_active = 0.0; + double t_twitch = 0.0; + + // Piecewise cosine parameters + double contract_start = 0.0; + double relax_start = 0.0; + double contract_duration = 0.0; + double relax_duration = 0.0; + + // Two hill parameters + double t_shift = 0.0; + double tau_1 = 0.0; + double tau_2 = 0.0; + double m1 = 0.0; + double m2 = 0.0; +}; + +/** + * @brief Base class for activation functions + * + * Activation functions compute the activation value (between 0 and 1) at a + * given time point within a cardiac cycle. These are used to modulate + * chamber elastance over time. + */ +class ActivationFunction { + public: + /** + * @brief Virtual destructor + */ + virtual ~ActivationFunction() = default; + + /** + * @brief Compute activation value at given time + * + * @param time Current time + * @return Activation value between 0 and 1 + */ + virtual double compute(double time) = 0; + + /** + * @brief Get the type of activation function + * + * @return Type of activation function + */ + virtual ActivationType get_type() const = 0; + + /** + * @brief Factory method to create activation function from parameters + * + * @param params Parameters for creating the activation function + * @return Unique pointer to the created activation function + */ + static std::unique_ptr create( + const ActivationFunctionParams& params); + + protected: + double cardiac_period_; ///< Cardiac cycle period +}; + +/** + * @brief Half cosine activation function + * + * This implements the activation function used in the original + * ChamberElastanceInductor. The activation follows a half cosine wave + * during the contraction period. + * + * \f[ + * A(t) = \begin{cases} + * -\frac{1}{2}\cos(2\pi t_{contract}/t_{twitch}) + \frac{1}{2}, & \text{if } t_{contract} \le t_{twitch} \\ + * 0, & \text{otherwise} + * \end{cases} + * \f] + * + * where \f$t_{contract} = \max(0, t_{in\_cycle} - t_{active})\f$ + */ +class HalfCosineActivation : public ActivationFunction { + public: + /** + * @brief Construct a new Half Cosine Activation object + * + * @param t_active Time when activation begins within cardiac cycle + * @param t_twitch Duration of the contraction twitch + * @param cardiac_period Cardiac cycle period + */ + HalfCosineActivation(double t_active, double t_twitch, double cardiac_period) + : t_active_(t_active), t_twitch_(t_twitch) { + cardiac_period_ = cardiac_period; + } + + /** + * @brief Compute activation value + * + * @param time Current time + * @return Activation value between 0 and 1 + */ + double compute(double time) override { + double t_in_cycle = std::fmod(time, cardiac_period_); + + double t_contract = 0.0; + if (t_in_cycle >= t_active_) { + t_contract = t_in_cycle - t_active_; + } + + double act = 0.0; + if (t_contract <= t_twitch_) { + act = -0.5 * std::cos(2.0 * M_PI * t_contract / t_twitch_) + 0.5; + } + + return act; + } + + /** + * @brief Get the type of activation function + * + * @return ActivationType::HALF_COSINE + */ + ActivationType get_type() const override { + return ActivationType::HALF_COSINE; + } + + private: + double t_active_; ///< Time when activation begins + double t_twitch_; ///< Duration of contraction twitch +}; + +/** + * @brief Piecewise cosine activation function + * + * This implements the activation function from the LinearElastanceChamber + * (Regazzoni chamber model). The activation consists of separate contraction + * and relaxation phases, each following a cosine curve. + * + * \f[ + * \phi(t, t_C, t_R, T_C, T_R) = \begin{cases} + * \frac{1}{2}\left[1 - \cos\left(\frac{\pi}{T_C} \bmod(t - t_C, T_{HB})\right)\right], + * & \text{if } 0 \le \bmod(t - t_C, T_{HB}) < T_C \\ + * \frac{1}{2}\left[1 + \cos\left(\frac{\pi}{T_R} \bmod(t - t_R, T_{HB})\right)\right], + * & \text{if } 0 \le \bmod(t - t_R, T_{HB}) < T_R \\ + * 0, & \text{otherwise} + * \end{cases} + * \f] + */ +class PiecewiseCosineActivation : public ActivationFunction { + public: + /** + * @brief Construct a new Piecewise Cosine Activation object + * + * @param contract_start Time when contraction starts + * @param relax_start Time when relaxation starts + * @param contract_duration Duration of contraction phase + * @param relax_duration Duration of relaxation phase + * @param cardiac_period Cardiac cycle period + */ + PiecewiseCosineActivation(double contract_start, double relax_start, + double contract_duration, double relax_duration, + double cardiac_period) + : contract_start_(contract_start), + relax_start_(relax_start), + contract_duration_(contract_duration), + relax_duration_(relax_duration) { + cardiac_period_ = cardiac_period; + } + + /** + * @brief Compute activation value + * + * @param time Current time + * @return Activation value between 0 and 1 + */ + double compute(double time) override { + double phi = 0.0; + + double piecewise_condition = std::fmod(time - contract_start_, cardiac_period_); + + if (0.0 <= piecewise_condition && piecewise_condition < contract_duration_) { + phi = 0.5 * (1.0 - std::cos((M_PI * piecewise_condition) / contract_duration_)); + } else { + piecewise_condition = std::fmod(time - relax_start_, cardiac_period_); + if (0.0 <= piecewise_condition && piecewise_condition < relax_duration_) { + phi = 0.5 * (1.0 + std::cos((M_PI * piecewise_condition) / relax_duration_)); + } + } + + return phi; + } + + /** + * @brief Get the type of activation function + * + * @return ActivationType::PIECEWISE_COSINE + */ + ActivationType get_type() const override { + return ActivationType::PIECEWISE_COSINE; + } + + private: + double contract_start_; ///< Start time of contraction + double relax_start_; ///< Start time of relaxation + double contract_duration_; ///< Duration of contraction + double relax_duration_; ///< Duration of relaxation +}; + +/** + * @brief Two hill activation function + * + * This implements the two-hill activation function which provides more + * flexible and physiologically realistic waveforms. See + * https://link.springer.com/article/10.1007/s10439-022-03047-3 + * + * The activation is computed as: + * \f[ + * A(t) = C \cdot \frac{g_1(t)}{1 + g_1(t)} \cdot \frac{1}{1 + g_2(t)} + * \f] + * + * where: + * \f[ + * g_1(t) = \left(\frac{t_{shifted}}{\tau_1}\right)^{m_1}, \quad + * g_2(t) = \left(\frac{t_{shifted}}{\tau_2}\right)^{m_2} + * \f] + * + * and \f$t_{shifted} = (t - t_{shift}) \bmod T_{cardiac}\f$, and \f$C\f$ is a + * normalization constant to ensure max activation is 1. + */ +class TwoHillActivation : public ActivationFunction { + public: + /** + * @brief Construct a new Two Hill Activation object + * + * @param t_shift Time shift parameter + * @param tau_1 Time constant for first hill + * @param tau_2 Time constant for second hill + * @param m1 Exponent for first hill + * @param m2 Exponent for second hill + * @param cardiac_period Cardiac cycle period (needed for normalization) + */ + TwoHillActivation(double t_shift, double tau_1, double tau_2, double m1, + double m2, double cardiac_period) + : t_shift_(t_shift), + tau_1_(tau_1), + tau_2_(tau_2), + m1_(m1), + m2_(m2), + normalization_factor_(1.0), + normalization_initialized_(false) { + cardiac_period_ = cardiac_period; + initialize_normalization(cardiac_period); + } + + /** + * @brief Compute activation value + * + * @param time Current time + * @return Activation value between 0 and 1 + */ + double compute(double time) override { + if (!normalization_initialized_) { + throw std::runtime_error("TwoHillActivation: Normalization not initialized"); + } + + double t_in_cycle = std::fmod(time, cardiac_period_); + + // Compute shifted time (handle negative modulo) + double t_shifted = std::fmod(t_in_cycle - t_shift_, cardiac_period_); + t_shifted = (t_shifted >= 0.0) ? t_shifted : t_shifted + cardiac_period_; + + // Compute hill functions + double g1 = std::pow(t_shifted / tau_1_, m1_); + double g2 = std::pow(t_shifted / tau_2_, m2_); + + // Compute activation with normalization + double act = normalization_factor_ * (g1 / (1.0 + g1)) * (1.0 / (1.0 + g2)); + + return act; + } + + /** + * @brief Get the type of activation function + * + * @return ActivationType::TWO_HILL + */ + ActivationType get_type() const override { + return ActivationType::TWO_HILL; + } + + private: + /** + * @brief Initialize normalization factor + * + * Computes the maximum value of the two-hill function over one cardiac + * cycle to normalize the activation to [0, 1]. + * + * @param cardiac_period Period of cardiac cycle + */ + void initialize_normalization(double cardiac_period) { + // Time step for numerical integration when finding maximum activation + // Value chosen to balance accuracy and computational cost + constexpr double NORMALIZATION_DT = 1e-5; + + double max_value = 0.0; + + for (double t_temp = 0.0; t_temp < cardiac_period; t_temp += NORMALIZATION_DT) { + double g1 = std::pow(t_temp / tau_1_, m1_); + double g2 = std::pow(t_temp / tau_2_, m2_); + double two_hill_val = (g1 / (1.0 + g1)) * (1.0 / (1.0 + g2)); + + max_value = std::max(max_value, two_hill_val); + } + + normalization_factor_ = 1.0 / max_value; + normalization_initialized_ = true; + } + + double t_shift_; ///< Time shift parameter + double tau_1_; ///< Time constant for first hill + double tau_2_; ///< Time constant for second hill + double m1_; ///< Exponent for first hill + double m2_; ///< Exponent for second hill + double normalization_factor_; ///< Normalization constant + bool normalization_initialized_; ///< Flag for normalization +}; + +/** + * @brief Factory method implementation for creating activation functions + */ +inline std::unique_ptr ActivationFunction::create( + const ActivationFunctionParams& params) { + switch (params.type) { + case ActivationType::HALF_COSINE: + return std::make_unique( + params.t_active, params.t_twitch, params.cardiac_period); + + case ActivationType::PIECEWISE_COSINE: + return std::make_unique( + params.contract_start, params.relax_start, + params.contract_duration, params.relax_duration, params.cardiac_period); + + case ActivationType::TWO_HILL: + return std::make_unique( + params.t_shift, params.tau_1, params.tau_2, + params.m1, params.m2, params.cardiac_period); + + default: + throw std::runtime_error( + "ActivationFunction::create: Invalid activation type " + + std::to_string(static_cast(params.type))); + } +} + +#endif // SVZERODSOLVER_MODEL_ACTIVATIONFUNCTION_HPP_ diff --git a/src/model/BlockType.h b/src/model/BlockType.h index a57972647..2284b33b5 100644 --- a/src/model/BlockType.h +++ b/src/model/BlockType.h @@ -30,8 +30,8 @@ enum class BlockType { chamber_elastance_inductor = 14, chamber_sphere = 15, blood_vessel_CRL = 16, - piecewise_cosine_chamber = 17, - piecewise_valve = 18 + piecewise_valve = 17, + linear_elastance_chamber = 18 }; /** diff --git a/src/model/CMakeLists.txt b/src/model/CMakeLists.txt index a76fc23e3..c6128e356 100644 --- a/src/model/CMakeLists.txt +++ b/src/model/CMakeLists.txt @@ -29,11 +29,12 @@ set(CXXSRCS ResistiveJunction.cpp ValveTanh.cpp WindkesselBC.cpp - PiecewiseCosineChamber.cpp + LinearElastanceChamber.cpp PiecewiseValve.cpp ) set(HDRS + ActivationFunction.h Block.h BlockType.h BloodVessel.h @@ -58,7 +59,7 @@ set(HDRS ResistiveJunction.h ValveTanh.h WindkesselBC.h - PiecewiseCosineChamber.h + LinearElastanceChamber.h PiecewiseValve.h ) diff --git a/src/model/ChamberElastanceInductor.cpp b/src/model/ChamberElastanceInductor.cpp index 440835e8a..a1ff39bc1 100644 --- a/src/model/ChamberElastanceInductor.cpp +++ b/src/model/ChamberElastanceInductor.cpp @@ -9,6 +9,11 @@ void ChamberElastanceInductor::setup_dofs(DOFHandler& dofhandler) { void ChamberElastanceInductor::update_constant( SparseSystem& system, std::vector& parameters) { + // Initialize activation function on first call + if (!activation_func_) { + initialize_activation_function(parameters); + } + double L = parameters[global_param_ids[ParamId::IMPEDANCE]]; // Eq 0: P_in - E(t)(Vc - Vrest) = 0 @@ -40,22 +45,24 @@ void ChamberElastanceInductor::get_elastance_values( double Emin = parameters[global_param_ids[ParamId::EMIN]]; double Vrd = parameters[global_param_ids[ParamId::VRD]]; double Vrs = parameters[global_param_ids[ParamId::VRS]]; - double t_active = parameters[global_param_ids[ParamId::TACTIVE]]; - double t_twitch = parameters[global_param_ids[ParamId::TTWITCH]]; - - auto T_cardiac = model->cardiac_cycle_period; - auto t_in_cycle = fmod(model->time, T_cardiac); - - double t_contract = 0; - if (t_in_cycle >= t_active) { - t_contract = t_in_cycle - t_active; - } - - double act = 0; - if (t_contract <= t_twitch) { - act = -0.5 * cos(2 * M_PI * t_contract / t_twitch) + 0.5; - } + + // Compute activation using the activation function + double act = activation_func_->compute(model->time); Vrest = (1.0 - act) * (Vrd - Vrs) + Vrs; Elas = (Emax - Emin) * act + Emin; } + +void ChamberElastanceInductor::initialize_activation_function( + std::vector& parameters) { + // Get required activation_type parameter + int activation_type_int = static_cast( + parameters[global_param_ids[ParamId::ACTIVATION_TYPE]]); + activation_params_.type = static_cast(activation_type_int); + + // Set cardiac period + activation_params_.cardiac_period = model->cardiac_cycle_period; + + // Use factory method to create activation function + activation_func_ = ActivationFunction::create(activation_params_); +} diff --git a/src/model/ChamberElastanceInductor.h b/src/model/ChamberElastanceInductor.h index 527947cc7..bf69e6ed6 100644 --- a/src/model/ChamberElastanceInductor.h +++ b/src/model/ChamberElastanceInductor.h @@ -8,7 +8,9 @@ #define SVZERODSOLVER_MODEL_CHAMBERELASTANCEINDUCTOR_HPP_ #include +#include +#include "ActivationFunction.h" #include "Block.h" #include "Model.h" #include "SparseSystem.h" @@ -152,9 +154,8 @@ class ChamberElastanceInductor : public Block { {"Emin", InputParameter()}, {"Vrd", InputParameter()}, {"Vrs", InputParameter()}, - {"t_active", InputParameter()}, - {"t_twitch", InputParameter()}, - {"Impedance", InputParameter()}}) {} + {"Impedance", InputParameter()}, + {"activation_type", InputParameter()}}) {} /** * @brief Local IDs of the parameters @@ -165,9 +166,8 @@ class ChamberElastanceInductor : public Block { EMIN = 1, VRD = 2, VRS = 3, - TACTIVE = 4, - TTWITCH = 5, - IMPEDANCE = 6 + IMPEDANCE = 4, + ACTIVATION_TYPE = 5 }; /** @@ -211,6 +211,26 @@ class ChamberElastanceInductor : public Block { private: double Elas; // Chamber Elastance double Vrest; // Rest Volume + std::unique_ptr activation_func_; // Activation function + ActivationFunctionParams activation_params_; // Activation function parameters + + public: + /** + * @brief Set activation function parameters from JSON configuration + * + * @param params Activation function parameters + */ + void set_activation_params(const ActivationFunctionParams& params) { + activation_params_ = params; + } + + private: + /** + * @brief Initialize the activation function based on parameters + * + * @param parameters Parameters of the model + */ + void initialize_activation_function(std::vector& parameters); /** * @brief Update the elastance functions which depend on time diff --git a/src/model/LinearElastanceChamber.cpp b/src/model/LinearElastanceChamber.cpp new file mode 100644 index 000000000..b1f73cf46 --- /dev/null +++ b/src/model/LinearElastanceChamber.cpp @@ -0,0 +1,69 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#include "LinearElastanceChamber.h" + +void LinearElastanceChamber::setup_dofs(DOFHandler& dofhandler) { + // Internal variable is chamber volume + Block::setup_dofs_(dofhandler, 3, {"Vc"}); +} + +void LinearElastanceChamber::update_constant(SparseSystem& system, + std::vector& parameters) { + // Initialize activation function on first call + if (!activation_func_) { + initialize_activation_function(parameters); + } + + // Eq 0: P_in - E(t)(Vc - Vrest) = 0 + system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = 1.0; + + // Eq 1: P_in - P_out = 0 + system.F.coeffRef(global_eqn_ids[1], global_var_ids[0]) = 1.0; + system.F.coeffRef(global_eqn_ids[1], global_var_ids[2]) = -1.0; + + // Eq 2: Q_in - Q_out - dVc = 0 + system.F.coeffRef(global_eqn_ids[2], global_var_ids[1]) = 1.0; + system.F.coeffRef(global_eqn_ids[2], global_var_ids[3]) = -1.0; + system.E.coeffRef(global_eqn_ids[2], global_var_ids[4]) = -1.0; +} + +void LinearElastanceChamber::update_time(SparseSystem& system, + std::vector& parameters) { + get_elastance_values(parameters); + + // Eq 0: P_in - E(t)(Vc - Vrest) = P_in - E(t)*Vc + E(t)*Vrest = 0 + system.F.coeffRef(global_eqn_ids[0], global_var_ids[4]) = -Elas; + system.C.coeffRef(global_eqn_ids[0]) = + Elas * parameters[global_param_ids[ParamId::VREST]]; +} + +void LinearElastanceChamber::get_elastance_values( + std::vector& parameters) { + double Emax = parameters[global_param_ids[ParamId::EMAX]]; + double Epass = parameters[global_param_ids[ParamId::EPASS]]; + + // Compute activation using the activation function + double phi = activation_func_->compute(model->time); + + Elas = Epass + Emax * phi; +} + +void LinearElastanceChamber::initialize_activation_function( + std::vector& parameters) { + // Check if activation_type parameter is provided (required parameter) + if (global_param_ids.size() <= ParamId::ACTIVATION_TYPE) { + throw std::runtime_error( + "LinearElastanceChamber: activation_type parameter is required"); + } + + int activation_type_int = static_cast( + parameters[global_param_ids[ParamId::ACTIVATION_TYPE]]); + activation_params_.type = static_cast(activation_type_int); + + // Set cardiac period + activation_params_.cardiac_period = model->cardiac_cycle_period; + + // Use factory method to create activation function + activation_func_ = ActivationFunction::create(activation_params_); +} diff --git a/src/model/PiecewiseCosineChamber.h b/src/model/LinearElastanceChamber.h similarity index 76% rename from src/model/PiecewiseCosineChamber.h rename to src/model/LinearElastanceChamber.h index b9a461503..0c1000546 100644 --- a/src/model/PiecewiseCosineChamber.h +++ b/src/model/LinearElastanceChamber.h @@ -2,25 +2,27 @@ // University of California, and others. SPDX-License-Identifier: BSD-3-Clause /** - * @file PiecewiseCosineChamber.h - * @brief model::PiecewiseCosineChamber source file + * @file LinearElastanceChamber.h + * @brief model::LinearElastanceChamber source file */ -#ifndef SVZERODSOLVER_MODEL_PIECEWISECOSINECHAMBER_HPP_ -#define SVZERODSOLVER_MODEL_PIECEWISECOSINECHAMBER_HPP_ +#ifndef SVZERODSOLVER_MODEL_LINEARELASTANCECHAMBER_HPP_ +#define SVZERODSOLVER_MODEL_LINEARELASTANCECHAMBER_HPP_ #include +#include +#include "ActivationFunction.h" #include "Block.h" #include "Model.h" #include "SparseSystem.h" #include "debug.h" /** - * @brief Cardiac chamber with piecewise elastance and inductor. + * @brief Cardiac chamber with linear elastance (no inductor). * * Models a cardiac chamber as a time-varying capacitor (elastance with - * specified resting volumes) and an inductor. See \cite Regazzoni2022 + * specified resting volumes) without inductance. See \cite Regazzoni2022 * * This chamber block can be connected to other blocks using junctions. * @@ -115,24 +117,21 @@ * * `Vc`: Chamber volume * */ -class PiecewiseCosineChamber : public Block { +class LinearElastanceChamber : public Block { public: /** - * @brief Construct a new BloodVessel object + * @brief Construct a new LinearElastanceChamber object * * @param id Global ID of the block * @param model The model to which the block belongs */ - PiecewiseCosineChamber(int id, Model* model) - : Block(id, model, BlockType::piecewise_cosine_chamber, + LinearElastanceChamber(int id, Model* model) + : Block(id, model, BlockType::linear_elastance_chamber, BlockClass::chamber, {{"Emax", InputParameter()}, {"Epass", InputParameter()}, {"Vrest", InputParameter()}, - {"contract_start", InputParameter()}, - {"relax_start", InputParameter()}, - {"contract_duration", InputParameter()}, - {"relax_duration", InputParameter()}}) {} + {"activation_type", InputParameter(false)}}) {} /** * @brief Local IDs of the parameters @@ -142,10 +141,7 @@ class PiecewiseCosineChamber : public Block { EMAX = 0, EPASS = 1, VREST = 2, - CSTART = 3, - RSTART = 4, - CDUR = 5, - RDUR = 6 + ACTIVATION_TYPE = 3 }; /** @@ -188,6 +184,26 @@ class PiecewiseCosineChamber : public Block { private: double Elas; // Chamber Elastance + std::unique_ptr activation_func_; // Activation function + ActivationFunctionParams activation_params_; // Activation function parameters + + public: + /** + * @brief Set activation function parameters from JSON configuration + * + * @param params Activation function parameters + */ + void set_activation_params(const ActivationFunctionParams& params) { + activation_params_ = params; + } + + private: + /** + * @brief Initialize the activation function based on parameters + * + * @param parameters Parameters of the model + */ + void initialize_activation_function(std::vector& parameters); /** * @brief Update the elastance functions which depend on time @@ -197,4 +213,4 @@ class PiecewiseCosineChamber : public Block { void get_elastance_values(std::vector& parameters); }; -#endif // SVZERODSOLVER_MODEL_PIECEWISECOSINECHAMBER_HPP_ +#endif // SVZERODSOLVER_MODEL_LINEARELASTANCECHAMBER_HPP_ diff --git a/src/model/Model.cpp b/src/model/Model.cpp index a527161a5..2344a6b60 100644 --- a/src/model/Model.cpp +++ b/src/model/Model.cpp @@ -31,7 +31,7 @@ Model::Model() { {"ChamberElastanceInductor", block_factory()}, {"BloodVesselCRL", block_factory()}, {"PiecewiseValve", block_factory()}, - {"PiecewiseCosineChamber", block_factory()}}; + {"LinearElastanceChamber", block_factory()}}; } Model::~Model() {} diff --git a/src/model/Model.h b/src/model/Model.h index 6a444e93a..c837af5e6 100644 --- a/src/model/Model.h +++ b/src/model/Model.h @@ -29,10 +29,10 @@ #include "DOFHandler.h" #include "FlowReferenceBC.h" #include "Junction.h" +#include "LinearElastanceChamber.h" #include "Node.h" #include "OpenLoopCoronaryBC.h" #include "Parameter.h" -#include "PiecewiseCosineChamber.h" #include "PiecewiseValve.h" #include "PressureReferenceBC.h" #include "ResistanceBC.h" diff --git a/src/model/PiecewiseCosineChamber.cpp b/src/model/PiecewiseCosineChamber.cpp deleted file mode 100644 index 010016576..000000000 --- a/src/model/PiecewiseCosineChamber.cpp +++ /dev/null @@ -1,62 +0,0 @@ -// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the -// University of California, and others. SPDX-License-Identifier: BSD-3-Clause - -#include "PiecewiseCosineChamber.h" - -void PiecewiseCosineChamber::setup_dofs(DOFHandler& dofhandler) { - // Internal variable is chamber volume - Block::setup_dofs_(dofhandler, 3, {"Vc"}); -} - -void PiecewiseCosineChamber::update_constant(SparseSystem& system, - std::vector& parameters) { - // Eq 0: P_in - E(t)(Vc - Vrest) = 0 - system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = 1.0; - - // Eq 1: P_in - P_out = 0 - system.F.coeffRef(global_eqn_ids[1], global_var_ids[0]) = 1.0; - system.F.coeffRef(global_eqn_ids[1], global_var_ids[2]) = -1.0; - - // Eq 2: Q_in - Q_out - dVc = 0 - system.F.coeffRef(global_eqn_ids[2], global_var_ids[1]) = 1.0; - system.F.coeffRef(global_eqn_ids[2], global_var_ids[3]) = -1.0; - system.E.coeffRef(global_eqn_ids[2], global_var_ids[4]) = -1.0; -} - -void PiecewiseCosineChamber::update_time(SparseSystem& system, - std::vector& parameters) { - get_elastance_values(parameters); - - // Eq 0: P_in - E(t)(Vc - Vrest) = P_in - E(t)*Vc + E(t)*Vrest = 0 - system.F.coeffRef(global_eqn_ids[0], global_var_ids[4]) = -1 * Elas; - system.C.coeffRef(global_eqn_ids[0]) = - Elas * parameters[global_param_ids[ParamId::VREST]]; -} - -void PiecewiseCosineChamber::get_elastance_values( - std::vector& parameters) { - double Emax = parameters[global_param_ids[ParamId::EMAX]]; - double Epass = parameters[global_param_ids[ParamId::EPASS]]; - double Vrest = parameters[global_param_ids[ParamId::VREST]]; - double contract_start = parameters[global_param_ids[ParamId::CSTART]]; - double relax_start = parameters[global_param_ids[ParamId::RSTART]]; - double contract_duration = parameters[global_param_ids[ParamId::CDUR]]; - double relax_duration = parameters[global_param_ids[ParamId::RDUR]]; - - auto T_HB = model->cardiac_cycle_period; - - double phi = 0; - - auto piecewise_condition = fmod(model->time - contract_start, T_HB); - - if (0 <= piecewise_condition && piecewise_condition < contract_duration) { - phi = 0.5 * (1 - cos((M_PI * piecewise_condition) / contract_duration)); - } else { - piecewise_condition = fmod(model->time - relax_start, T_HB); - if (0 <= piecewise_condition && piecewise_condition < relax_duration) { - phi = 0.5 * (1 + cos((M_PI * piecewise_condition) / relax_duration)); - } - } - - Elas = Epass + Emax * phi; -} diff --git a/src/solve/SimulationParameters.cpp b/src/solve/SimulationParameters.cpp index b914c401d..dac52629d 100644 --- a/src/solve/SimulationParameters.cpp +++ b/src/solve/SimulationParameters.cpp @@ -79,6 +79,9 @@ int generate_block(Model& model, const nlohmann::json& block_params_json, } } } else { + // First pass: collect all missing parameters + std::vector missing_params; + for (const auto& block_param : block->input_params) { // Time parameter is read at the same time as time-dependent value if (block_param.first.compare("t") == 0) { @@ -90,42 +93,85 @@ int generate_block(Model& model, const nlohmann::json& block_params_json, continue; } - // Get vector parameter + // Check vector parameter if (block_param.second.is_array) { - // Get parameter vector + // Check parameter vector std::vector val; err = get_param_vector(block_params_json, block_param.first, block_param.second, val); if (err) { - throw std::runtime_error("Array parameter " + block_param.first + - " is mandatory in " + block_type + - " block " + static_cast(name)); + missing_params.push_back(block_param.first + " (array)"); } - // Get time vector + // Check time vector InputParameter t_param{false, true}; std::vector time; err = get_param_vector(block_params_json, "t", t_param, time); if (err) { - throw std::runtime_error("Array parameter t is mandatory in " + - block_type + " block " + - static_cast(name)); + // Only add "t" if not already in the list + if (std::find(missing_params.begin(), missing_params.end(), "t (array)") == missing_params.end()) { + missing_params.push_back("t (array)"); + } } - - // Add parameters to model - new_id = model.add_parameter(time, val, periodic); } - // Get scalar parameter + // Check scalar parameter else { double val; err = get_param_scalar(block_params_json, block_param.first, block_param.second, val); if (err) { - throw std::runtime_error("Scalar parameter " + block_param.first + - " is mandatory in " + block_type + - " block " + static_cast(name)); + missing_params.push_back(block_param.first); } + } + } + + // If any parameters are missing, throw error with complete list + if (!missing_params.empty()) { + std::string error_msg = "Missing required parameter(s) in " + block_type + + " block " + static_cast(name) + ": "; + for (size_t i = 0; i < missing_params.size(); ++i) { + error_msg += missing_params[i]; + if (i < missing_params.size() - 1) { + error_msg += ", "; + } + } + throw std::runtime_error(error_msg); + } + + // Second pass: actually read and store the parameters + for (const auto& block_param : block->input_params) { + // Time parameter is read at the same time as time-dependent value + if (block_param.first.compare("t") == 0) { + continue; + } + + // Skip reading parameters that are not a number + if (!block_param.second.is_number) { + continue; + } + + // Get vector parameter + if (block_param.second.is_array) { + // Get parameter vector + std::vector val; + get_param_vector(block_params_json, block_param.first, + block_param.second, val); + + // Get time vector + InputParameter t_param{false, true}; + std::vector time; + get_param_vector(block_params_json, "t", t_param, time); + + // Add parameters to model + new_id = model.add_parameter(time, val, periodic); + } + + // Get scalar parameter + else { + double val; + get_param_scalar(block_params_json, block_param.first, + block_param.second, val); // Add parameter to model new_id = model.add_parameter(val); @@ -544,7 +590,141 @@ void create_chambers( const auto& chamber_config = JsonWrapper(config, component, "name", i); std::string chamber_type = chamber_config["type"]; std::string chamber_name = chamber_config["name"]; - generate_block(model, chamber_config["values"], chamber_type, chamber_name); + + // Create a mutable copy of the values to handle activation function parameters + nlohmann::json chamber_values = chamber_config["values"]; + + // Prepare activation function parameters + ActivationFunctionParams act_params; + act_params.cardiac_period = 0.0; // Will be set later from model + + // Check if activation function type is specified - it's required + if (!chamber_values.contains("activation_function_type")) { + throw std::runtime_error("Missing required parameter 'activation_function_type' in chamber " + + chamber_name + ". Must be one of: 'half_cosine', 'piecewise_cosine', 'two_hill'"); + } + + if (!chamber_values["activation_function_type"].is_string()) { + throw std::runtime_error("Parameter 'activation_function_type' must be a string in chamber " + + chamber_name + ". Must be one of: 'half_cosine', 'piecewise_cosine', 'two_hill'"); + } + + std::string act_type_str = chamber_values["activation_function_type"]; + + // Map string to activation type + if (act_type_str == "half_cosine") { + act_params.type = ActivationType::HALF_COSINE; + } else if (act_type_str == "piecewise_cosine") { + act_params.type = ActivationType::PIECEWISE_COSINE; + } else if (act_type_str == "two_hill") { + act_params.type = ActivationType::TWO_HILL; + } else { + throw std::runtime_error("Unknown activation_function_type: '" + act_type_str + + "' in chamber " + chamber_name + + ". Must be one of: 'half_cosine', 'piecewise_cosine', 'two_hill'"); + } + + // Check if activation function values are specified in nested dict + if (!chamber_values.contains("activation_function_values")) { + throw std::runtime_error("Missing required parameter 'activation_function_values' in chamber " + + chamber_name + ". This nested dictionary must contain parameters for the " + + act_type_str + " activation function."); + } + + if (!chamber_values["activation_function_values"].is_object()) { + throw std::runtime_error("Parameter 'activation_function_values' must be an object/dictionary in chamber " + + chamber_name); + } + + const auto& act_values = chamber_values["activation_function_values"]; + std::vector missing_params; + + // Extract and validate activation function parameters based on type + if (act_params.type == ActivationType::HALF_COSINE) { + // Required parameters for half_cosine + if (!act_values.contains("t_active")) missing_params.push_back("t_active"); + if (!act_values.contains("t_twitch")) missing_params.push_back("t_twitch"); + + if (!missing_params.empty()) { + std::string error_msg = "Missing required activation_function_values for 'half_cosine' in chamber " + + chamber_name + ": "; + for (size_t i = 0; i < missing_params.size(); i++) { + error_msg += missing_params[i]; + if (i < missing_params.size() - 1) error_msg += ", "; + } + throw std::runtime_error(error_msg); + } + + act_params.t_active = act_values["t_active"]; + act_params.t_twitch = act_values["t_twitch"]; + + } else if (act_params.type == ActivationType::PIECEWISE_COSINE) { + // Required parameters for piecewise_cosine + if (!act_values.contains("contract_start")) missing_params.push_back("contract_start"); + if (!act_values.contains("relax_start")) missing_params.push_back("relax_start"); + if (!act_values.contains("contract_duration")) missing_params.push_back("contract_duration"); + if (!act_values.contains("relax_duration")) missing_params.push_back("relax_duration"); + + if (!missing_params.empty()) { + std::string error_msg = "Missing required activation_function_values for 'piecewise_cosine' in chamber " + + chamber_name + ": "; + for (size_t i = 0; i < missing_params.size(); i++) { + error_msg += missing_params[i]; + if (i < missing_params.size() - 1) error_msg += ", "; + } + throw std::runtime_error(error_msg); + } + + act_params.contract_start = act_values["contract_start"]; + act_params.relax_start = act_values["relax_start"]; + act_params.contract_duration = act_values["contract_duration"]; + act_params.relax_duration = act_values["relax_duration"]; + + } else if (act_params.type == ActivationType::TWO_HILL) { + // Required parameters for two_hill + if (!act_values.contains("t_shift")) missing_params.push_back("t_shift"); + if (!act_values.contains("tau_1")) missing_params.push_back("tau_1"); + if (!act_values.contains("tau_2")) missing_params.push_back("tau_2"); + if (!act_values.contains("m1")) missing_params.push_back("m1"); + if (!act_values.contains("m2")) missing_params.push_back("m2"); + + if (!missing_params.empty()) { + std::string error_msg = "Missing required activation_function_values for 'two_hill' in chamber " + + chamber_name + ": "; + for (size_t i = 0; i < missing_params.size(); i++) { + error_msg += missing_params[i]; + if (i < missing_params.size() - 1) error_msg += ", "; + } + throw std::runtime_error(error_msg); + } + + act_params.t_shift = act_values["t_shift"]; + act_params.tau_1 = act_values["tau_1"]; + act_params.tau_2 = act_values["tau_2"]; + act_params.m1 = act_values["m1"]; + act_params.m2 = act_values["m2"]; + } + + // Remove the nested dict (don't flatten to avoid parameter conflicts) + chamber_values.erase("activation_function_values"); + + // Create the chamber block + int block_id = generate_block(model, chamber_values, chamber_type, chamber_name); + + // Set activation parameters on the chamber + auto chamber_block = model.get_block(chamber_name); + if (chamber_type == "ChamberElastanceInductor") { + auto chamber = dynamic_cast(chamber_block); + if (chamber) { + chamber->set_activation_params(act_params); + } + } else if (chamber_type == "LinearElastanceChamber") { + auto chamber = dynamic_cast(chamber_block); + if (chamber) { + chamber->set_activation_params(act_params); + } + } + DEBUG_MSG("Created chamber " << chamber_name); } } diff --git a/src/solve/SimulationParameters.h b/src/solve/SimulationParameters.h index 512a5aaac..34e09d1d1 100644 --- a/src/solve/SimulationParameters.h +++ b/src/solve/SimulationParameters.h @@ -12,6 +12,9 @@ #include #include +#include "ActivationFunction.h" +#include "ChamberElastanceInductor.h" +#include "LinearElastanceChamber.h" #include "Model.h" #include "State.h" #include "debug.h" diff --git a/src/solve/csv_writer.cpp b/src/solve/csv_writer.cpp index 4e521597a..8dbdda49b 100644 --- a/src/solve/csv_writer.cpp +++ b/src/solve/csv_writer.cpp @@ -45,7 +45,7 @@ std::string to_vessel_csv(const std::vector& times, if (dynamic_cast(block) == nullptr && dynamic_cast(block) == nullptr && - dynamic_cast(block) == nullptr && + dynamic_cast(block) == nullptr && dynamic_cast(block) == nullptr && dynamic_cast(block) == nullptr) { continue; diff --git a/tests/cases/chamber_elastance_inductor.json b/tests/cases/chamber_elastance_inductor.json index 9b756054e..bac17b8ad 100644 --- a/tests/cases/chamber_elastance_inductor.json +++ b/tests/cases/chamber_elastance_inductor.json @@ -101,9 +101,12 @@ "Emin": 0.091, "Vrd": 26.1, "Vrs": 18.0, - "t_active": 0.2, - "t_twitch": 0.3, - "Impedance": 0.000351787 + "Impedance": 0.000351787, + "activation_function_type": "half_cosine", + "activation_function_values": { + "t_active": 0.2, + "t_twitch": 0.3 + } } } ], diff --git a/tests/cases/chamber_two_hill_activation.json b/tests/cases/chamber_two_hill_activation.json new file mode 100644 index 000000000..d4c2fbb54 --- /dev/null +++ b/tests/cases/chamber_two_hill_activation.json @@ -0,0 +1,119 @@ +{ + "description": { + "description of test case": "pulsatile pressure -> valve -> chamber (two hill activation) -> valve -> RC -> R", + "analytical results": [ + "Boundary conditions:", + "inlet:", + "pressure = sinusoid from 0 to 10 over 0.17s, then 0 for the rest of cycle", + "outlet:", + "Resistance: Rd = 1000, Pd = 64" + ] + }, + "simulation_parameters": { + "number_of_cardiac_cycles": 30, + "number_of_time_pts_per_cardiac_cycle": 101, + "absolute_tolerance": 1e-9, + "output_all_cycles": false, + "output_variable_based": true + }, + "boundary_conditions": [ + { + "bc_name": "INLET", + "bc_type": "PRESSURE", + "bc_values": { + "P": [ + 0.0, 0.003414699, 0.4075651, 1.431949, 2.938218, 4.722943, 6.545085, + 8.158555, 9.345445, 9.945458, 9.87756, 9.15092, 7.863676, 6.189676, + 4.355004, 2.607442, 1.183009, 0.274081, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0 + ], + "t": [ + 0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.11, + 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.23, + 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3, 0.31, 0.32, 0.33, 0.34, 0.35, + 0.36, 0.37, 0.38, 0.39, 0.4, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, + 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, + 0.6, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.7, 0.71, + 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.8, 0.81, 0.82, 0.83, + 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, + 0.96, 0.97, 0.98, 0.99, 1.0 + ] + } + }, + { + "bc_name": "OUTLET", + "bc_type": "RESISTANCE", + "bc_values": { + "Pd": 64.0, + "R": 1000.0 + } + } + ], + "junctions": [], + "vessels": [ + { + "boundary_conditions": { + "outlet": "OUTLET" + }, + "vessel_id": 0, + "vessel_length": 10.0, + "vessel_name": "vessel", + "zero_d_element_type": "BloodVessel", + "zero_d_element_values": { + "C": 0.0001, + "R_poiseuille": 100.0 + } + } + ], + "valves": [ + { + "type": "ValveTanh", + "name": "valve0", + "params": { + "Rmax": 100000.0, + "Rmin": 100.0, + "Steepness": 100.0, + "upstream_block": "INLET", + "downstream_block": "ventricle" + } + }, + { + "type": "ValveTanh", + "name": "valve1", + "params": { + "Rmax": 100000.0, + "Rmin": 100.0, + "Steepness": 100.0, + "upstream_block": "ventricle", + "downstream_block": "vessel" + } + } + ], + "chambers": [ + { + "type": "ChamberElastanceInductor", + "name": "ventricle", + "values": { + "Emax": 1.057, + "Emin": 0.091, + "Vrd": 26.1, + "Vrs": 18.0, + "Impedance": 0.000351787, + "activation_function_type": "two_hill", + "activation_function_values": { + "t_shift": 0.15, + "tau_1": 0.25, + "tau_2": 0.35, + "m1": 1.5, + "m2": 10.0 + } + } + } + ], + "initial_condition": { + "Vc:ventricle": 96.07 + } +} diff --git a/tests/cases/linear_chamber_two_hill_activation.json b/tests/cases/linear_chamber_two_hill_activation.json new file mode 100644 index 000000000..5b0f6bf8f --- /dev/null +++ b/tests/cases/linear_chamber_two_hill_activation.json @@ -0,0 +1,142 @@ +{ + "description": { + "description of test case": "LinearElastanceChamber with two hill activation function", + "analytical results": [ + "Testing LinearElastanceChamber with two hill activation" + ] + }, + "simulation_parameters": { + "number_of_cardiac_cycles": 30, + "number_of_time_pts_per_cardiac_cycle": 689, + "output_variable_based": true, + "output_all_cycles": false, + "cardiac_period": 0.68900516753 + }, + "vessels": [ + { + "vessel_id": 1, + "vessel_length": 10.0, + "vessel_name": "pul_artery", + "zero_d_element_type": "BloodVessel", + "zero_d_element_values": { + "C": 0.0, + "L": 0.0, + "R_poiseuille": 0.0 + } + }, + { + "vessel_id": 5, + "vessel_length": 10.0, + "vessel_name": "sys_artery", + "zero_d_element_type": "BloodVessel", + "zero_d_element_values": { + "C": 0.001138164597527442, + "R_poiseuille": 596.82, + "L": 6.665 + } + }, + { + "vessel_id": 7, + "vessel_length": 10.0, + "vessel_name": "sys_vein", + "zero_d_element_type": "BloodVessel", + "zero_d_element_values": { + "C": 0.045011252813952737, + "R_poiseuille": 249.42, + "L": 0.6665 + } + } + ], + "junctions": [ + { + "inlet_blocks": ["sys_vein"], + "junction_name": "J1", + "junction_type": "NORMAL_JUNCTION", + "outlet_blocks": ["right_atrium"] + }, + { + "inlet_vessels": [5], + "junction_name": "J2", + "junction_type": "NORMAL_JUNCTION", + "outlet_vessels": [7] + } + ], + "boundary_conditions": [], + "chambers": [ + { + "type": "LinearElastanceChamber", + "name": "right_atrium", + "values": { + "Emax": 199.95, + "Epass": 89.80375933024839, + "Vrest": 41.680015938842274, + "activation_function_type": "two_hill", + "activation_function_values": { + "t_shift": 0.025, + "tau_1": 0.06, + "tau_2": 0.15, + "m1": 2.0, + "m2": 8.0 + } + } + }, + { + "type": "LinearElastanceChamber", + "name": "right_ventricle", + "values": { + "Emax": 1662.0158240056528, + "Epass": 40.85565535747109, + "Vrest": 72.05452710344869, + "activation_function_type": "two_hill", + "activation_function_values": { + "t_shift": 0.207, + "tau_1": 0.08, + "tau_2": 0.22, + "m1": 2.5, + "m2": 12.0 + } + } + } + ], + "valves": [ + { + "type": "PiecewiseValve", + "name": "tricuspid", + "params": { + "Rmax": 6665, + "Rmin": 6.665, + "upstream_block": "right_atrium", + "downstream_block": "right_ventricle" + } + }, + { + "type": "PiecewiseValve", + "name": "pulmonary", + "params": { + "Rmax": 6665, + "Rmin": 6.665, + "upstream_block": "right_ventricle", + "downstream_block": "pul_artery" + } + }, + { + "type": "PiecewiseValve", + "name": "aortic", + "params": { + "Rmax": 6665, + "Rmin": 6.665, + "upstream_block": "right_ventricle", + "downstream_block": "sys_artery" + } + } + ], + "initial_condition": { + "Vc:right_ventricle": 128.58981029386334, + "Vc:right_atrium": 76.8340776729488, + "pressure:aortic:sys_artery": 84604.18388331511, + "pressure:J2:sys_vein": 31311.64989129829, + "pressure:pulmonary:pul_artery": 20525.08438550143, + "flow:sys_artery:J2": 91.00177508885831, + "flow:sys_vein:J1": 112.86832799795421 + } +} diff --git a/tests/cases/piecewise_Chamber_and_Valve.json b/tests/cases/piecewise_Chamber_and_Valve.json index e4ec38b93..1efe3d2a0 100644 --- a/tests/cases/piecewise_Chamber_and_Valve.json +++ b/tests/cases/piecewise_Chamber_and_Valve.json @@ -126,55 +126,67 @@ "boundary_conditions": [], "chambers": [ { - "type": "PiecewiseCosineChamber", + "type": "LinearElastanceChamber", "name": "right_atrium", "values": { "Emax": 199.95, "Epass": 89.80375933024839, "Vrest": 41.680015938842274, - "contract_start": 0.025, - "relax_start": 0.08625, - "contract_duration": 0.06125, - "relax_duration": 0.18375 + "activation_function_type": "piecewise_cosine", + "activation_function_values": { + "contract_start": 0.025, + "relax_start": 0.08625, + "contract_duration": 0.06125, + "relax_duration": 0.18375 + } } }, { - "type": "PiecewiseCosineChamber", + "type": "LinearElastanceChamber", "name": "right_ventricle", "values": { "Emax": 1662.0158240056528, "Epass": 40.85565535747109, "Vrest": 72.05452710344869, - "contract_start": 0.207, - "relax_start": 0.29625, - "contract_duration": 0.08925, - "relax_duration": 0.26975 + "activation_function_type": "piecewise_cosine", + "activation_function_values": { + "contract_start": 0.207, + "relax_start": 0.29625, + "contract_duration": 0.08925, + "relax_duration": 0.26975 + } } }, { - "type": "PiecewiseCosineChamber", + "type": "LinearElastanceChamber", "name": "left_atrium", "values": { "Emax": 199.95, "Epass": 260.59064494602634, "Vrest": 26.23534455334443, - "contract_start": 0.025, - "relax_start": 0.08625, - "contract_duration": 0.06125, - "relax_duration": 0.18375 + "activation_function_type": "piecewise_cosine", + "activation_function_values": { + "contract_start": 0.025, + "relax_start": 0.08625, + "contract_duration": 0.06125, + "relax_duration": 0.18375 + } } }, { - "type": "PiecewiseCosineChamber", + "type": "LinearElastanceChamber", "name": "left_ventricle", "values": { "Emax": 17837.499965252825, "Epass": 122.82901158252674, "Vrest": 32.41857776349184, - "contract_start": 0.207, - "relax_start": 0.29625, - "contract_duration": 0.08925, - "relax_duration": 0.26975 + "activation_function_type": "piecewise_cosine", + "activation_function_values": { + "contract_start": 0.207, + "relax_start": 0.29625, + "contract_duration": 0.08925, + "relax_duration": 0.26975 + } } } ],