Split ChamberElastanceInductor into linear and exponential variants

Extract the exponential passive P-V relation (for atrial chambers) into
a separate ChamberElastanceInductorExponential class that inherits from
ChamberElastanceInductor. This removes if/else branching on Kxp==0 and
makes the two physical models (linear vs exponential passive P-V)
explicit through the type system.

- ChamberElastanceInductor: linear E(t)*(V-Vrest) model
- ChamberElastanceInductorExponential: adds Kxp, Kxv, Vaso parameters
  for exponential passive pressure (Sankaran 2012, Menon 2023)

Base class provides a protected constructor for derived classes and
virtual get_elastance_values() for the elastance computation.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: mrp089

src/model/BlockType.h

@@ -35,7 +35,8 @@ enum class BlockType {
   open_loop_coronary_var_res_bc = 19,
   open_loop_coronary_detailed_bc = 20,
   blood_vessel_rc = 21,
-  closed_loop_heart_pulmonary_smooth = 22
+  closed_loop_heart_pulmonary_smooth = 22,
+  chamber_elastance_inductor_exponential = 23
 };
 
 /**

src/model/CMakeLists.txt

@@ -14,6 +14,7 @@ set(CXXSRCS
   BloodVesselJunction.cpp
   ChamberSphere.cpp
   ChamberElastanceInductor.cpp
+  ChamberElastanceInductorExponential.cpp
   ClosedLoopCoronaryBC.cpp
   ClosedLoopCoronaryLeftBC.cpp
   ClosedLoopCoronaryRightBC.cpp
@@ -46,6 +47,7 @@ set(HDRS
   BloodVesselJunction.h
   ChamberSphere.h
   ChamberElastanceInductor.h
+  ChamberElastanceInductorExponential.h
   ClosedLoopCoronaryBC.h
   ClosedLoopCoronaryLeftBC.h
   ClosedLoopCoronaryRightBC.h

src/model/ChamberElastanceInductor.cpp

@@ -2,8 +2,9 @@
 // University of California, and others. SPDX-License-Identifier: BSD-3-Clause
 #include "ChamberElastanceInductor.h"
 
+#include "Model.h"
+
 void ChamberElastanceInductor::setup_dofs(DOFHandler& dofhandler) {
-  // Internal variable is chamber volume
   Block::setup_dofs_(dofhandler, 3, {"Vc"});
 }
 
@@ -29,37 +30,9 @@ void ChamberElastanceInductor::update_time(SparseSystem& system,
                                            std::vector<double>& 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;
-
-  // In exponential passive mode (Kxp > 0), C[0] is set in update_solution.
-  // In standard linear mode, set C[0] here.
-  double Kxp = parameters[global_param_ids[ParamId::KXP]];
-  if (Kxp <= 0.0) {
-    system.C.coeffRef(global_eqn_ids[0]) = Elas * Vrest;
-  }
-}
-
-void ChamberElastanceInductor::update_solution(
-    SparseSystem& system, std::vector<double>& parameters,
-    const Eigen::Matrix<double, Eigen::Dynamic, 1>& y,
-    const Eigen::Matrix<double, Eigen::Dynamic, 1>& dy) {
-  // --- Exponential passive P-V (atrial mode, Kxp > 0) ---
-  double Kxp = parameters[global_param_ids[ParamId::KXP]];
-  if (Kxp <= 0.0) return;
-
-  double Kxv = parameters[global_param_ids[ParamId::KXV]];
-  double Vaso = parameters[global_param_ids[ParamId::VASO]];
-  double Emax = parameters[global_param_ids[ParamId::EMAX]];
-  double Vc = y[global_var_ids[4]];
-
-  double exp_term = std::exp(Kxv * (Vc - Vaso));
-  double psi = Kxp * (exp_term - 1.0);
-  double psi_d = Kxp * Kxv * exp_term;
-
-  system.C(global_eqn_ids[0]) = act_ * Emax * Vaso + psi * (act_ - 1.0);
-  system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[4]) =
-      psi_d * (act_ - 1.0);
+  // Eq 0: F[0][4] = -E(t), C[0] = E(t)*Vrest
+  system.F.coeffRef(global_eqn_ids[0], global_var_ids[4]) = -Elas;
+  system.C.coeffRef(global_eqn_ids[0]) = Elas * Vrest;
 }
 
 void ChamberElastanceInductor::get_elastance_values(
@@ -70,15 +43,8 @@ void ChamberElastanceInductor::get_elastance_values(
   double Vrs = parameters[global_param_ids[ParamId::VRS]];
 
   act_ = activation_func_->compute(model->time);
-
-  double Kxp = parameters[global_param_ids[ParamId::KXP]];
-  if (Kxp > 0.0) {
-    Elas = act_ * Emax;
-    Vrest = parameters[global_param_ids[ParamId::VASO]];
-  } else {
-    Vrest = (1.0 - act_) * (Vrd - Vrs) + Vrs;
-    Elas = (Emax - Emin) * act_ + Emin;
-  }
+  Vrest = (1.0 - act_) * (Vrd - Vrs) + Vrs;
+  Elas = (Emax - Emin) * act_ + Emin;
 }
 
 void ChamberElastanceInductor::set_activation_function(

src/model/ChamberElastanceInductor.h

@@ -13,12 +13,11 @@
 
 #include "ActivationFunction.h"
 #include "Block.h"
-#include "Model.h"
 #include "SparseSystem.h"
 #include "debug.h"
 
 /**
- * @brief Cardiac chamber with elastance and inductor.
+ * @brief Cardiac chamber with linear elastance and inductor.
  *
  * Models a cardiac chamber as a time-varying capacitor (elastance with
  * specified resting volumes) and an inductor. See \cite kerckhoffs2007coupling
@@ -55,91 +54,26 @@
  * Q_{in}-Q_{out}-\dot{V}_c=0
  * \f]
  *
- * ### Local contributions
- *
- * \f[
- * \mathbf{y}^{e}=\left[\begin{array}{lllll}P_{in} & Q_{in} &
- * P_{out} & Q_{out} & V_c\end{array}\right]^{T} \f]
- *
- * \f[
- * \mathbf{E}^{e}=\left[\begin{array}{ccccc}
- * 0 & 0 & 0 & 0 & 0\\
- * 0 & 0 & 0 & -L & 0\\
- * 0 & 0 & 0 & 0 & -1
- * \end{array}\right]
- * \f]
- *
- * \f[
- * \mathbf{F}^{e}=\left[\begin{array}{ccccc}
- * 1 & 0 &  0 & 0  & E(t) \\
- * 1 & 0 & -1 & 0  & 0 \\
- * 0 & 1 &  0 & -1 & 0
- * \end{array}\right]
- * \f]
- *
- * \f[
- * \mathbf{c}^{e}=\left[\begin{array}{c}
- * E(t)V_{rest} \\
- * 0 \\
- * 0
- * \end{array}\right]
- * \f]
- *
- * In the above equations,
+ * where
  *
  * \f[
  * V_{rest}(t)= \{1-A(t)\}(V_{rd}-V_{rs})+V_{rs}
  * \f]
  *
  * \f[
- * A(t)=-\frac{1}{2}cos(2 \pi T_{contract}/T_{twitch})
- * \f]
- *
- * \f[
  * E(t)=(E_{max}-E_{min})A(t) + E_{min}
  * \f]
  *
- *
  * ### Parameters
  *
- * Parameter sequence for constructing this block
- *
  * * `0` Emax: Maximum elastance
  * * `1` Emin: Minimum elastance
  * * `2` Vrd: Rest diastolic volume
  * * `3` Vrs: Rest systolic volume
- * * `4` t_active: Activation time
- * * `5` t_twitch: Twitch time
- * * `6` Impedance: Impedance of the outflow
- *
- * ### Usage in json configuration file
- *
- *     "chambers": [
- *         {
- *             "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",
- *                 "t_active": 0.2,
- *                 "t_twitch": 0.3
- *             }
- *         }
- *     ],
- *     "initial_condition": {
- *         "Vc:ventricle": 96.07
- *     }
+ * * `4` Impedance: Impedance of the outflow
  *
  * ### Internal variables
  *
- * Names of internal variables in this block's output:
- *
  * * `Vc`: Chamber volume
  *
  */
@@ -158,93 +92,46 @@ class ChamberElastanceInductor : public Block {
                {"Emin", InputParameter()},
                {"Vrd", InputParameter()},
                {"Vrs", InputParameter()},
-               {"Impedance", InputParameter()},
-               {"Kxp", InputParameter(true)},
-               {"Kxv", InputParameter(true)},
-               {"Vaso", InputParameter(true)}}) {}
+               {"Impedance", InputParameter()}}) {}
 
   /**
    * @brief Local IDs of the parameters
-   *
    */
   enum ParamId {
     EMAX = 0,
     EMIN = 1,
     VRD = 2,
     VRS = 3,
     IMPEDANCE = 4,
-    KXP = 5,  ///< Passive pressure scaling (optional, 0 = linear mode)
-    KXV = 6,  ///< Passive volume scaling (optional)
-    VASO = 7  ///< Passive resting volume (optional)
   };
 
-  /**
-   * @brief Set up the degrees of freedom (DOF) of the block
-   *
-   * Set global_var_ids and global_eqn_ids of the element based on the
-   * number of equations and the number of internal variables of the
-   * element.
-   *
-   * @param dofhandler Degree-of-freedom handler to register variables and
-   * equations at
-   */
   void setup_dofs(DOFHandler& dofhandler);
-
-  /**
-   * @brief Update the constant contributions of the element in a sparse
-   system
-   *
-   * @param system System to update contributions at
-   * @param parameters Parameters of the model
-   */
   void update_constant(SparseSystem& system, std::vector<double>& parameters);
-
-  /**
-   * @brief Update the time-dependent contributions of the element in a sparse
-   * system
-   *
-   * @param system System to update contributions at
-   * @param parameters Parameters of the model
-   */
   void update_time(SparseSystem& system, std::vector<double>& parameters);
 
-  /**
-   * @brief Update the solution-dependent contributions of the element in a
-   * sparse system. Only active when Kxp > 0 (exponential passive atrial mode).
-   */
-  void update_solution(SparseSystem& system, std::vector<double>& parameters,
-                       const Eigen::Matrix<double, Eigen::Dynamic, 1>& y,
-                       const Eigen::Matrix<double, Eigen::Dynamic, 1>& dy);
-
-  /**
-   * @brief Number of triplets of element
-   *
-   * Number of triplets that the element contributes to the global system
-   * (relevant for sparse memory reservation)
-   */
-  TripletsContributions num_triplets{6, 2, 1};
+  TripletsContributions num_triplets{6, 2, 0};
 
- private:
-  double Elas;        // Chamber Elastance
-  double Vrest;       // Rest Volume
-  double act_ = 0.0;  // Last computed activation
-  std::unique_ptr<ActivationFunction> activation_func_;  // Activation function
+  void set_activation_function(std::unique_ptr<ActivationFunction> af) override;
 
- public:
+ protected:
   /**
-   * @brief Set the activation function (takes ownership)
-   *
-   * @param af Unique pointer to the activation function
+   * @brief Construct a ChamberElastanceInductor with custom block type and
+   * parameters. Used by derived classes.
    */
-  void set_activation_function(std::unique_ptr<ActivationFunction> af) override;
+  ChamberElastanceInductor(
+      int id, Model* model, BlockType block_type,
+      std::vector<std::pair<std::string, InputParameter>> params)
+      : Block(id, model, block_type, BlockClass::chamber, params) {}
+
+  double Elas = 0.0;   ///< Current chamber elastance
+  double Vrest = 0.0;   ///< Current rest volume
+  double act_ = 0.0;   ///< Last computed activation
+  std::unique_ptr<ActivationFunction> activation_func_;
 
- private:
   /**
-   * @brief Update the elastance functions which depend on time
-   *
-   * @param parameters Parameters of the model
+   * @brief Compute elastance and rest volume from activation and parameters.
    */
-  void get_elastance_values(std::vector<double>& parameters);
+  virtual void get_elastance_values(std::vector<double>& parameters);
 };
 
 #endif  // SVZERODSOLVER_MODEL_CHAMBERELASTANCEINDUCTOR_HPP_

src/model/ChamberElastanceInductorExponential.cpp

@@ -0,0 +1,41 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
+// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
+#include "ChamberElastanceInductorExponential.h"
+
+#include "Model.h"
+
+void ChamberElastanceInductorExponential::get_elastance_values(
+    std::vector<double>& parameters) {
+  double Emax = parameters[global_param_ids[ParamId::EMAX]];
+
+  act_ = activation_func_->compute(model->time);
+  Elas = act_ * Emax;
+  Vrest = parameters[global_param_ids[ExponentialParamId::VASO]];
+}
+
+void ChamberElastanceInductorExponential::update_time(
+    SparseSystem& system, std::vector<double>& parameters) {
+  get_elastance_values(parameters);
+
+  // Eq 0: F[0][4] = -E(t). C[0] is set in update_solution (nonlinear).
+  system.F.coeffRef(global_eqn_ids[0], global_var_ids[4]) = -Elas;
+}
+
+void ChamberElastanceInductorExponential::update_solution(
+    SparseSystem& system, std::vector<double>& parameters,
+    const Eigen::Matrix<double, Eigen::Dynamic, 1>& y,
+    const Eigen::Matrix<double, Eigen::Dynamic, 1>& dy) {
+  double Kxp = parameters[global_param_ids[ExponentialParamId::KXP]];
+  double Kxv = parameters[global_param_ids[ExponentialParamId::KXV]];
+  double Vaso = parameters[global_param_ids[ExponentialParamId::VASO]];
+  double Emax = parameters[global_param_ids[ParamId::EMAX]];
+  double Vc = y[global_var_ids[4]];
+
+  double exp_term = std::exp(Kxv * (Vc - Vaso));
+  double psi = Kxp * (exp_term - 1.0);
+  double psi_d = Kxp * Kxv * exp_term;
+
+  system.C(global_eqn_ids[0]) = act_ * Emax * Vaso + psi * (act_ - 1.0);
+  system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[4]) =
+      psi_d * (act_ - 1.0);
+}