Fix OpenLoopCoronary equations (SimVascular#125)

.gitignore

@@ -28,7 +28,7 @@ Release*/
 Debug*/
 externals/
 *.so
-build
+build*/
 
 # IDE
 .cproject

CMakeLists.txt

@@ -148,6 +148,7 @@ add_executable(svzerodsolver applications/svzerodsolver.cpp
 )
 
 add_executable(svzerodcalibrator applications/svzerodcalibrator.cpp
+ $<TARGET_OBJECTS:svzero_algebra_library> 
  $<TARGET_OBJECTS:svzero_optimize_library> 
  $<TARGET_OBJECTS:svzero_model_library> 
 )

applications/svzerodsolver.cpp

@@ -61,6 +61,7 @@ int main(int argc, char* argv[]) {
  }
 
  std::string input_file_name = argv[1];
+ std::string output_file_path;
  std::string output_file_name;
 
  if (argc == 3) {
@@ -73,12 +74,15 @@ int main(int argc, char* argv[]) {
  if (end_of_path == std::string::npos) {
  end_of_path = input_file_name.rfind("\\"); // For Windows paths (?)
 
+ // If <path to .json> is still not found, use current directory
  if (end_of_path == std::string::npos) {
- throw std::runtime_error("Error: No output file path provided. Tried to create a default output file but could not find the simulation directory from the input JSON file path.");
+ output_file_path = ".";
  }
+ } else {
+ output_file_path = input_file_name.substr(0, end_of_path);
  }
 
- output_file_name = input_file_name.substr(0, end_of_path) + "/output.csv";
+ output_file_name = output_file_path + "/output.csv";
  std::cout << "[svzerodsolver] Output will be written to '" << output_file_name << "'." << std::endl;;
  }
 

src/model/Block.cpp

@@ -72,6 +72,9 @@ void Block::setup_dofs(DOFHandler &dofhandler) {}
 
 void Block::setup_model_dependent_params() {}
 
+void Block::setup_initial_state_dependent_params(
+ State initial_state, std::vector<double> &parameters) {}
+
 void Block::update_constant(SparseSystem &system,
  std::vector<double> &parameters) {}
 

src/model/Block.h

@@ -43,6 +43,7 @@
 #include "DOFHandler.h"
 #include "Parameter.h"
 #include "SparseSystem.h"
+#include "State.h"
 
 /**
  * @brief The number of triplets that the element contributes
@@ -227,6 +228,15 @@ class Block {
  */
  virtual void setup_model_dependent_params();
 
+ /**
+ * @brief Setup parameters that depend on the initial state
+ *
+ * @param initial_state The initial state of the system
+ * @param parameters The parameter values vector (at time 0)
+ */
+ virtual void setup_initial_state_dependent_params(
+ State initial_state, std::vector<double> &parameters);
+
  /**
  * @brief Update the constant contributions of the element in a sparse system
  *

src/model/Model.cpp

@@ -254,6 +254,7 @@ void Model::to_steady() {
  param.to_steady();
  }
 
+ // Special handling for time-varying capacitance
  for (size_t i = 0; i < get_num_blocks(true); i++) {
  get_block(i)->steady = true;
  if ((block_types[i] == BlockType::windkessel_bc) ||
@@ -289,6 +290,14 @@ TripletsContributions Model::get_num_triplets() const {
  return triplets_sum;
 }
 
+void Model::setup_initial_state_dependent_parameters(State initial_state) {
+ DEBUG_MSG("Setup initial state dependent parameters");
+ for (auto &block : blocks) {
+ block->setup_initial_state_dependent_params(initial_state,
+ parameter_values);
+ }
+}
+
 void Model::update_has_windkessel_bc(bool has_windkessel) {
  has_windkessel_bc = has_windkessel;
 }

src/model/Model.h

@@ -60,6 +60,7 @@
 #include "PressureReferenceBC.h"
 #include "ResistanceBC.h"
 #include "ResistiveJunction.h"
+#include "State.h"
 #include "ValveTanh.h"
 #include "WindkesselBC.h"
 #include "debug.h"
@@ -334,6 +335,12 @@ class Model {
  * @return double Largest Windkessel time constant of model
  */
  double get_largest_windkessel_time_constant();
+ /**
+ * @brief Setup model parameters that depend on the initial state
+ *
+ * @param initial_state The initial state vector
+ */
+ void setup_initial_state_dependent_parameters(State initial_state);
 
  private:
  int block_count = 0;

src/model/OpenLoopCoronaryBC.cpp

@@ -43,10 +43,10 @@ void OpenLoopCoronaryBC::update_constant(SparseSystem &system,
  auto Cim = parameters[global_param_ids[4]];
 
  if (steady) {
- // Different assmembly for steady block to avoid singular system
- // and solve for the internal variable V_im inherently
- system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = -Cim;
- system.F.coeffRef(global_eqn_ids[0], global_var_ids[1]) = Cim * (Ra + Ram);
+ // Different equations for steady initial condition
+ // Equations are:
+ // -P_in + (Ra+Ram+Rv)Q_in + Pv = 0
+ // V_im = 0
  system.F.coeffRef(global_eqn_ids[0], global_var_ids[2]) = 1.0;
  system.F.coeffRef(global_eqn_ids[1], global_var_ids[0]) = -1.0;
  system.F.coeffRef(global_eqn_ids[1], global_var_ids[1]) = Ra + Ram + Rv;
@@ -74,10 +74,32 @@ void OpenLoopCoronaryBC::update_time(SparseSystem &system,
  auto Pv = parameters[global_param_ids[6]];
 
  if (steady) {
- system.C(global_eqn_ids[0]) = -Cim * Pim;
  system.C(global_eqn_ids[1]) = Pv;
  } else {
- system.C(global_eqn_ids[0]) = Cim * (-Pim + Pv);
- system.C(global_eqn_ids[1]) = -Cim * (Rv + Ram) * Pim + Ram * Cim * Pv;
+ system.C(global_eqn_ids[0]) =
+ Cim * (-Pim + Pv + this->Pim_0 - this->P_Cim_0);
+ system.C(global_eqn_ids[1]) =
+ (Ram * Cim * Pv) -
+ Cim * (Rv + Ram) * (Pim + this->P_Cim_0 - this->Pim_0);
  }
 }
+
+void OpenLoopCoronaryBC::setup_initial_state_dependent_params(
+ State initial_state, std::vector<double> &parameters) {
+ auto P_in = initial_state.y[global_var_ids[0]];
+ auto Q_in = initial_state.y[global_var_ids[1]];
+ auto P_in_dot = initial_state.ydot[global_var_ids[0]];
+ auto Q_in_dot = initial_state.ydot[global_var_ids[1]];
+ auto Ra = parameters[global_param_ids[0]];
+ auto Ram = parameters[global_param_ids[1]];
+ auto Ca = parameters[global_param_ids[3]];
+ // Pressure proximal to Ca and distal to Ra
+ auto P_Ca = P_in - Ra * Q_in;
+ auto P_Ca_dot = P_in_dot - Ra * Q_in_dot;
+ // Flow into Ram (inflow minus flow into Ca)
+ auto Q_am = Q_in - Ca * P_Ca_dot;
+ // Pressure proximal to Cim/Vim and distal to Ram
+ this->P_Cim_0 = P_Ca - Ram * Q_am;
+ // Initial intramyocardial pressure
+ this->Pim_0 = parameters[global_param_ids[5]];
+}

src/model/OpenLoopCoronaryBC.h

@@ -45,8 +45,9 @@
  * \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_a$, *-] (3,0)
+ * to [R, l=$R_a$, *-*] (3,0)
  * to [R, l=$R_{am}$, -] (5,0)
+ * node[anchor=south]{$P_{cim}$}
  * to [R, l=$R_v$, *-*] (7,0)
  * node[anchor=south]{$P_{v}$}
  * (5,0) to [C, l=$C_{im} \;V_{im}$, -*] (5,-1.5)
@@ -58,16 +59,18 @@
  *
  * ### Governing equations
  *
- * \f[
- * C_{i m} R_{v} Q_{in}-V_{i m}-C_{i m} P_{i m}+C_{i m} P_{v}-C_{i m} R_{v}
- * \frac{d V_{i m}}{d t}-C_{a} C_{i m} R_{v} \frac{d P_{in}}{d t}+R_{a} C_{a}
- * C_{i m} R_{v} \frac{d Q_{in}}{d t}+C_{a} C_{i m} R_{v} \frac{d P_{a}}{d
- * t}=0 \f]
+ * \f{eqnarray*}{
+ * &C_{i m} R_{v} Q_{in}-V_{i m}+C_{i m} \left(-P_{c i m}(0)+P_{i m}(0)-P_{i
+ * m}+P_{v}\right)-C_{i m} R_{v} \frac{d V_{i m}}{d t}-C_{a} C_{i m} R_{v}
+ * \frac{d P_{in}}{d t}+R_{a} C_{a} C_{i m} R_{v} \frac{d Q_{in}}{d t}\\ &+C_{a}
+ * C_{i m} R_{v} \frac{d P_{a}}{d t}=0
+ * \f}
  *
- * \f[
- * C_{i m} R_v P^{e}-C_{i m} R_{v} R_{a} Q^{e}-R_{v} V_{i m}^{e}-C_{i m} R_{v}
- * P_{i m}-C_{i m} R_{v} R_{a m} \frac{d V_{i m}^{e}}{d t}-R_{a m} V_{i
- * m}^{e}-C_{i m} R_{a m} P_{i m}+R_{a m} C_{i m} P_{v}=0 \f]
+ * \f{eqnarray*}{
+ * &C_{i m} R_v P_{in}-C_{i m} R_{v} R_{a} Q_{in}-R_{v} V_{i m}-C_{i
+ * m}\left(R_{v}+R_{a m}\right) \left(P_{c i m}(0)-P_{i m}(0)+P_{i
+ * m}\right)-C_{i m} R_{v} R_{a m} \frac{d V_{i m}}{d t} \\ &-R_{a m} V_{i
+ * m}+R_{a m} C_{i m} P_{v}=0 \f}
  *
  * ### Local contributions
  *
@@ -85,9 +88,10 @@
  * & -C_{i m} R_{v} R_{a} & -\left(R_{v}+R_{a m}\right)\end{array}\right] \f]
  *
  * \f[
- * \mathbf{c}^{e}=\left[\begin{array}{c}C_{i m}\left(-P_{i m}+P_{v}\right)+C_{a}
- * C_{i m} R_{v} \frac{d P_{a}}{d t} \\-C_{i m}\left(R_{v}+R_{a m}\right) P_{i
- * m}+R_{a m} C_{i m} P_{v}\end{array}\right] \f]
+ * \mathbf{c}^{e}=\left[\begin{array}{c}C_{i m}\left(-P_{i m}+P_{i m}(0)-P_{c i
+ * m}(0)+P_{v}\right)+C_{a} C_{i m} R_{v} \frac{d P_{a}}{d t} \\-C_{i
+ * m}\left(R_{v} + R_{a m}\right)\left(P_{cim}(0)-P_{i m}(0)+P_{i m}\right)+R_{a
+ * m} C_{i m} P_{v}\end{array}\right] \f]
  *
  * Assume \f$P_a=0\f$.
  *
@@ -143,6 +147,15 @@ class OpenLoopCoronaryBC : public Block {
  */
  void setup_dofs(DOFHandler &dofhandler);
 
+ /**
+ * @brief Setup parameters that depend on the initial state
+ *
+ * @param initial_state The initial state of the system
+ * @param parameters The parameter values vector (at time 0)
+ */
+ void setup_initial_state_dependent_params(State initial_state,
+ std::vector<double> &parameters);
+
  /**
  * @brief Update the constant contributions of the element in a sparse system
  *
@@ -167,6 +180,10 @@ class OpenLoopCoronaryBC : public Block {
  * (relevant for sparse memory reservation)
  */
  TripletsContributions num_triplets{5, 4, 0};
+
+ protected:
+ double P_Cim_0 = 0; ///< Pressure proximal to Cim/Vim at initial state
+ double Pim_0 = 0; ///< Pim at initial state
 };
 
 #endif // SVZERODSOLVER_MODEL_OPENLOOPCORONARYBC_HPP_

src/model/WindkesselBC.h

@@ -58,11 +58,11 @@
  * ### Governing equations
  *
  * \f[
- * R_{d} Q^{e}-P_{c}^{e}+P_{r e f}-R_{d} C \frac{d P_{c}^{e}}{d t}=0
+ * R_{d} Q_{in}-P_{c}+P_{r e f}-R_{d} C \frac{d P_{c}}{d t}=0
  * \f]
  *
  * \f[
- * P^{e}-P_{c}^{e}-R_{p} Q^{e}=0
+ * P_{in}-P_{c}-R_{p} Q_{in}=0
  * \f]
  *
  * ### Local contributions

src/solve/Solver.cpp

@@ -40,7 +40,7 @@ Solver::Solver(const nlohmann::json& config) {
 void Solver::run() {
  auto state = initial_state;
 
- // Create steady initial
+ // Create steady initial condition
  if (simparams.sim_steady_initial) {
  DEBUG_MSG("Calculate steady initial condition");
  double time_step_size_steady = this->model->cardiac_cycle_period / 10.0;
@@ -57,6 +57,10 @@ void Solver::run() {
  this->model->to_unsteady();
  }
 
+ // Use the initial condition (steady or user-provided) to set up parameters
+ // which depend on the initial condition
+ this->model->setup_initial_state_dependent_parameters(state);
+
  // Set-up integrator
  DEBUG_MSG("Setup time integration");
  Integrator integrator(this->model.get(), simparams.sim_time_step_size,