Add Impedance boundary condition

.gitignore

@@ -29,6 +29,7 @@ Debug*/
 externals/
 *.so
 build*/
+.npm-cache/
 
 # IDE
 .cproject

README.md

@@ -19,3 +19,37 @@ You can find more information under the following links:
 * [**Bug Reports**](https://github.com/simvascular/svZeroDSolver/issues)
 * [**Forum**](https://github.com/simvascular/svZeroDSolver/discussions)
 * [**About SimVascular**](https://simvascular.github.io)
+
+## Impedance BC (Olufsen periodic memory)
+
+`IMPEDANCE` is a boundary condition that enforces
+
+`P = Pd + z0 * Q + sum_{m=1..Nk-1} z[m] * Q_lag[m]`
+
+with one-cycle memory stored in a ring buffer and coupling-safe trial/accept
+semantics.
+
+Required JSON fields in `bc_values`:
+
+* `z`: time-domain impedance kernel array.
+
+Optional fields:
+
+* `Pd`: distal/reference pressure (default `0.0`).
+* `convolution_mode`: `exact` (default) or `truncated`.
+* `num_kernel_terms`: required when `convolution_mode` is `truncated`.
+
+Simulation timing requirements:
+
+* Non-coupled configs must set
+  `simulation_parameters.number_of_time_pts_per_cardiac_cycle = z.size() + 1`.
+* 3D-coupled configs must keep `simulation_parameters.number_of_time_pts = 2`
+  for the external interface step. The canonical coupled simulation parameters
+  are sufficient; if
+  `simulation_parameters.number_of_time_pts_per_cardiac_cycle` is omitted, the
+  impedance cycle discretization is inferred from
+  `simulation_parameters.cardiac_period / dt`.
+* In all cases, `simulation_parameters.cardiac_period / dt` must match
+  `z.size()`.
+
+For `truncated`, runtime cost is `O(num_kernel_terms)` per accepted 0D step.

src/algebra/Integrator.cpp

@@ -17,6 +17,7 @@ Integrator::Integrator(Model* model, double time_step_size, double rho,
   size = model->dofhandler.size();
   system = SparseSystem(size);
   this->time_step_size = time_step_size;
+  this->model->set_time_step_size(time_step_size);
   this->atol = atol;
   this->max_iter = max_iter;
 
@@ -27,7 +28,7 @@ Integrator::Integrator(Model* model, double time_step_size, double rho,
   system.reserve(model);
 }
 
-// Must declare default constructord and dedtructor
+// Must declare default constructors and dedtructor
 // because of Eigen.
 Integrator::Integrator() {}
 Integrator::~Integrator() {}
@@ -40,6 +41,7 @@ void Integrator::clean() {
 
 void Integrator::update_params(double time_step_size) {
   this->time_step_size = time_step_size;
+  model->set_time_step_size(time_step_size);
   y_coeff = gamma * time_step_size;
   y_coeff_jacobian = alpha_f * y_coeff;
   model->update_constant(system);
@@ -103,7 +105,7 @@ State Integrator::step(const State& old_state, double time) {
     // Count total number of nonlinear iterations
     n_nonlin_iter++;
   }
-
+  model->accept_timestep(new_state.y);
   return new_state;
 }
 

src/interface/interface.cpp

@@ -11,6 +11,12 @@
 int SolverInterface::problem_id_count_ = 0;
 std::map<int, SolverInterface*> SolverInterface::interface_list_;
 
+namespace {
+void commit_persistent_state_snapshot(SolverInterface* interface) {
+  interface->committed_block_states_ = interface->model_->get_persistent_states();
+}
+}  // namespace
+
 //-----------------
 // SolverInterface
 //-----------------
@@ -89,6 +95,7 @@ void initialize(std::string input_file_arg, int& problem_id, int& pts_per_cycle,
   auto simparams = load_simulation_params(config);
 
   auto model = std::shared_ptr<Model>(new Model());
+  model->impedance_points_per_cycle = simparams.sim_impedance_pts_per_cycle;
 
   load_simulation_model(config, *model.get());
   auto state = load_initial_condition(config, *model.get());
@@ -175,6 +182,8 @@ void initialize(std::string input_file_arg, int& problem_id, int& pts_per_cycle,
       state = integrator_steady.step(state, time_step_size_steady * double(i));
     }
     model_steady->to_unsteady();
+    // Reset dt-dependent persistent memory after steady initialization.
+    model_steady->clear_persistent_states();
   }
 
   interface->state_ = state;
@@ -187,6 +196,7 @@ void initialize(std::string input_file_arg, int& problem_id, int& pts_per_cycle,
   interface->integrator_ =
       Integrator(model.get(), interface->time_step_size_, interface->rho_infty_,
                  interface->absolute_tolerance_, interface->max_nliter_);
+  commit_persistent_state_snapshot(interface);
 
   DEBUG_MSG("[initialize] Done");
 }
@@ -321,6 +331,9 @@ void get_block_node_IDs(int problem_id, std::string block_name,
 /**
  * @brief Return the y state vector.
  *
+ * This also marks the current internal block persistent state as accepted for
+ * subsequent coupling rollback via update_state().
+ *
  * @param problem_id The ID used to identify the 0D problem.
  * @param y The state vector containing all state.y degrees-of-freedom.
  */
@@ -337,11 +350,15 @@ void return_y(int problem_id, std::vector<double>& y) {
   for (int i = 0; i < system_size; i++) {
     y[i] = state.y[i];
   }
+  commit_persistent_state_snapshot(interface);
 }
 
 /**
  * @brief Return the ydot state vector.
  *
+ * This also marks the current internal block persistent state as accepted for
+ * subsequent coupling rollback via update_state().
+ *
  * @param problem_id The ID used to identify the 0D problem.
  * @param ydot The state vector containing all state.ydot degrees-of-freedom.
  */
@@ -358,11 +375,15 @@ void return_ydot(int problem_id, std::vector<double>& ydot) {
   for (int i = 0; i < system_size; i++) {
     ydot[i] = state.ydot[i];
   }
+  commit_persistent_state_snapshot(interface);
 }
 
 /**
  * @brief Update the state vector.
  *
+ * For coupling safety, this restores the last accepted persistent block state
+ * snapshot before applying the new y/ydot vectors.
+ *
  * @param problem_id The ID used to identify the 0D problem.
  * @param new_state_y The new state vector containing all state.y
  * degrees-of-freedom.
@@ -380,6 +401,12 @@ void update_state(int problem_id, std::vector<double> new_state_y,
         "ERROR: State vector size is wrong in update_state().");
   }
 
+  if (interface->committed_block_states_.empty()) {
+    commit_persistent_state_snapshot(interface);
+  } else {
+    interface->model_->set_persistent_states(interface->committed_block_states_);
+  }
+
   auto state = interface->state_;
   for (int i = 0; i < system_size; i++) {
     state.y[i] = new_state_y[i];

src/interface/interface.h

@@ -116,6 +116,10 @@ class SolverInterface {
    * @brief The current 0D state vector
    */
   State state_;
+  /**
+   * @brief Last accepted persistent block states for coupling rollback safety.
+   */
+  std::vector<std::vector<double>> committed_block_states_;
   /**
    * @brief Vector to store solution times
    */

src/model/Block.cpp

@@ -5,7 +5,9 @@
 
 #include "Model.h"
 
-std::string Block::get_name() { return this->model->get_block_name(this->id); }
+std::string Block::get_name() const {
+  return this->model->get_block_name(this->id);
+}
 
 void Block::update_vessel_type(VesselType type) { vessel_type = type; }
 
@@ -61,6 +63,13 @@ void Block::update_solution(
 
 void Block::post_solve(Eigen::Matrix<double, Eigen::Dynamic, 1>& y) {}
 
+void Block::accept_timestep(const Eigen::Matrix<double, Eigen::Dynamic, 1>& y) {
+}
+
+std::vector<double> Block::get_persistent_state() const { return {}; }
+
+void Block::set_persistent_state(const std::vector<double>& state) {}
+
 void Block::update_gradient(Eigen::SparseMatrix<double>& jacobian,
                             Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
                             Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha,

src/model/Block.h

@@ -154,7 +154,7 @@ class Block {
    *
    * @return std::string Name of the block
    */
-  std::string get_name();
+  std::string get_name() const;
 
   /**
    * @brief Update vessel type of the block
@@ -251,6 +251,34 @@ class Block {
    */
   virtual void post_solve(Eigen::Matrix<double, Eigen::Dynamic, 1>& y);
 
+  /**
+   * @brief Accept a converged time-step.
+   *
+   * This hook is called once per accepted integrator step and can be used by
+   * blocks that maintain persistent memory (for example convolution history).
+   *
+   * @param y Accepted state vector
+   */
+  virtual void accept_timestep(
+      const Eigen::Matrix<double, Eigen::Dynamic, 1>& y);
+
+  /**
+   * @brief Get persistent internal state of the block.
+   *
+   * The returned vector must be deterministic and sufficient to restore the
+   * block memory using @ref set_persistent_state.
+   *
+   * @return Persistent state payload
+   */
+  virtual std::vector<double> get_persistent_state() const;
+
+  /**
+   * @brief Restore persistent internal state of the block.
+   *
+   * @param state Persistent state payload from @ref get_persistent_state
+   */
+  virtual void set_persistent_state(const std::vector<double>& state);
+
   /**
    * @brief Set the gradient of the block contributions with respect to the
    * parameters

src/model/BlockType.h

@@ -34,7 +34,8 @@ enum class BlockType {
   linear_elastance_chamber = 18,
   open_loop_coronary_var_res_bc = 19,
   blood_vessel_rc = 20,
-  chamber_elastance_inductor_exponential = 21
+  chamber_elastance_inductor_exponential = 21,
+  impedance_bc = 22
 };
 
 /**

src/model/CMakeLists.txt

@@ -22,6 +22,7 @@ set(CXXSRCS
   ClosedLoopRCRBC.cpp
   DOFHandler.cpp
   FlowReferenceBC.cpp
+  ImpedanceBC.cpp
   Junction.cpp
   Model.cpp
   Node.cpp
@@ -55,6 +56,7 @@ set(HDRS
   ClosedLoopRCRBC.h
   DOFHandler.h
   FlowReferenceBC.h
+  ImpedanceBC.h
   Junction.h
   Model.h
   Node.h
@@ -80,4 +82,3 @@ target_include_directories( ${lib} PUBLIC
 
 target_link_libraries( ${lib} Eigen3::Eigen )
 target_link_libraries( ${lib} nlohmann_json::nlohmann_json )
-