Fixed capacitance blood vessel

CMakeLists.txt

@@ -38,7 +38,7 @@ if(ENABLE_COVERAGE)
   # add coverage target
   add_custom_target(coverage
     # gather data
-    COMMAND ${LCOV} --directory . --capture --output-file coverage.info --ignore-errors gcov --exclude '/usr/include/*' --exclude '/usr/lib/*' --exclude '*/_deps/*'
+    COMMAND ${LCOV} --directory . --capture --output-file coverage.info --ignore-errors gcov --ignore-errors mismatch --exclude '/usr/include/*' --exclude '/usr/lib/*' --exclude '*/_deps/*'
     # generate report
     COMMAND ${GENHTML} --demangle-cpp -o coverage coverage.info
     WORKING_DIRECTORY ${CMAKE_BINARY_DIR})

src/model/BloodVesselFC.cpp

@@ -0,0 +1,152 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
+// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
+
+#include "BloodVesselFC.h"
+
+#include <cmath>
+
+#include "Model.h"
+
+void BloodVesselFC::setup_dofs(DOFHandler& dofhandler) {
+  Block::setup_dofs_(dofhandler, 2, {});
+}
+
+void BloodVesselFC::update_constant(SparseSystem& system,
+                                    std::vector<double>& parameters) {
+  // Get parameters
+  double inductance = parameters[global_param_ids[ParamId::INDUCTANCE]];
+  double resistance = parameters[global_param_ids[ParamId::RESISTANCE]];
+
+  // Get fixed capacitance from model by looking up this block's name
+  double capacitance = 0.0;
+  std::string block_name = "";
+  for (size_t i = 0; i < model->get_num_blocks(); i++) {
+    if (model->get_block(i) == this) {
+      block_name = model->get_block_name(i);
+      break;
+    }
+  }
+  if (!block_name.empty() && model->fixed_capacitance.count(block_name)) {
+    capacitance = model->fixed_capacitance.at(block_name);
+  }
+
+  // Set element contributions
+  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.E.coeffRef(global_eqn_ids[1], global_var_ids[1]) =
+      capacitance * resistance;
+  system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = 1.0;
+  system.F.coeffRef(global_eqn_ids[0], global_var_ids[1]) = -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 BloodVesselFC::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 fixed capacitance from model
+  std::string block_name = "";
+  for (size_t i = 0; i < model->get_num_blocks(); i++) {
+    if (model->get_block(i) == this) {
+      block_name = model->get_block_name(i);
+      break;
+    }
+  }
+  double capacitance = 0.0;
+  if (!block_name.empty() && model->fixed_capacitance.count(block_name)) {
+    capacitance = model->fixed_capacitance.at(block_name);
+  }
+
+  double stenosis_coeff = 0.0;
+  if (global_param_ids.size() >
+      static_cast<size_t>(ParamId::STENOSIS_COEFFICIENT)) {
+    stenosis_coeff =
+        parameters[global_param_ids[ParamId::STENOSIS_COEFFICIENT]];
+  }
+  double q_in = y[global_var_ids[1]];
+  double dq_in = dy[global_var_ids[1]];
+  double stenosis_resistance = stenosis_coeff * fabs(q_in);
+
+  // Set element contributions
+  system.C(global_eqn_ids[0]) = stenosis_resistance * -q_in;
+  system.C(global_eqn_ids[1]) = stenosis_resistance * 2.0 * capacitance * dq_in;
+
+  double sgn_q_in = (0.0 < q_in) - (q_in < 0.0);
+  system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[1]) =
+      stenosis_coeff * sgn_q_in * -2.0 * q_in;
+  system.dC_dy.coeffRef(global_eqn_ids[1], global_var_ids[1]) =
+      stenosis_coeff * sgn_q_in * 2.0 * capacitance * dq_in;
+
+  system.dC_dydot.coeffRef(global_eqn_ids[1], global_var_ids[1]) =
+      stenosis_resistance * 2.0 * capacitance;
+}
+
+void BloodVesselFC::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) {
+  // Check if all required observations are available (not NaN)
+  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]];
+
+  // Skip residual computation if any required observation is missing (NaN)
+  if (std::isnan(y0) || std::isnan(y1) || std::isnan(y2) || std::isnan(y3) ||
+      std::isnan(dy0) || std::isnan(dy1) || std::isnan(dy3)) {
+    return;
+  }
+
+  auto resistance = alpha[global_param_ids[ParamId::RESISTANCE]];
+  auto inductance = alpha[global_param_ids[ParamId::INDUCTANCE]];
+  double stenosis_coeff = 0.0;
+
+  if (global_param_ids.size() > 2) {
+    stenosis_coeff = alpha[global_param_ids[ParamId::STENOSIS_COEFFICIENT]];
+  }
+
+  // Get fixed capacitance from model
+  std::string block_name = "";
+  for (size_t i = 0; i < model->get_num_blocks(); i++) {
+    if (model->get_block(i) == this) {
+      block_name = model->get_block_name(i);
+      break;
+    }
+  }
+  double capacitance = 0.0;
+  if (!block_name.empty() && model->fixed_capacitance.count(block_name)) {
+    capacitance = model->fixed_capacitance.at(block_name);
+  }
+
+  auto stenosis_resistance = stenosis_coeff * fabs(y1);
+
+  // Jacobian entries for R (param 0) - same as BloodVessel
+  jacobian.coeffRef(global_eqn_ids[0], global_param_ids[0]) = -y1;
+  jacobian.coeffRef(global_eqn_ids[1], global_param_ids[0]) = capacitance * dy1;
+
+  // Jacobian entries for L (param 1) - same as BloodVessel param 2
+  jacobian.coeffRef(global_eqn_ids[0], global_param_ids[1]) = -dy3;
+
+  // Note: No Jacobian entries for C since it's fixed!
+
+  if (global_param_ids.size() > 2) {
+    // Jacobian entries for stenosis (param 2) - same as BloodVessel param 3
+    jacobian.coeffRef(global_eqn_ids[0], global_param_ids[2]) = -fabs(y1) * y1;
+    jacobian.coeffRef(global_eqn_ids[1], global_param_ids[2]) =
+        2.0 * capacitance * fabs(y1) * dy1;
+  }
+
+  // Residual - same as BloodVessel but using fixed capacitance
+  residual(global_eqn_ids[0]) =
+      y0 - (resistance + stenosis_resistance) * y1 - y2 - inductance * dy3;
+  residual(global_eqn_ids[1]) =
+      (y1 - y3 - capacitance * dy0 +
+       capacitance * (resistance + 2.0 * stenosis_resistance) * dy1);
+}

src/model/BloodVesselFC.h

@@ -0,0 +1,108 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
+// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
+/**
+ * @file BloodVesselFC.h
+ * @brief model::BloodVesselFC source file - Blood vessel with fixed capacitance
+ */
+#ifndef SVZERODSOLVER_MODEL_BLOODVESSELFC_HPP_
+#define SVZERODSOLVER_MODEL_BLOODVESSELFC_HPP_
+
+#include <math.h>
+
+#include "Block.h"
+#include "SparseSystem.h"
+
+/**
+ * @brief Resistor-capacitor-inductor blood vessel with fixed capacitance
+ *
+ * Same as BloodVessel, but capacitance is not an optimization parameter.
+ * Instead, capacitance is read from the model's fixed_capacitance map.
+ * This is used in the calibrator when capacitance should not be optimized.
+ *
+ * ### Parameters
+ *
+ * Parameter sequence for constructing this block (note: no capacitance)
+ *
+ * * `0` Poiseuille resistance
+ * * `1` Inductance
+ * * `2` Stenosis coefficient
+ *
+ */
+class BloodVesselFC : public Block {
+ public:
+  /**
+   * @brief Local IDs of the parameters
+   *
+   */
+  enum ParamId {
+    RESISTANCE = 0,
+    INDUCTANCE = 1,
+    STENOSIS_COEFFICIENT = 2,
+  };
+
+  /**
+   * @brief Construct a new BloodVesselFC object
+   *
+   * @param id Global ID of the block
+   * @param model The model to which the block belongs
+   */
+  BloodVesselFC(int id, Model* model)
+      : Block(id, model, BlockType::blood_vessel, BlockClass::vessel,
+              {{"R_poiseuille", InputParameter()},
+               {"L", InputParameter(true)},
+               {"stenosis_coefficient", InputParameter(true)}}) {}
+
+  /**
+   * @brief Set up the degrees of freedom (DOF) of the block
+   *
+   * @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_BLOODVESSELFC_HPP_

src/model/CMakeLists.txt

@@ -10,6 +10,7 @@ set(CXXSRCS
   Block.cpp
   BloodVessel.cpp
   BloodVesselCRL.cpp
+  BloodVesselFC.cpp
   BloodVesselJunction.cpp
   ChamberSphere.cpp
   ChamberElastanceInductor.cpp
@@ -41,6 +42,7 @@ set(HDRS
   BlockType.h
   BloodVessel.h
   BloodVesselCRL.h
+  BloodVesselFC.h
   BloodVesselJunction.h
   ChamberSphere.h
   ChamberElastanceInductor.h

src/model/Model.cpp

@@ -31,6 +31,7 @@ Model::Model() {
       {"ValveTanh", block_factory<ValveTanh>()},
       {"ChamberElastanceInductor", block_factory<ChamberElastanceInductor>()},
       {"BloodVesselCRL", block_factory<BloodVesselCRL>()},
+      {"BloodVesselFC", block_factory<BloodVesselFC>()},
       {"PiecewiseValve", block_factory<PiecewiseValve>()},
       {"LinearElastanceChamber", block_factory<LinearElastanceChamber>()}};
 }

src/model/Model.h

@@ -19,6 +19,7 @@
 #include "BlockFactory.h"
 #include "BloodVessel.h"
 #include "BloodVesselCRL.h"
+#include "BloodVesselFC.h"
 #include "BloodVesselJunction.h"
 #include "ChamberElastanceInductor.h"
 #include "ChamberSphere.h"
@@ -72,6 +73,9 @@ class Model {
   double cardiac_cycle_period = -1.0;  ///< Cardiac cycle period
   double time = 0.0;                   ///< Current time
 
+  /// Fixed capacitance values for BloodVesselFC blocks (vessel_name -> C)
+  std::map<std::string, double> fixed_capacitance;
+
   /**
    * @brief Create a new block
    *