CRL Vessel

src/model/BlockType.h

@@ -29,6 +29,7 @@ enum class BlockType {
   valve_tanh = 13,
   chamber_elastance_inductor = 14,
   chamber_sphere = 15,
+  blood_vessel_CRL = 16,
   piecewise_cosine_chamber = 17,
   piecewise_valve = 18
 };

src/model/BloodVesselCRL.cpp

@@ -0,0 +1,87 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
+// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
+
+#include "BloodVesselCRL.h"
+
+void BloodVesselCRL::setup_dofs(DOFHandler& dofhandler) {
+  Block::setup_dofs_(dofhandler, 2, {});
+}
+
+void BloodVesselCRL::update_constant(SparseSystem& system,
+                                     std::vector<double>& parameters) {
+  // Get parameters
+  double capacitance = parameters[global_param_ids[ParamId::CAPACITANCE]];
+  double inductance = parameters[global_param_ids[ParamId::INDUCTANCE]];
+  double resistance = parameters[global_param_ids[ParamId::RESISTANCE]];
+
+  // Set element contributions
+  // coeffRef args are the indices (i,j) of the matrix
+  // global_eqn_ids: number of rows in the matrix, set in setup_dofs
+  // global_var_ids: number of columns, organized as pressure and flow of all
+  // inlets and then all outlets of the block
+  system.E.coeffRef(global_eqn_ids[0], global_var_ids[3]) = -inductance;
+  system.E.coeffRef(global_eqn_ids[1], global_var_ids[0]) = -capacitance;
+  system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = 1.0;
+  system.F.coeffRef(global_eqn_ids[0], global_var_ids[3]) = -resistance;
+  system.F.coeffRef(global_eqn_ids[0], global_var_ids[2]) = -1.0;
+  system.F.coeffRef(global_eqn_ids[1], global_var_ids[1]) = 1.0;
+  system.F.coeffRef(global_eqn_ids[1], global_var_ids[3]) = -1.0;
+}
+
+void BloodVesselCRL::update_solution(
+    SparseSystem& system, std::vector<double>& parameters,
+    const Eigen::Matrix<double, Eigen::Dynamic, 1>& y,
+    const Eigen::Matrix<double, Eigen::Dynamic, 1>& dy) {
+  // Get parameters
+  double capacitance = parameters[global_param_ids[ParamId::CAPACITANCE]];
+  double stenosis_coeff =
+      parameters[global_param_ids[ParamId::STENOSIS_COEFFICIENT]];
+  double q_out = y[global_var_ids[3]];
+  double dq_out = dy[global_var_ids[3]];
+  double stenosis_resistance = stenosis_coeff * fabs(q_out);
+
+  // Set element contributions
+  system.C(global_eqn_ids[0]) = stenosis_resistance * -q_out;
+
+  double sgn_q_out = (0.0 < q_out) - (q_out < 0.0);
+  system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[1]) =
+      stenosis_coeff * sgn_q_out * -2.0 * q_out;
+}
+
+void BloodVesselCRL::update_gradient(
+    Eigen::SparseMatrix<double>& jacobian,
+    Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
+    Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha, std::vector<double>& y,
+    std::vector<double>& dy) {
+  auto y0 = y[global_var_ids[0]];
+  auto y1 = y[global_var_ids[1]];
+  auto y2 = y[global_var_ids[2]];
+  auto y3 = y[global_var_ids[3]];
+
+  auto dy0 = dy[global_var_ids[0]];
+  auto dy1 = dy[global_var_ids[1]];
+  auto dy3 = dy[global_var_ids[3]];
+
+  auto resistance = alpha[global_param_ids[ParamId::RESISTANCE]];
+  auto capacitance = alpha[global_param_ids[ParamId::CAPACITANCE]];
+  auto inductance = alpha[global_param_ids[ParamId::INDUCTANCE]];
+  double stenosis_coeff = 0.0;
+
+  if (global_param_ids.size() > 3) {
+    stenosis_coeff = alpha[global_param_ids[ParamId::STENOSIS_COEFFICIENT]];
+  }
+  auto stenosis_resistance = stenosis_coeff * fabs(y3);
+
+  jacobian.coeffRef(global_eqn_ids[0], global_param_ids[0]) = -y3;
+  jacobian.coeffRef(global_eqn_ids[0], global_param_ids[2]) = -dy3;
+
+  if (global_param_ids.size() > 3) {
+    jacobian.coeffRef(global_eqn_ids[0], global_param_ids[3]) = -fabs(y3) * y3;
+  }
+
+  jacobian.coeffRef(global_eqn_ids[1], global_param_ids[1]) = -dy0;
+
+  residual(global_eqn_ids[0]) =
+      y0 - (resistance + stenosis_resistance) * y3 - y2 - inductance * dy3;
+  residual(global_eqn_ids[1]) = y1 - y3 - capacitance * dy0;
+}

src/model/BloodVesselCRL.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 BloodVesselCRL.h
+ * @brief model::BloodVesselCRL source file
+ */
+#ifndef SVZERODSOLVER_MODEL_BLOODVESSELCRL_HPP_
+#define SVZERODSOLVER_MODEL_BLOODVESSELCRL_HPP_
+
+#include <math.h>
+
+#include "Block.h"
+#include "SparseSystem.h"
+
+/**
+ * @brief Capacitor-resistor-inductor blood vessel with optional stenosis
+ *
+ * Models the mechanical behavior of a bloodvesselCRL with optional stenosis.
+ *
+ * \f[
+ * \begin{circuitikz} \draw
+ * node[left] {$Q_{in}$} [-latex] (0,0) -- (0.8,0);
+ * \draw (1,0) node[anchor=south]{$P_{in}$}
+ * to [R, l=$R$, *-] (3,0)
+ * to [R, l=$S$, -] (5,0)
+ * (5,0) to [L, l=$L$, -*] (7,0)
+ * node[anchor=south]{$P_{out}$}
+ * (1,0) to [C, l=$C$, -] (1,-1.5)
+ * node[ground]{};
+ * \draw [-latex] (7.2,0) -- (8,0) node[right] {$Q_{out}$};
+ * \end{circuitikz}
+ * \f]
+ *
+ * ### Governing equations
+ *
+ * \f[
+ * P_\text{in}-P_\text{out} - (R + S|Q_\text{out}|) Q_\text{out}-L
+ * \dot{Q}_\text{out}=0 \f]
+ *
+ * \f[
+ * Q_\text{in}-Q_\text{out} - C \dot{P}_\text{in}=0 \f]
+ *
+ * ### Local contributions
+ *
+ * \f[
+ * \mathbf{y}^{e}=\left[\begin{array}{llll}P_{i n} & Q_{in} &
+ * P_{out} & Q_{out}\end{array}\right]^\text{T} \f]
+ *
+ * \f[
+ * \mathbf{F}^{e}=\left[\begin{array}{cccc}
+ * 1 & 0 & -1 &  -R \\
+ * 0 &  1 &  0 & -1
+ * \end{array}\right]
+ * \f]
+ *
+ * \f[
+ * \mathbf{E}^{e}=\left[\begin{array}{cccc}
+ *  0 &  0 & 0 & -L \\
+ * -C & 0 & 0 &  0
+ * \end{array}\right]
+ * \f]
+ *
+ * \f[
+ * \mathbf{c}^{e} = S|Q_\text{out}|
+ * \left[\begin{array}{c}
+ * -Q_\text{out} \\
+ * 0
+ * \end{array}\right]
+ * \f]
+ *
+ * \f[
+ * \left(\frac{\partial\mathbf{c}}{\partial\mathbf{y}}\right)^{e} =
+ *  S \text{sgn} (Q_\text{in})
+ * \left[\begin{array}{cccc}
+ * 0 & -2Q_\text{out}        & 0 & 0 \\
+ * 0 & 0 & 0 & 0
+ * \end{array}\right]
+ * \f]
+ *
+ * \f[
+ * \left(\frac{\partial\mathbf{c}}{\partial\dot{\mathbf{y}}}\right)^{e} =
+ *  S|Q_\text{out}|
+ * \left[\begin{array}{cccc}
+ * 0 &  0 & 0 & 0 \\
+ * 0 & 0 & 0 & 0
+ * \end{array}\right]
+ * \f]
+ *
+ * with the stenosis resistance \f$ S=K_{t} \frac{\rho}{2
+ * A_{o}^{2}}\left(\frac{A_{o}}{A_{s}}-1\right)^{2} \f$.
+ * \f$R\f$, \f$C\f$, and \f$L\f$ refer to
+ * Poisieuille resistance, capacitance and inductance, respectively.
+ *
+ * ### Gradient
+ *
+ * Gradient of the equations with respect to the parameters:
+ *
+ * \f[
+ * \mathbf{J}^{e} = \left[\begin{array}{cccc}
+ * -y_3 & 0 & -\dot{y}_3 & -|y_3|y_3 \\
+ * 0 & 0 & -\dot{y}_0 & 0 \\
+ * \end{array}\right]
+ * \f]
+ *
+ * ### Parameters
+ *
+ * Parameter sequence for constructing this block
+ *
+ * * `0` Poiseuille resistance
+ * * `1` Capacitance
+ * * `2` Inductance
+ * * `3` Stenosis coefficient
+ *
+ */
+class BloodVesselCRL : public Block {
+ public:
+  /**
+   * @brief Local IDs of the parameters
+   *
+   */
+  enum ParamId {
+    RESISTANCE = 0,
+    CAPACITANCE = 1,
+    INDUCTANCE = 2,
+    STENOSIS_COEFFICIENT = 3,
+  };
+
+  /**
+   * @brief Construct a new BloodVesselCRL object
+   *
+   * @param id Global ID of the block
+   * @param model The model to which the block belongs
+   */
+  BloodVesselCRL(int id, Model* model)
+      : Block(id, model, BlockType::blood_vessel_CRL, BlockClass::vessel,
+              {{"R_poiseuille", InputParameter()},
+               {"C", InputParameter(true)},
+               {"L", InputParameter(true)},
+               {"stenosis_coefficient", InputParameter(true)}}) {}
+
+  /**
+   * @brief Set up the degrees of freedom (DOF) of the block
+   *
+   * Set \ref global_var_ids and \ref 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 solution-dependent contributions of the element in a
+   * sparse system
+   *
+   * @param system System to update contributions at
+   * @param parameters Parameters of the model
+   * @param y Current solution
+   * @param dy Current derivate of the solution
+   */
+  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 Set the gradient of the block contributions with respect to the
+   * parameters
+   *
+   * @param jacobian Jacobian with respect to the parameters
+   * @param alpha Current parameter vector
+   * @param residual Residual with respect to the parameters
+   * @param y Current solution
+   * @param dy Time-derivative of the current solution
+   */
+  void update_gradient(Eigen::SparseMatrix<double>& jacobian,
+                       Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
+                       Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha,
+                       std::vector<double>& y, std::vector<double>& 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{5, 3, 2};
+};
+
+#endif  // SVZERODSOLVER_MODEL_BLOODVESSELCRL_HPP_

src/model/BloodVesselJunction.h

@@ -9,6 +9,7 @@
 
 #include "Block.h"
 #include "BloodVessel.h"
+#include "BloodVesselCRL.h"
 #include "SparseSystem.h"
 
 /**

src/model/CMakeLists.txt

@@ -8,6 +8,7 @@ set(lib svzero_model_library)
 set(CXXSRCS 
   Block.cpp
   BloodVessel.cpp
+  BloodVesselCRL.cpp
   BloodVesselJunction.cpp
   ChamberSphere.cpp
   ChamberElastanceInductor.cpp
@@ -36,6 +37,7 @@ set(HDRS
   Block.h
   BlockType.h
   BloodVessel.h
+  BloodVesselCRL.h
   BloodVesselJunction.h
   ChamberSphere.h
   ChamberElastanceInductor.h

src/model/Model.cpp

@@ -29,6 +29,7 @@ Model::Model() {
       {"resistive_junction", block_factory<ResistiveJunction>()},
       {"ValveTanh", block_factory<ValveTanh>()},
       {"ChamberElastanceInductor", block_factory<ChamberElastanceInductor>()},
+      {"BloodVesselCRL", block_factory<BloodVesselCRL>()},
       {"PiecewiseValve", block_factory<PiecewiseValve>()},
       {"PiecewiseCosineChamber", block_factory<PiecewiseCosineChamber>()}};
 }

src/model/Model.h

@@ -18,6 +18,7 @@
 #include "Block.h"
 #include "BlockFactory.h"
 #include "BloodVessel.h"
+#include "BloodVesselCRL.h"
 #include "BloodVesselJunction.h"
 #include "ChamberElastanceInductor.h"
 #include "ChamberSphere.h"

src/solve/csv_writer.cpp

@@ -46,7 +46,8 @@ std::string to_vessel_csv(const std::vector<double>& times,
     if (dynamic_cast<const BloodVessel*>(block) == nullptr &&
         dynamic_cast<const ChamberSphere*>(block) == nullptr &&
         dynamic_cast<const PiecewiseCosineChamber*>(block) == nullptr &&
-        dynamic_cast<const PiecewiseValve*>(block) == nullptr) {
+        dynamic_cast<const PiecewiseValve*>(block) == nullptr &&
+        dynamic_cast<const BloodVesselCRL*>(block) == nullptr) {
       continue;
     }
 

tests/cases/double_pulsatileFlow_CRL.json → tests/cases/pulsatileFlow_CRL.json

@@ -18,7 +18,7 @@
       "bc_type": "PRESSURE",
       "bc_values": {
         "P": [
-          0.01047121, 0.020937827, 0.03139526, 0.041838922, 0.052264232,
+          0.0, 0.01047121, 0.020937827, 0.03139526, 0.041838922, 0.052264232,
           0.062666617, 0.073041514, 0.083384373, 0.093690657, 0.103955845,
           0.114175435, 0.124344944, 0.13445991, 0.144515898, 0.154508497,
           0.164433323, 0.174286024, 0.184062276, 0.193757793, 0.203368322,
@@ -80,7 +80,7 @@
           -0.010471211, -5.89793e-10
         ],
         "t": [
-          0.020943951, 0.041887902, 0.062831853, 0.083775804, 0.104719755,
+          0.0, 0.020943951, 0.041887902, 0.062831853, 0.083775804, 0.104719755,
           0.125663706, 0.146607657, 0.167551608, 0.188495559, 0.20943951,
           0.230383461, 0.251327412, 0.272271363, 0.293215314, 0.314159265,
           0.335103216, 0.356047167, 0.376991118, 0.397935069, 0.41887902,
@@ -148,7 +148,7 @@
       "bc_type": "FLOW",
       "bc_values": {
         "Q": [
-          0.99912283, 0.996492859, 0.992114701, 0.985996037, 0.978147601,
+          1.0, 0.99912283, 0.996492859, 0.992114701, 0.985996037, 0.978147601,
           0.968583161, 0.957319498, 0.94437637, 0.929776486, 0.913545458,
           0.89571176, 0.87630668, 0.85536426, 0.832921241, 0.809016994,
           0.783693457, 0.756995056, 0.728968628, 0.699663341, 0.669130606,
@@ -207,10 +207,10 @@
           0.756995054, 0.783693456, 0.809016993, 0.832921239, 0.855364259,
           0.876306679, 0.895711759, 0.913545457, 0.929776485, 0.944376369,
           0.957319497, 0.968583161, 0.9781476, 0.985996037, 0.992114701,
-          0.996492859, 0.99912283, 1
+          0.996492859, 0.99912283, 1.0
         ],
         "t": [
-          0.020943951, 0.041887902, 0.062831853, 0.083775804, 0.104719755,
+          0.0, 0.020943951, 0.041887902, 0.062831853, 0.083775804, 0.104719755,
           0.125663706, 0.146607657, 0.167551608, 0.188495559, 0.20943951,
           0.230383461, 0.251327412, 0.272271363, 0.293215314, 0.314159265,
           0.335103216, 0.356047167, 0.376991118, 0.397935069, 0.41887902,
@@ -277,8 +277,10 @@
 
   "simulation_parameters": {
     "number_of_cardiac_cycles": 10,
-    "number_of_time_pts_per_cardiac_cycle": 100,
-    "output_variable_based": false
+    "number_of_time_pts_per_cardiac_cycle": 300,
+    "output_variable_based": false,
+    "output_all_cycles": false,
+    "cardiac_period": 6.283185306
   },
   "vessels": [
     {