Simulate an open-loop 0D model until it converges

src/model/Block.cpp

@@ -34,6 +34,8 @@
 
 std::string Block::get_name() { return this->model->get_block_name(this->id); }
 
+void Block::update_vessel_type(VesselType type) { vessel_type = type; }
+
 Block::~Block() {}
 
 void Block::setup_params_(const std::vector<int> &param_ids) {

src/model/Block.h

@@ -103,6 +103,7 @@ class Block {
   const Model *model;            ///< The model to which the block belongs
   const BlockType block_type;    ///< Type of this block
   const BlockClass block_class;  ///< Class of this block
+  VesselType vessel_type = VesselType::neither;  ///< Vessel type of this block
   const std::vector<std::pair<std::string, InputParameter>>
       input_params;  ///< Map from name to input parameter
 
@@ -179,6 +180,13 @@ class Block {
    */
   std::string get_name();
 
+  /**
+   * @brief Update vessel type of the block
+   *
+   * @param type Type of vessel
+   */
+  void update_vessel_type(VesselType type);
+
   /**
    * @brief Setup parameter IDs for the block
    * @param param_ids Global IDs of the block parameters

src/model/BlockType.h

@@ -71,4 +71,9 @@ enum class BlockClass {
   chamber = 6
 };
 
+/**
+ * @brief The types of vessel blocks supported.
+ */
+enum class VesselType { inlet = 0, outlet = 1, both = 2, neither = 3 };
+
 #endif

src/model/Model.cpp

@@ -288,3 +288,17 @@ TripletsContributions Model::get_num_triplets() const {
 
   return triplets_sum;
 }
+
+void Model::update_has_windkessel_bc(bool has_windkessel) {
+  has_windkessel_bc = has_windkessel;
+}
+
+void Model::update_largest_windkessel_time_constant(double time_constant) {
+  largest_windkessel_time_constant = time_constant;
+}
+
+bool Model::get_has_windkessel_bc() { return has_windkessel_bc; }
+
+double Model::get_largest_windkessel_time_constant() {
+  return largest_windkessel_time_constant;
+}

src/model/Model.h

@@ -305,6 +305,36 @@ class Model {
    */
   int get_num_blocks(bool internal = false) const;
 
+  /**
+   * @brief Specify if model has at least one Windkessel boundary condition
+   *
+   * @param has_windkessel Toggle if model has at least one Windkessel boundary
+   * condition
+   */
+  void update_has_windkessel_bc(bool has_windkessel);
+
+  /**
+   * @brief Update model with largest time constant among all Windkessel
+   * boundary conditions present in model
+   *
+   * @param time_constant Largest Windkessel time constant
+   */
+  void update_largest_windkessel_time_constant(double time_constant);
+
+  /**
+   * @brief Check if model has at least one Windkessel boundary condition
+   *
+   * @return bool True if model has at least one Windkessel boundary condition
+   */
+  bool get_has_windkessel_bc();
+
+  /**
+   * @brief Get largest Windkessel time constant in model
+   *
+   * @return double Largest Windkessel time constant of model
+   */
+  double get_largest_windkessel_time_constant();
+
  private:
   int block_count = 0;
   int node_count = 0;
@@ -333,6 +363,9 @@ class Model {
                          ///< to blocks to set up the linear system in
                          ///< `update_constant`, `update_time` and
                          ///< `update_solution`.
+
+  bool has_windkessel_bc = false;
+  double largest_windkessel_time_constant = 0.0;
 };
 
 #endif  // SVZERODSOLVER_MODEL_MODEL_HPP_

src/solve/SimulationParameters.cpp

@@ -188,6 +188,14 @@ SimulationParameters load_simulation_params(const nlohmann::json& config) {
         sim_config["number_of_time_pts_per_cardiac_cycle"];
     sim_params.sim_num_time_steps =
         (sim_params.sim_pts_per_cycle - 1) * sim_params.sim_num_cycles + 1;
+    sim_params.use_cycle_to_cycle_error =
+        sim_config.value("use_cycle_to_cycle_error", false);
+    if (sim_params.use_cycle_to_cycle_error) {
+      assert(sim_params.sim_num_cycles >=
+             2);  // need at least two cycles to compute cycle-to-cycle error
+      sim_params.sim_cycle_to_cycle_error =
+          sim_config.value("sim_cycle_to_cycle_percent_error", 1.0) / 100;
+    }
     sim_params.sim_external_step_size = 0.0;
 
   } else {
@@ -306,9 +314,15 @@ void create_vessels(
       const auto& vessel_bc_config = vessel_config["boundary_conditions"];
       if (vessel_bc_config.contains("inlet")) {
         connections.push_back({vessel_bc_config["inlet"], vessel_name});
+        if (vessel_bc_config.contains("outlet")) {
+          model.get_block(vessel_name)->update_vessel_type(VesselType::both);
+        } else {
+          model.get_block(vessel_name)->update_vessel_type(VesselType::inlet);
+        }
       }
       if (vessel_bc_config.contains("outlet")) {
         connections.push_back({vessel_name, vessel_bc_config["outlet"]});
+        model.get_block(vessel_name)->update_vessel_type(VesselType::outlet);
       }
     }
   }
@@ -329,6 +343,15 @@ void create_boundary_conditions(Model& model, const nlohmann::json& config,
     // Keep track of closed loop blocks
     Block* block = model.get_block(block_id);
 
+    if (block->block_type == BlockType::windkessel_bc) {
+      model.update_has_windkessel_bc(true);
+      double Rd = bc_values["Rd"];
+      double C = bc_values["C"];
+      double time_constant = Rd * C;
+      model.update_largest_windkessel_time_constant(std::max(
+          model.get_largest_windkessel_time_constant(), time_constant));
+    }
+
     if (block->block_type == BlockType::closed_loop_rcr_bc) {
       if (bc_values["closed_loop_outlet"] == true) {
         closed_loop_bcs.push_back(bc_name);

src/solve/SimulationParameters.h

@@ -53,9 +53,16 @@ struct SimulationParameters {
   double sim_time_step_size{0.0};  ///< Simulation time step size
   double sim_abs_tol{0.0};         ///< Absolute tolerance for simulation
 
-  int sim_num_cycles{0};      ///< Number of cardiac cycles to simulate
-  int sim_pts_per_cycle{0};   ///< Number of time steps per cardiac cycle
-  int sim_num_time_steps{0};  ///< Total number of time steps
+  int sim_num_cycles{0};     ///< Number of cardiac cycles to simulate
+  int sim_pts_per_cycle{0};  ///< Number of time steps per cardiac cycle
+  bool use_cycle_to_cycle_error{
+      false};  ///< If model does not have RCR boundary conditions, simulate
+               ///< model to convergence (based on cycle-to-cycle error of last
+               ///< two cardiac cycles); if it does, update number of cardiac
+               ///< cycles to simulate to be value estimated from equation 21 of
+               ///< Pfaller 2021
+  double sim_cycle_to_cycle_error{0};  ///< Cycle-to-cycle error
+  int sim_num_time_steps{0};           ///< Total number of time steps
   int sim_nliter{0};  ///< Maximum number of non-linear iterations in time
                       ///< integration
   double sim_rho_infty{0.0};  ///< Spectral radius of generalized-alpha

src/solve/Solver.cpp

@@ -14,6 +14,17 @@ Solver::Solver(const nlohmann::json& config) {
 
   DEBUG_MSG("Cardiac cycle period " << this->model->cardiac_cycle_period);
 
+  if (!simparams.sim_coupled && simparams.use_cycle_to_cycle_error &&
+      this->model->get_has_windkessel_bc()) {
+    simparams.sim_num_cycles =
+        int(ceil(-1 * this->model->get_largest_windkessel_time_constant() /
+                 this->model->cardiac_cycle_period *
+                 log(simparams.sim_cycle_to_cycle_error)));  // equation 21 of
+                                                             // Pfaller 2021
+    simparams.sim_num_time_steps =
+        (simparams.sim_pts_per_cycle - 1) * simparams.sim_num_cycles + 1;
+  }
+
   // Calculate time step size
   if (!simparams.sim_coupled) {
     simparams.sim_time_step_size = this->model->cardiac_cycle_period /
@@ -81,8 +92,29 @@ void Solver::run() {
     states.push_back(std::move(state));
   }
 
+  int num_time_pts_in_two_cycles;
+  std::vector<State> states_last_two_cycles;
+  int last_two_cycles_time_pt_counter = 0;
+  if (simparams.use_cycle_to_cycle_error) {
+    num_time_pts_in_two_cycles = 2 * (simparams.sim_pts_per_cycle - 1) + 1;
+    states_last_two_cycles =
+        std::vector<State>(num_time_pts_in_two_cycles, state);
+  }
   for (int i = 1; i < simparams.sim_num_time_steps; i++) {
+    if (simparams.use_cycle_to_cycle_error) {
+      if (i == simparams.sim_num_time_steps - num_time_pts_in_two_cycles + 1) {
+        // add first state
+        states_last_two_cycles[last_two_cycles_time_pt_counter] = state;
+        last_two_cycles_time_pt_counter +=
+            1;  // last_two_cycles_time_pt_counter becomes 1
+      }
+    }
     state = integrator.step(state, time);
+    if (simparams.use_cycle_to_cycle_error &&
+        last_two_cycles_time_pt_counter > 0) {
+      states_last_two_cycles[last_two_cycles_time_pt_counter] = state;
+      last_two_cycles_time_pt_counter += 1;
+    }
     interval_counter += 1;
     time = simparams.sim_time_step_size * double(i);
 
@@ -96,6 +128,68 @@ void Solver::run() {
     }
   }
 
+  if (simparams.use_cycle_to_cycle_error) {
+    std::vector<std::pair<int, int>> vessel_caps_dof_indices =
+        get_vessel_caps_dof_indices();
+
+    if (!(this->model->get_has_windkessel_bc())) {
+      assert(last_two_cycles_time_pt_counter == num_time_pts_in_two_cycles);
+      double converged = check_vessel_cap_convergence(states_last_two_cycles,
+                                                      vessel_caps_dof_indices);
+      int extra_num_cycles = 0;
+
+      while (!converged) {
+        std::rotate(
+            states_last_two_cycles.begin(),
+            states_last_two_cycles.begin() + simparams.sim_pts_per_cycle - 1,
+            states_last_two_cycles.end());
+        last_two_cycles_time_pt_counter = simparams.sim_pts_per_cycle;
+        for (size_t i = 1; i < simparams.sim_pts_per_cycle; i++) {
+          state = integrator.step(state, time);
+          states_last_two_cycles[last_two_cycles_time_pt_counter] = state;
+          last_two_cycles_time_pt_counter += 1;
+          interval_counter += 1;
+          time = simparams.sim_time_step_size * double(i);
+
+          if ((interval_counter == simparams.output_interval) ||
+              (!simparams.output_all_cycles && (i == start_last_cycle))) {
+            if (simparams.output_all_cycles || (i >= start_last_cycle)) {
+              times.push_back(time);
+              states.push_back(std::move(state));
+            }
+            interval_counter = 0;
+          }
+        }
+        extra_num_cycles++;
+        converged = check_vessel_cap_convergence(states_last_two_cycles,
+                                                 vessel_caps_dof_indices);
+        assert(last_two_cycles_time_pt_counter == num_time_pts_in_two_cycles);
+      }
+      std::cout << "Ran simulation for " << extra_num_cycles
+                << " more cycles to converge flow and pressures at caps"
+                << std::endl;
+    } else {
+      for (const std::pair<int, int>& dof_indices : vessel_caps_dof_indices) {
+        std::pair<double, double> cycle_to_cycle_errors_in_flow_and_pressure =
+            get_cycle_to_cycle_errors_in_flow_and_pressure(
+                states_last_two_cycles, dof_indices);
+        double cycle_to_cycle_error_flow =
+            cycle_to_cycle_errors_in_flow_and_pressure.first;
+        double cycle_to_cycle_error_pressure =
+            cycle_to_cycle_errors_in_flow_and_pressure.second;
+        std::cout << "Percent error between last two simulated cardiac cycles "
+                     "for dof index "
+                  << dof_indices.first
+                  << " (mean flow)    : " << cycle_to_cycle_error_flow * 100.0
+                  << std::endl;
+        std::cout << "Percent error between last two simulated cardiac cycles "
+                     "for dof index "
+                  << dof_indices.second << " (mean pressure): "
+                  << cycle_to_cycle_error_pressure * 100.0 << std::endl;
+      }
+    }
+  }
+
   DEBUG_MSG("Avg. number of nonlinear iterations per time step: "
             << integrator.avg_nonlin_iter());
 
@@ -108,6 +202,94 @@ void Solver::run() {
   }
 }
 
+std::vector<std::pair<int, int>> Solver::get_vessel_caps_dof_indices() {
+  std::vector<std::pair<int, int>> vessel_caps_dof_indices;
+
+  for (size_t i = 0; i < this->model->get_num_blocks(); i++) {
+    auto block = this->model->get_block(i);
+
+    if (block->block_class == BlockClass::vessel) {
+      if ((block->vessel_type == VesselType::inlet) ||
+          (block->vessel_type == VesselType::both)) {
+        int inflow_dof = block->inlet_nodes[0]->flow_dof;
+        int inpres_dof = block->inlet_nodes[0]->pres_dof;
+        std::pair<int, int> dofs{inflow_dof, inpres_dof};
+        vessel_caps_dof_indices.push_back(dofs);
+      } else if ((block->vessel_type == VesselType::outlet) ||
+                 (block->vessel_type == VesselType::both)) {
+        int outflow_dof = block->outlet_nodes[0]->flow_dof;
+        int outpres_dof = block->outlet_nodes[0]->pres_dof;
+        std::pair<int, int> dofs{outflow_dof, outpres_dof};
+        vessel_caps_dof_indices.push_back(dofs);
+      }
+    }
+  }
+
+  return vessel_caps_dof_indices;
+}
+
+bool Solver::check_vessel_cap_convergence(
+    const std::vector<State>& states_last_two_cycles,
+    const std::vector<std::pair<int, int>>& vessel_caps_dof_indices) {
+  double converged = true;
+  for (const std::pair<int, int>& dof_indices : vessel_caps_dof_indices) {
+    std::pair<double, double> cycle_to_cycle_errors_in_flow_and_pressure =
+        get_cycle_to_cycle_errors_in_flow_and_pressure(states_last_two_cycles,
+                                                       dof_indices);
+    double cycle_to_cycle_error_flow =
+        cycle_to_cycle_errors_in_flow_and_pressure.first;
+    double cycle_to_cycle_error_pressure =
+        cycle_to_cycle_errors_in_flow_and_pressure.second;
+
+    if (cycle_to_cycle_error_flow > simparams.sim_cycle_to_cycle_error ||
+        cycle_to_cycle_error_pressure > simparams.sim_cycle_to_cycle_error) {
+      converged = false;
+      break;
+    }
+  }
+
+  return converged;
+}
+
+std::pair<double, double>
+Solver::get_cycle_to_cycle_errors_in_flow_and_pressure(
+    const std::vector<State>& states_last_two_cycles,
+    const std::pair<int, int>& dof_indices) {
+  double mean_flow_second_to_last_cycle = 0.0;
+  double mean_pressure_second_to_last_cycle = 0.0;
+  double mean_flow_last_cycle = 0.0;
+  double mean_pressure_last_cycle = 0.0;
+
+  for (size_t i = 0; i < simparams.sim_pts_per_cycle; i++) {
+    mean_flow_second_to_last_cycle +=
+        states_last_two_cycles[i].y[dof_indices.first];
+    mean_pressure_second_to_last_cycle +=
+        states_last_two_cycles[i].y[dof_indices.second];
+    mean_flow_last_cycle +=
+        states_last_two_cycles[simparams.sim_pts_per_cycle - 1 + i]
+            .y[dof_indices.first];
+    mean_pressure_last_cycle +=
+        states_last_two_cycles[simparams.sim_pts_per_cycle - 1 + i]
+            .y[dof_indices.second];
+  }
+  mean_flow_second_to_last_cycle /= simparams.sim_pts_per_cycle;
+  mean_pressure_second_to_last_cycle /= simparams.sim_pts_per_cycle;
+  mean_flow_last_cycle /= simparams.sim_pts_per_cycle;
+  mean_pressure_last_cycle /= simparams.sim_pts_per_cycle;
+
+  double cycle_to_cycle_error_flow =
+      abs((mean_flow_last_cycle - mean_flow_second_to_last_cycle) /
+          mean_flow_second_to_last_cycle);
+  double cycle_to_cycle_error_pressure =
+      abs((mean_pressure_last_cycle - mean_pressure_second_to_last_cycle) /
+          mean_pressure_second_to_last_cycle);
+
+  std::pair<double, double> cycle_to_cycle_errors_in_flow_and_pressure{
+      cycle_to_cycle_error_flow, cycle_to_cycle_error_pressure};
+
+  return cycle_to_cycle_errors_in_flow_and_pressure;
+}
+
 std::vector<double> Solver::get_times() const { return times; }
 
 std::string Solver::get_full_result() const {