Piecewise Chamber and Valve

docs/references.bib

@@ -125,3 +125,18 @@ @article{caruel13
  volume = {13},
  year = {2013}
 }
+
+@article{Regazzoni2022,
+title = {A cardiac electromechanical model coupled with a lumped-parameter model for closed-loop blood circulation},
+journal = {Journal of Computational Physics},
+volume = {457},
+pages = {111083},
+year = {2022},
+issn = {0021-9991},
+doi = {https://doi.org/10.1016/j.jcp.2022.111083},
+url = {https://www.sciencedirect.com/science/article/pii/S0021999122001450},
+author = {F. Regazzoni and M. Salvador and P.C. Africa and M. Fedele and L. Dedè and A. Quarteroni},
+keywords = {Mathematical modeling, Multiphysics models, Multiscale models, Cardiac modeling, Cardiac electromechanics},
+abstract = {We propose a novel mathematical and numerical model for cardiac electromechanics, wherein biophysically detailed core models describe the different physical processes concurring to the cardiac function. The core models, once suitably approximated, are coupled by a computationally efficient strategy. Our model is based on: (1) the combination of implicit-explicit (IMEX) schemes to solve the different core cardiac models, (2) an Artificial Neural Network based model, that surrogates a biophysically detailed but computationally demanding microscale model of active force generation and (3) appropriate partitioned schemes to couple the different models in this multiphysics setting. We employ a flexible and scalable intergrid transfer operator, which allows to interpolate Finite Element functions between nested meshes and, possibly, among arbitrary Finite Element spaces for the different core models. Our core 3D electromechanical model of the left ventricle is coupled with a closed-loop 0D model of the vascular network (and the other cardiac chambers) by an approach that is energy preserving. More precisely, we derive a balance law for the mechanical energy of the whole circulatory network. This provides a quantitative insight into the energy utilization, dissipation and transfer among the different compartments of the cardiovascular network and during different stages of the heartbeat. On this ground, a new tool is proposed to validate some energy indicators adopted in the daily clinical practice. A further contribution of this paper is the proposition of a robust algorithm for the reconstruction of the stress-free reference configuration. This feature is fundamental to correctly initialize our electromechanical simulations. As a matter of fact, the geometry acquired from medical imaging typically refers to a configuration affected by residual internal stresses, whereas the elastodynamics equations that govern the mechanics core model are related to a stress-free configuration. To prove the biophysical accuracy of our computational model, we address different scenarios of clinical interest, namely by varying preload, afterload and contractility.}
+}
+

src/model/BlockType.h

@@ -28,7 +28,9 @@ enum class BlockType {
   closed_loop_heart_pulmonary = 12,
   valve_tanh = 13,
   chamber_elastance_inductor = 14,
-  chamber_sphere = 15
+  chamber_sphere = 15,
+  piecewise_cosine_chamber = 17,
+  piecewise_valve = 18
 };
 
 /**

src/model/CMakeLists.txt

@@ -28,6 +28,8 @@ set(CXXSRCS
   ResistiveJunction.cpp
   ValveTanh.cpp
   WindkesselBC.cpp
+  PiecewiseCosineChamber.cpp
+  PiecewiseValve.cpp
 )
 
 set(HDRS 
@@ -54,6 +56,8 @@ set(HDRS
   ResistiveJunction.h
   ValveTanh.h
   WindkesselBC.h
+  PiecewiseCosineChamber.h
+  PiecewiseValve.h
 )
 
 add_library(${lib} OBJECT ${CXXSRCS} ${HDRS})

src/model/Model.cpp

@@ -28,7 +28,9 @@ Model::Model() {
       {"RESISTANCE", block_factory<ResistanceBC>()},
       {"resistive_junction", block_factory<ResistiveJunction>()},
       {"ValveTanh", block_factory<ValveTanh>()},
-      {"ChamberElastanceInductor", block_factory<ChamberElastanceInductor>()}};
+      {"ChamberElastanceInductor", block_factory<ChamberElastanceInductor>()},
+      {"PiecewiseValve", block_factory<PiecewiseValve>()},
+      {"PiecewiseCosineChamber", block_factory<PiecewiseCosineChamber>()}};
 }
 
 Model::~Model() {}
@@ -116,7 +118,7 @@ std::string Model::get_block_name(int block_id) const {
 int Model::add_node(const std::vector<Block*>& inlet_eles,
                     const std::vector<Block*>& outlet_eles,
                     const std::string_view& name) {
-  // DEBUG_MSG("Adding node " << name);
+  DEBUG_MSG("Adding node " << name);
   auto node = std::shared_ptr<Node>(
       new Node(node_count, inlet_eles, outlet_eles, this));
   nodes.push_back(node);

src/model/Model.h

@@ -31,6 +31,8 @@
 #include "Node.h"
 #include "OpenLoopCoronaryBC.h"
 #include "Parameter.h"
+#include "PiecewiseCosineChamber.h"
+#include "PiecewiseValve.h"
 #include "PressureReferenceBC.h"
 #include "ResistanceBC.h"
 #include "ResistiveJunction.h"

src/model/OpenLoopCoronaryBC.h

@@ -73,7 +73,7 @@
  * Parameter sequence for constructing this block
  *
  * * `0` Ra: Small artery resistance
- * * `1` Ram: Microvascualar resistance
+ * * `1` Ram: Microvascular resistance
  * * `2` Rv: Venous resistance
  * * `3` Ca: Small artery capacitance
  * * `4` Cim: Intramyocardial capacitance

src/model/PiecewiseCosineChamber.cpp

@@ -0,0 +1,62 @@
+// 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<double>& 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<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;
+  system.C.coeffRef(global_eqn_ids[0]) =
+      Elas * parameters[global_param_ids[ParamId::VREST]];
+}
+
+void PiecewiseCosineChamber::get_elastance_values(
+    std::vector<double>& 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;
+}

src/model/PiecewiseCosineChamber.h

@@ -0,0 +1,200 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
+// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
+
+/**
+ * @file PiecewiseCosineChamber.h
+ * @brief model::PiecewiseCosineChamber source file
+ */
+
+#ifndef SVZERODSOLVER_MODEL_PIECEWISECOSINECHAMBER_HPP_
+#define SVZERODSOLVER_MODEL_PIECEWISECOSINECHAMBER_HPP_
+
+#include <math.h>
+
+#include "Block.h"
+#include "Model.h"
+#include "SparseSystem.h"
+#include "debug.h"
+
+/**
+ * @brief Cardiac chamber with piecewise elastance and inductor.
+ *
+ * Models a cardiac chamber as a time-varying capacitor (elastance with
+ * specified resting volumes) and an inductor. See \cite Regazzoni2022
+ *
+ * This chamber block can be connected to other blocks using junctions.
+ *
+ * \f[
+ * \begin{circuitikz}
+ * \draw node[left] {$Q_{in}$} [-latex] (0,0) -- (0.8,0);
+ * \draw (1,0) node[anchor=south]{$P_{in}$}
+ * to[short, *-*] (3,0) node[anchor=south]{$P_{out}$}
+ * (1,0) to [vC, l=$E$, *-] (1,-1.5)
+ * node[ground]{};
+ * \draw (3.2,0) -- (4.0,0) node[right] {$Q_{out}$} [-latex];
+ * \end{circuitikz}
+ * \f]
+ *
+ * ### Governing equations
+ *
+ * \f[
+ * P_{in}-E(t)(V_c-V_{rest})=0
+ * \f]
+ *
+ * \f[
+ * P_{in}-P_{out}=0
+ * \f]
+ *
+ * \f[
+ * 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 & 0 & 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,
+ *
+ * \f[
+ * E_i(t) = E_i^{\text{pass}} + E_i^{\text{act,max}} \,
+ * \phi\!\left(t, t_C^i, t_R^i, T_C^i, T_R^i\right),
+ * \f]
+ *
+ * \f[
+ * \phi(t, t_C, t_R, T_C, T_R) =
+ * \begin{cases}
+ * \frac{1}{2}\left[1 - \cos\!\left(\frac{\pi}{T_C} \operatorname{mod}(t - t_C,
+ * T_{\mathrm{HB}})\right)\right], & \text{if } 0 \le \operatorname{mod}(t -
+ * t_C, T_{\mathrm{HB}}) < T_C, \\[1.2em] \frac{1}{2}\left[1 +
+ * \cos\!\left(\frac{\pi}{T_R} \operatorname{mod}(t - t_R,
+ * T_{\mathrm{HB}})\right)\right], & \text{if } 0 \le \operatorname{mod}(t -
+ * t_R, T_{\mathrm{HB}}) < T_R, \\[1.2em] 0, & \text{otherwise.} \end{cases} \f]
+ *
+ * ### Parameters
+ *
+ * Parameter sequence for constructing this block
+ *
+ * * `0` Emax: Maximum elastance
+ * * `1` Epass: Passive elastance
+ * * `2` Vrest: Rest diastolic volume
+ * * `3` contract_start: Contract start time
+ * * `4` relax_start: Relax start time
+ * * `5` contract_duration: Contract duration
+ * * `6` relax_duration: Relaxation duration
+ *
+ * ### Internal variables
+ *
+ * Names of internal variables in this block's output:
+ *
+ * * `Vc`: Chamber volume
+ *
+ */
+class PiecewiseCosineChamber : public Block {
+ public:
+  /**
+   * @brief Construct a new BloodVessel 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,
+              BlockClass::chamber,
+              {{"Emax", InputParameter()},
+               {"Epass", InputParameter()},
+               {"Vrest", InputParameter()},
+               {"contract_start", InputParameter()},
+               {"relax_start", InputParameter()},
+               {"contract_duration", InputParameter()},
+               {"relax_duration", InputParameter()}}) {}
+
+  /**
+   * @brief Local IDs of the parameters
+   *
+   */
+  enum ParamId {
+    EMAX = 0,
+    EPASS = 1,
+    VREST = 2,
+    CSTART = 3,
+    RSTART = 4,
+    CDUR = 5,
+    RDUR = 6
+  };
+
+  /**
+   * @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 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, 0};
+
+ private:
+  double Elas;  // Chamber Elastance
+
+  /**
+   * @brief Update the elastance functions which depend on time
+   *
+   * @param parameters Parameters of the model
+   */
+  void get_elastance_values(std::vector<double>& parameters);
+};
+
+#endif  // SVZERODSOLVER_MODEL_PIECEWISECOSINECHAMBER_HPP_