diff --git a/.gitignore b/.gitignore index cd14e84b2..dba1dd985 100644 --- a/.gitignore +++ b/.gitignore @@ -29,6 +29,7 @@ Debug*/ externals/ *.so build*/ +.npm-cache/ # IDE .cproject diff --git a/README.md b/README.md index e0866c3fe..5eb0ff317 100755 --- a/README.md +++ b/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. diff --git a/src/algebra/Integrator.cpp b/src/algebra/Integrator.cpp index ab01c100d..f130262bf 100644 --- a/src/algebra/Integrator.cpp +++ b/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; } diff --git a/src/interface/interface.cpp b/src/interface/interface.cpp index bfd25d977..ab917039e 100644 --- a/src/interface/interface.cpp +++ b/src/interface/interface.cpp @@ -11,6 +11,12 @@ int SolverInterface::problem_id_count_ = 0; std::map 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(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& 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& 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 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]; diff --git a/src/interface/interface.h b/src/interface/interface.h index 14674cc86..1ce4fdcf4 100644 --- a/src/interface/interface.h +++ b/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> committed_block_states_; /** * @brief Vector to store solution times */ diff --git a/src/model/Block.cpp b/src/model/Block.cpp index 7800fa10d..3ee63cdbc 100644 --- a/src/model/Block.cpp +++ b/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& y) {} +void Block::accept_timestep(const Eigen::Matrix& y) { +} + +std::vector Block::get_persistent_state() const { return {}; } + +void Block::set_persistent_state(const std::vector& state) {} + void Block::update_gradient(Eigen::SparseMatrix& jacobian, Eigen::Matrix& residual, Eigen::Matrix& alpha, diff --git a/src/model/Block.h b/src/model/Block.h index 752a5c0c0..916b183d7 100644 --- a/src/model/Block.h +++ b/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& 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& 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 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& state); + /** * @brief Set the gradient of the block contributions with respect to the * parameters diff --git a/src/model/BlockType.h b/src/model/BlockType.h index 7767f0372..245de517b 100644 --- a/src/model/BlockType.h +++ b/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 }; /** diff --git a/src/model/CMakeLists.txt b/src/model/CMakeLists.txt index 88ba98d5f..4bc8564b7 100644 --- a/src/model/CMakeLists.txt +++ b/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 ) - diff --git a/src/model/ImpedanceBC.cpp b/src/model/ImpedanceBC.cpp new file mode 100644 index 000000000..9b60cc4a0 --- /dev/null +++ b/src/model/ImpedanceBC.cpp @@ -0,0 +1,304 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause +#include "ImpedanceBC.h" + +#include +#include +#include +#include +#include + +#include "Model.h" + +void ImpedanceBC::setup_dofs(DOFHandler& dofhandler) { + Block::setup_dofs_(dofhandler, 1, {}); +} + +void ImpedanceBC::update_constant(SparseSystem& system, + std::vector& parameters) { + system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = 1.0; +} + +void ImpedanceBC::update_time(SparseSystem& system, + std::vector& parameters) { + ensure_initialized(); + + double conv_sum = 0.0; + // Match Olufsen startup semantics: do not evaluate lagged terms until one + // full cycle of accepted data has elapsed. + if (num_accepted_steps >= num_period_steps) { + for (int m = 1; m < num_kernel_terms; m++) { + conv_sum += kernel[m] * lagged_flow(m); + } + } + + // Kernel entries are discrete-time convolution coefficients (Olufsen + // periodic kernel), so no extra dt scaling is applied in the algebraic BC. + system.F.coeffRef(global_eqn_ids[0], global_var_ids[1]) = -kernel[0]; + system.C(global_eqn_ids[0]) = -(pd + conv_sum); +} + +void ImpedanceBC::accept_timestep( + const Eigen::Matrix& y) { + ensure_initialized(); + + const double flow_accepted = y[global_var_ids[1]]; + head = (head + 1) % num_period_steps; + flow_history[head] = flow_accepted; + committed_samples = std::min(committed_samples + 1, num_period_steps); + num_accepted_steps += 1; +} + +std::vector ImpedanceBC::get_persistent_state() const { + if (!initialized) { + return {}; + } + + std::vector state; + state.reserve(6 + flow_history.size()); + state.push_back(dt); + state.push_back(static_cast(num_period_steps)); + state.push_back(static_cast(num_kernel_terms)); + state.push_back(static_cast(head)); + state.push_back(static_cast(committed_samples)); + state.push_back(static_cast(num_accepted_steps)); + state.insert(state.end(), flow_history.begin(), flow_history.end()); + return state; +} + +void ImpedanceBC::set_persistent_state(const std::vector& state) { + if (state.empty()) { + initialized = false; + dt = -1.0; + num_period_steps = -1; + num_kernel_terms = -1; + head = -1; + committed_samples = 0; + num_accepted_steps = 0; + flow_history.clear(); + return; + } + + if (state.size() < 5) { + throw std::runtime_error("Invalid persistent state for IMPEDANCE block " + + get_name() + "."); + } + + const double new_dt = state[0]; + const int new_num_period_steps = static_cast(std::llround(state[1])); + const int new_num_kernel_terms = static_cast(std::llround(state[2])); + const int new_head = static_cast(std::llround(state[3])); + const int new_committed_samples = static_cast(std::llround(state[4])); + const bool has_num_accepted_steps = + (state.size() == static_cast(6 + new_num_period_steps)); + const long long new_num_accepted_steps = + has_num_accepted_steps + ? static_cast(std::llround(state[5])) + : static_cast(new_committed_samples); + + if ((new_dt <= 0.0) || (new_num_period_steps <= 0) || + (new_num_kernel_terms <= 0) || + (new_num_kernel_terms > new_num_period_steps) || (new_head < 0) || + (new_head >= new_num_period_steps) || (new_committed_samples < 0) || + (new_committed_samples > new_num_period_steps) || + (new_num_accepted_steps < 0)) { + throw std::runtime_error("Corrupted persistent state for IMPEDANCE block " + + get_name() + "."); + } + + const size_t expected_size = has_num_accepted_steps + ? static_cast(6 + new_num_period_steps) + : static_cast(5 + new_num_period_steps); + if (state.size() != expected_size) { + throw std::runtime_error("Persistent state size mismatch for IMPEDANCE " + "block " + + get_name() + "."); + } + + if (configured && (kernel.size() != static_cast(new_num_period_steps))) { + throw std::runtime_error( + "Persistent IMPEDANCE state is incompatible with configured kernel " + "size in block " + + get_name() + "."); + } + + dt = new_dt; + num_period_steps = new_num_period_steps; + num_kernel_terms = new_num_kernel_terms; + head = new_head; + committed_samples = new_committed_samples; + num_accepted_steps = new_num_accepted_steps; + flow_history.assign(state.begin() + (has_num_accepted_steps ? 6 : 5), + state.end()); + initialized = true; +} + +void ImpedanceBC::configure(const std::vector& z, double pd, + const std::string& convolution_mode, + int num_kernel_terms) { + if (z.empty()) { + throw std::runtime_error("IMPEDANCE block " + get_name() + + " requires non-empty kernel `z`."); + } + for (size_t i = 0; i < z.size(); i++) { + if (!std::isfinite(z[i])) { + throw std::runtime_error("IMPEDANCE block " + get_name() + + " has non-finite kernel value at index " + + std::to_string(i) + "."); + } + } + + std::string mode_lower = convolution_mode; + std::transform(mode_lower.begin(), mode_lower.end(), mode_lower.begin(), + [](unsigned char c) { return std::tolower(c); }); + if (mode_lower.empty()) { + mode_lower = "exact"; + } + + if (mode_lower == "exact") { + mode = ConvolutionMode::exact; + } else if (mode_lower == "truncated") { + mode = ConvolutionMode::truncated; + } else { + throw std::runtime_error("IMPEDANCE block " + get_name() + + " has invalid `convolution_mode` value `" + + convolution_mode + "`. Supported values are " + "`exact` and `truncated`."); + } + + if ((mode == ConvolutionMode::truncated) && (num_kernel_terms <= 0)) { + throw std::runtime_error("IMPEDANCE block " + get_name() + + " in truncated mode requires " + "`num_kernel_terms > 0`."); + } + + kernel = z; + this->pd = pd; + num_kernel_terms_input = num_kernel_terms; + configured = true; + + // Reset runtime state; it will be re-initialized using current dt. + initialized = false; + dt = -1.0; + num_period_steps = -1; + this->num_kernel_terms = -1; + head = -1; + committed_samples = 0; + num_accepted_steps = 0; + flow_history.clear(); +} + +void ImpedanceBC::ensure_initialized() { + if (!configured) { + throw std::runtime_error("IMPEDANCE block " + get_name() + + " was not configured. Missing or invalid " + "impedance kernel input."); + } + + const double model_dt = model->get_time_step_size(); + if (model_dt <= 0.0) { + throw std::runtime_error("IMPEDANCE block " + get_name() + + " requires a positive solver time-step size."); + } + const double period = model->cardiac_cycle_period; + if (period <= 0.0) { + throw std::runtime_error("IMPEDANCE block " + get_name() + + " requires `simulation_parameters.cardiac_period " + "> 0`."); + } + + if (!initialized) { + const double period_over_dt = period / model_dt; + const long rounded = std::lround(period_over_dt); + if (rounded <= 0) { + throw std::runtime_error("IMPEDANCE block " + get_name() + + " has invalid period-to-dt ratio."); + } + + const double mismatch = + std::abs(period_over_dt - static_cast(rounded)); + if (mismatch > 1.0e-8) { + throw std::runtime_error("IMPEDANCE block " + get_name() + + " requires period/dt to be an integer. Got " + "period=" + + std::to_string(period) + + ", dt=" + std::to_string(model_dt) + "."); + } + + const int dt_period_steps = static_cast(rounded); + if (model->impedance_points_per_cycle > 0) { + const int configured_period_steps = model->impedance_points_per_cycle - 1; + if (configured_period_steps <= 0) { + throw std::runtime_error( + "IMPEDANCE block " + get_name() + + " requires `simulation_parameters.number_of_time_pts_per_cardiac_cycle > 1`."); + } + if (configured_period_steps != dt_period_steps) { + throw std::runtime_error( + "IMPEDANCE block " + get_name() + + " requires `simulation_parameters.number_of_time_pts_per_cardiac_cycle - 1` " + "to match cardiac_period / dt. Expected " + + std::to_string(dt_period_steps + 1) + ", got " + + std::to_string(model->impedance_points_per_cycle) + "."); + } + num_period_steps = configured_period_steps; + } else { + num_period_steps = dt_period_steps; + } + if (kernel.size() != static_cast(num_period_steps)) { + throw std::runtime_error( + "IMPEDANCE block " + get_name() + + " requires kernel length `z.size()` to match one period in time " + "steps. Expected " + + std::to_string(num_period_steps) + ", got " + + std::to_string(kernel.size()) + "."); + } + + dt = model_dt; + if (mode == ConvolutionMode::exact) { + num_kernel_terms = num_period_steps; + } else { + num_kernel_terms = num_kernel_terms_input; + } + if ((num_kernel_terms <= 0) || (num_kernel_terms > num_period_steps)) { + throw std::runtime_error( + "IMPEDANCE block " + get_name() + + " has invalid `num_kernel_terms`. It must be in [1, " + + std::to_string(num_period_steps) + "]."); + } + + flow_history.assign(num_period_steps, 0.0); + head = num_period_steps - 1; + committed_samples = 0; + num_accepted_steps = 0; + initialized = true; + return; + } + + const double denom = + std::max({1.0, std::abs(model_dt), std::abs(dt)}); + const double rel_change = std::abs(model_dt - dt) / denom; + if (rel_change > 1.0e-12) { + throw std::runtime_error( + "IMPEDANCE block " + get_name() + + " does not support changing dt after initialization. " + "Reinitialize the model when external step size changes."); + } +} + +double ImpedanceBC::lagged_flow(int lag) const { + if (lag <= 0) { + throw std::runtime_error("Internal error: lag must be positive in " + "IMPEDANCE block " + + get_name() + "."); + } + if (lag > committed_samples) { + return 0.0; + } + + int idx = (head - (lag - 1)) % num_period_steps; + if (idx < 0) { + idx += num_period_steps; + } + return flow_history[idx]; +} diff --git a/src/model/ImpedanceBC.h b/src/model/ImpedanceBC.h new file mode 100644 index 000000000..a7a2ae239 --- /dev/null +++ b/src/model/ImpedanceBC.h @@ -0,0 +1,151 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause +/** + * @file ImpedanceBC.h + * @brief model::ImpedanceBC source file + */ +#ifndef SVZERODSOLVER_MODEL_IMPEDANCEBC_HPP_ +#define SVZERODSOLVER_MODEL_IMPEDANCEBC_HPP_ + +#include +#include + +#include "Block.h" +#include "SparseSystem.h" + +/** + * @brief Structured-tree impedance boundary condition with periodic memory. + * + * This boundary condition applies an Olufsen-style discrete convolution using + * one-cycle flow history: + * + * \f[ + * P_{n+1}=P_d+z_0Q_{n+1}+\sum_{m=1}^{N_k-1}z_mQ_{n+1-m} + * \f] + * + * The `z` entries are discrete-time convolution weights. The `m=0` term is + * treated implicitly through the linear system matrix. + * History terms (`m>=1`) are evaluated from committed (accepted) timesteps only + * via a fixed-size ring buffer of one cardiac period. + * + * ### Usage in json configuration file + * + * "boundary_conditions": [ + * { + * "bc_name": "OUT", + * "bc_type": "IMPEDANCE", + * "bc_values": { + * "Pd": 0.0, + * "z": [ ... ], + * "convolution_mode": "exact", + * "num_kernel_terms": 128 + * } + * } + * ] + * + * `convolution_mode` may be `exact` (default) or `truncated`. + * In `truncated` mode, only the first `num_kernel_terms` kernel entries are + * used. + */ +class ImpedanceBC : public Block { + public: + /** + * @brief Construct a new ImpedanceBC object + * + * @param id Global ID of the block + * @param model The model to which the block belongs + */ + ImpedanceBC(int id, Model* model) + : Block(id, model, BlockType::impedance_bc, BlockClass::boundary_condition, + {{"z", InputParameter(false, false, false)}, + {"Pd", InputParameter(true, false, false)}, + {"convolution_mode", InputParameter(true, false, false)}, + {"num_kernel_terms", InputParameter(true, false, false)}}) {} + + /** + * @brief Set up the algebraic degree of freedom for the impedance BC. + * + * @param dofhandler Global degree-of-freedom handler + */ + void setup_dofs(DOFHandler& dofhandler) override; + + /** + * @brief Assemble constant Jacobian contributions. + * + * @param system Global sparse system + * @param parameters Current parameter values + */ + void update_constant(SparseSystem& system, + std::vector& parameters) override; + + /** + * @brief Assemble time-dependent residual contributions. + * + * @param system Global sparse system + * @param parameters Current parameter values + */ + void update_time(SparseSystem& system, + std::vector& parameters) override; + + /** + * @brief Commit the accepted outlet flow for the next convolution update. + * + * @param y Accepted global state vector + */ + void accept_timestep( + const Eigen::Matrix& y) override; + + /** + * @brief Serialize the convolution history needed for restart/coupling. + * + * @return Persistent block state + */ + std::vector get_persistent_state() const override; + + /** + * @brief Restore a previously serialized convolution history state. + * + * @param state Persistent block state + */ + void set_persistent_state(const std::vector& state) override; + + /** + * @brief Configure the impedance kernel and convolution controls. + * + * @param z Time-domain impedance kernel + * @param pd Distal/reference pressure + * @param convolution_mode Convolution mode (`exact` or `truncated`) + * @param num_kernel_terms Number of kernel terms for truncated mode + */ + void configure(const std::vector& z, double pd, + const std::string& convolution_mode, int num_kernel_terms); + + /** + * @brief Nonzero triplet counts contributed by this block. + */ + TripletsContributions num_triplets{2, 0, 0}; + + private: + enum class ConvolutionMode { exact = 0, truncated = 1 }; + + void ensure_initialized(); + double lagged_flow(int lag) const; + + std::vector kernel; + double pd = 0.0; + ConvolutionMode mode = ConvolutionMode::exact; + int num_kernel_terms_input = -1; + bool configured = false; + + // Runtime-discretized state. + bool initialized = false; + double dt = -1.0; + int num_period_steps = -1; + int num_kernel_terms = -1; + int head = -1; + int committed_samples = 0; + long long num_accepted_steps = 0; + std::vector flow_history; +}; + +#endif // SVZERODSOLVER_MODEL_IMPEDANCEBC_HPP_ diff --git a/src/model/Model.cpp b/src/model/Model.cpp index 5b7fb5066..234231ca2 100644 --- a/src/model/Model.cpp +++ b/src/model/Model.cpp @@ -23,6 +23,7 @@ Model::Model() { {"CORONARY", block_factory()}, {"CORONARY_VAR_RES", block_factory()}, {"FLOW", block_factory()}, + {"IMPEDANCE", block_factory()}, {"NORMAL_JUNCTION", block_factory()}, {"PRESSURE", block_factory()}, {"RCR", block_factory()}, @@ -225,6 +226,47 @@ void Model::post_solve(Eigen::Matrix& y) { } } +void Model::accept_timestep(const Eigen::Matrix& y) { + for (size_t i = 0; i < get_num_blocks(true); i++) { + get_block(i)->accept_timestep(y); + } +} + +std::vector> Model::get_persistent_states() const { + std::vector> states; + states.reserve(get_num_blocks(true)); + + for (size_t i = 0; i < get_num_blocks(true); i++) { + states.push_back(get_block(i)->get_persistent_state()); + } + + return states; +} + +void Model::set_persistent_states( + const std::vector>& states) { + if (states.size() != static_cast(get_num_blocks(true))) { + throw std::runtime_error( + "Persistent-state block count mismatch in Model::set_persistent_states"); + } + + for (size_t i = 0; i < states.size(); i++) { + get_block(i)->set_persistent_state(states[i]); + } +} + +void Model::clear_persistent_states() { + for (size_t i = 0; i < get_num_blocks(true); i++) { + get_block(i)->set_persistent_state({}); + } +} + +void Model::set_time_step_size(double time_step_size) { + this->time_step_size = time_step_size; +} + +double Model::get_time_step_size() const { return time_step_size; } + void Model::to_steady() { for (auto& param : parameters) { param.to_steady(); diff --git a/src/model/Model.h b/src/model/Model.h index 9412e932b..9ea15a283 100644 --- a/src/model/Model.h +++ b/src/model/Model.h @@ -30,6 +30,7 @@ #include "ClosedLoopRCRBC.h" #include "DOFHandler.h" #include "FlowReferenceBC.h" +#include "ImpedanceBC.h" #include "Junction.h" #include "LinearElastanceChamber.h" #include "Node.h" @@ -73,6 +74,7 @@ class Model { double cardiac_cycle_period = -1.0; ///< Cardiac cycle period double time = 0.0; ///< Current time + int impedance_points_per_cycle = 0; ///< One-cycle sample count for IMPEDANCE BCs /** * @brief Create a new block @@ -254,6 +256,46 @@ class Model { */ void post_solve(Eigen::Matrix& y); + /** + * @brief Accept one converged time-step for all blocks. + * + * @param y Accepted state vector + */ + void accept_timestep(const Eigen::Matrix& y); + + /** + * @brief Get persistent state payloads for all blocks (including internal). + * + * @return Per-block persistent states + */ + std::vector> get_persistent_states() const; + + /** + * @brief Restore persistent states for all blocks. + * + * @param states Per-block persistent states + */ + void set_persistent_states(const std::vector>& states); + + /** + * @brief Clear persistent state for all blocks. + */ + void clear_persistent_states(); + + /** + * @brief Set model time step size used by blocks with dt-dependent logic. + * + * @param time_step_size Model time step size + */ + void set_time_step_size(double time_step_size); + + /** + * @brief Get model time step size used by blocks with dt-dependent logic. + * + * @return double Model time step size + */ + double get_time_step_size() const; + /** * @brief Convert the blocks to a steady behavior * @@ -353,6 +395,7 @@ class Model { bool has_windkessel_bc = false; double largest_windkessel_time_constant = 0.0; + double time_step_size = -1.0; }; #endif // SVZERODSOLVER_MODEL_MODEL_HPP_ diff --git a/src/solve/SimulationParameters.cpp b/src/solve/SimulationParameters.cpp index 036609994..b1870206b 100644 --- a/src/solve/SimulationParameters.cpp +++ b/src/solve/SimulationParameters.cpp @@ -2,6 +2,10 @@ // University of California, and others. SPDX-License-Identifier: BSD-3-Clause #include "SimulationParameters.h" +#include + +#include "ImpedanceBC.h" + bool get_param_scalar(const nlohmann::json& data, const std::string& name, const InputParameter& param, double& val) { if (data.contains(name)) { @@ -207,16 +211,31 @@ void validate_input(const nlohmann::json& config) { } } +bool config_has_impedance_bc(const nlohmann::json& config) { + if (!config.contains("boundary_conditions") || + !config["boundary_conditions"].is_array()) { + return false; + } + for (const auto& bc_config : config["boundary_conditions"]) { + if (bc_config.value("bc_type", "") == "IMPEDANCE") { + return true; + } + } + return false; +} + SimulationParameters load_simulation_params(const nlohmann::json& config) { DEBUG_MSG("Loading simulation parameters"); SimulationParameters sim_params; const auto& sim_config = config["simulation_parameters"]; sim_params.sim_coupled = sim_config.value("coupled_simulation", false); + const bool has_impedance = config_has_impedance_bc(config); if (!sim_params.sim_coupled) { sim_params.sim_num_cycles = sim_config["number_of_cardiac_cycles"]; sim_params.sim_pts_per_cycle = sim_config["number_of_time_pts_per_cardiac_cycle"]; + sim_params.sim_impedance_pts_per_cycle = sim_params.sim_pts_per_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 = @@ -232,8 +251,23 @@ SimulationParameters load_simulation_params(const nlohmann::json& config) { sim_params.sim_num_cycles = 1; sim_params.sim_num_time_steps = sim_config["number_of_time_pts"]; sim_params.sim_pts_per_cycle = sim_params.sim_num_time_steps; + sim_params.sim_impedance_pts_per_cycle = + sim_config.value("number_of_time_pts_per_cardiac_cycle", 0); sim_params.sim_external_step_size = sim_config.value("external_step_size", 0.1); + if (has_impedance) { + if (sim_params.sim_num_time_steps != 2) { + throw std::runtime_error( + "Coupled configs with IMPEDANCE boundary conditions require " + "`simulation_parameters.number_of_time_pts = 2`."); + } + if (sim_config.contains("number_of_time_pts_per_cardiac_cycle") && + sim_params.sim_impedance_pts_per_cycle <= 1) { + throw std::runtime_error( + "Coupled configs with IMPEDANCE boundary conditions require " + "`simulation_parameters.number_of_time_pts_per_cardiac_cycle > 1`."); + } + } } sim_params.sim_abs_tol = sim_config.value("absolute_tolerance", 1e-8); sim_params.sim_nliter = sim_config.value("maximum_nonlinear_iterations", 30); @@ -383,6 +417,31 @@ void create_boundary_conditions(Model& model, const nlohmann::json& config, model.get_largest_windkessel_time_constant(), time_constant)); } + if (block->block_type == BlockType::impedance_bc) { + auto* impedance_block = dynamic_cast(block); + if (impedance_block == nullptr) { + throw std::runtime_error( + "Internal error: IMPEDANCE block creation failed for " + bc_name + + "."); + } + + if (!bc_values.contains("z")) { + throw std::runtime_error("IMPEDANCE block " + bc_name + + " is missing required `z` kernel."); + } + if (!bc_values["z"].is_array()) { + throw std::runtime_error("IMPEDANCE block " + bc_name + + " requires `z` to be an array."); + } + std::vector z = bc_values["z"].get>(); + + double pd = bc_values.value("Pd", 0.0); + std::string convolution_mode = bc_values.value("convolution_mode", "exact"); + int num_kernel_terms = bc_values.value("num_kernel_terms", -1); + + impedance_block->configure(z, pd, convolution_mode, num_kernel_terms); + } + if (block->block_type == BlockType::closed_loop_rcr_bc) { if (bc_values["closed_loop_outlet"] == true) { closed_loop_bcs.push_back(bc_name); @@ -446,6 +505,7 @@ void create_external_coupling( if (coupling_loc == "inlet") { std::vector possible_types = {"RESISTANCE", "RCR", + "IMPEDANCE", "ClosedLoopRCR", "SimplifiedRCR", "CORONARY", diff --git a/src/solve/SimulationParameters.h b/src/solve/SimulationParameters.h index b84b2845a..755505104 100644 --- a/src/solve/SimulationParameters.h +++ b/src/solve/SimulationParameters.h @@ -43,6 +43,8 @@ struct SimulationParameters { double sim_cardiac_period{-1.0}; ///< Cardiac period 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_impedance_pts_per_cycle{ + 0}; ///< Number of sample points spanning one impedance 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 diff --git a/src/solve/Solver.cpp b/src/solve/Solver.cpp index 7f9bf06a7..3ae89e4f9 100644 --- a/src/solve/Solver.cpp +++ b/src/solve/Solver.cpp @@ -11,6 +11,7 @@ Solver::Solver(const nlohmann::json& config) { simparams = load_simulation_params(config); DEBUG_MSG("Load model"); this->model = std::shared_ptr(new Model()); + this->model->impedance_points_per_cycle = simparams.sim_impedance_pts_per_cycle; load_simulation_model(config, *this->model.get()); // If period isn't specified anywhere, set to 1 @@ -74,6 +75,7 @@ void Solver::setup_initial() { } this->model->to_unsteady(); + this->model->clear_persistent_states(); } // Use the initial condition (steady or user-provided) to set up parameters diff --git a/tests/cases/dirgraph-results/pulsatileFlow_R_impedance_directed_graph.dot b/tests/cases/dirgraph-results/pulsatileFlow_R_impedance_directed_graph.dot new file mode 100644 index 000000000..93378161d --- /dev/null +++ b/tests/cases/dirgraph-results/pulsatileFlow_R_impedance_directed_graph.dot @@ -0,0 +1,7 @@ +strict digraph { +BC0_inlet; +V0; +BC0_outlet; +BC0_inlet -> V0; +V0 -> BC0_outlet; +} diff --git a/tests/cases/pulsatileFlow_R_impedance.json b/tests/cases/pulsatileFlow_R_impedance.json new file mode 100644 index 000000000..9e2a80d5c --- /dev/null +++ b/tests/cases/pulsatileFlow_R_impedance.json @@ -0,0 +1,152 @@ +{ + "description": { + "description of test case": "pulsatile flow -> R -> IMPEDANCE (periodic convolution)" + }, + "boundary_conditions": [ + { + "bc_name": "INFLOW", + "bc_type": "FLOW", + "bc_values": { + "Q": [ + 2.2, 2.273601511491127, 2.347025710494341, 2.4200957116830426, + 2.4926354830241926, 2.564470269854896, 2.6354270158816937, + 2.7053347800883305, 2.7740251485476346, 2.8413326401454233, + 2.9070951052389966, 2.9711541162898327, 3.0333553495294034, + 3.0935489567386503, 3.1515899262454683, 3.207338432270528, + 3.260660171779821, 3.311426688032439, 3.3595156800441055, + 3.4048112972209674, 3.447204418453818, 3.4865929150004082, + 3.5228818965225326, 3.555983939685165, 3.58581929876693, + 3.6123160977745314, 3.6354105035983135, 3.655046879791816, + 3.671177920604846, 3.6837647649471714, 3.6927770900082955, + 3.6981931843077587, 3.7, 3.6981931843077587, 3.6927770900082955, + 3.6837647649471714, 3.671177920604846, 3.655046879791816, + 3.6354105035983135, 3.6123160977745314, 3.58581929876693, + 3.5559839396851656, 3.5228818965225326, 3.4865929150004082, + 3.4472044184538184, 3.4048112972209674, 3.359515680044106, + 3.311426688032439, 3.2606601717798216, 3.207338432270528, + 3.1515899262454683, 3.0935489567386503, 3.0333553495294034, + 2.9711541162898327, 2.907095105238997, 2.8413326401454233, + 2.774025148547635, 2.7053347800883305, 2.6354270158816937, + 2.5644702698548962, 2.492635483024193, 2.420095711683043, + 2.3470257104943415, 2.273601511491127, 2.2, 2.1263984885088734, + 2.0529742895056593, 1.9799042883169578, 1.9073645169758078, + 1.8355297301451046, 1.764572984118307, 1.69466521991167, + 1.6259748514523658, 1.5586673598545775, 1.4929048947610037, + 1.4288458837101679, 1.3666446504705974, 1.3064510432613503, + 1.2484100737545323, 1.1926615677294725, 1.139339828220179, + 1.0885733119675618, 1.0404843199558953, 0.9951887027790325, + 0.9527955815461824, 0.9134070849995921, 0.8771181034774678, + 0.8440160603148354, 0.8141807012330704, 0.7876839022254689, + 0.7645894964016868, 0.7449531202081843, 0.7288220793951548, + 0.7162352350528289, 0.7072229099917049, 0.7018068156922417, + 0.7000000000000002, 0.7018068156922417, 0.7072229099917049, + 0.7162352350528289, 0.7288220793951545, 0.7449531202081843, + 0.7645894964016868, 0.7876839022254689, 0.8141807012330702, + 0.8440160603148352, 0.8771181034774675, 0.9134070849995919, + 0.9527955815461819, 0.9951887027790323, 1.0404843199558949, + 1.0885733119675616, 1.1393398282201788, 1.192661567729472, + 1.2484100737545312, 1.3064510432613503, 1.366644650470597, + 1.4288458837101674, 1.4929048947610033, 1.5586673598545764, + 1.6259748514523646, 1.6946652199116703, 1.7645729841183064, + 1.835529730145104, 1.9073645169758071, 1.9799042883169566, + 2.0529742895056593, 2.126398488508873, 2.1999999999999997 + ], + "t": [ + 0.0, 0.0078125, 0.015625, 0.0234375, 0.03125, 0.0390625, 0.046875, + 0.0546875, 0.0625, 0.0703125, 0.078125, 0.0859375, 0.09375, 0.1015625, + 0.109375, 0.1171875, 0.125, 0.1328125, 0.140625, 0.1484375, 0.15625, + 0.1640625, 0.171875, 0.1796875, 0.1875, 0.1953125, 0.203125, + 0.2109375, 0.21875, 0.2265625, 0.234375, 0.2421875, 0.25, 0.2578125, + 0.265625, 0.2734375, 0.28125, 0.2890625, 0.296875, 0.3046875, 0.3125, + 0.3203125, 0.328125, 0.3359375, 0.34375, 0.3515625, 0.359375, + 0.3671875, 0.375, 0.3828125, 0.390625, 0.3984375, 0.40625, 0.4140625, + 0.421875, 0.4296875, 0.4375, 0.4453125, 0.453125, 0.4609375, 0.46875, + 0.4765625, 0.484375, 0.4921875, 0.5, 0.5078125, 0.515625, 0.5234375, + 0.53125, 0.5390625, 0.546875, 0.5546875, 0.5625, 0.5703125, 0.578125, + 0.5859375, 0.59375, 0.6015625, 0.609375, 0.6171875, 0.625, 0.6328125, + 0.640625, 0.6484375, 0.65625, 0.6640625, 0.671875, 0.6796875, 0.6875, + 0.6953125, 0.703125, 0.7109375, 0.71875, 0.7265625, 0.734375, + 0.7421875, 0.75, 0.7578125, 0.765625, 0.7734375, 0.78125, 0.7890625, + 0.796875, 0.8046875, 0.8125, 0.8203125, 0.828125, 0.8359375, 0.84375, + 0.8515625, 0.859375, 0.8671875, 0.875, 0.8828125, 0.890625, 0.8984375, + 0.90625, 0.9140625, 0.921875, 0.9296875, 0.9375, 0.9453125, 0.953125, + 0.9609375, 0.96875, 0.9765625, 0.984375, 0.9921875, 1.0 + ] + } + }, + { + "bc_name": "OUT", + "bc_type": "IMPEDANCE", + "bc_values": { + "Pd": 80000.0, + "z": [ + 915.0, 752.7, 619.5659999999998, 510.3490799999999, 420.7451063999999, + 347.22467083199984, 286.8936399945599, 237.37880650961267, + 196.73412261767453, 163.3638118006894, 135.95932030567766, + 113.44761738718579, 94.94880149929185, 79.74133736638282, + 67.23355037465817, 56.94025196675946, 48.46357245909172, + 41.4772439458772, 35.713712274452845, 30.953568859108216, + 27.016884762644473, 23.756104637580695, 21.05021975238415, + 18.79998986753163, 16.924025168541025, 15.355573445825426, + 14.0398855770462, 12.932055217617847, 11.9952473419978, + 11.199245642718338, 10.519261392863577, 9.93495670866598, + 9.429643620293596, 8.98962730544449, 8.60366753653215, + 8.262537063102675, 7.958659481227581, 7.685812282300335, + 7.438883349026119, 7.213671278190467, 7.006721641465449, + 6.81519271548395, 6.636745376789475, 6.469452812058153, + 6.311726476916653, 6.162255378680043, 6.019956284773855, + 5.883932890285658, 5.753442332069913, 5.62786772709229, + 5.506695650714739, 5.389497665795166, 5.2759151735169345, + 5.165646988097487, 5.0584391451372674, 4.954076541611942, + 4.852376077869186, 4.7531810313251786, 4.656356440209644, + 4.561785315604446, 4.469365532735532, 4.379007279304225, + 4.290630960641332, 4.204165479505522, 4.119546823138565, + 4.03671690231878, 3.955622597099653, 3.8762149720760033, + 3.798448630707525, 3.7222811837132688, 3.6476728110473164, + 3.5745859006531875, 3.502984750218114, 3.432835320627703, + 3.3641050318545895, 3.296762593681838, 3.230777865028961, + 3.166121736769405, 3.102766033847656, 3.0406834332578865, + 2.9798473950642195, 2.9202321041495574, 2.861812420795597, + 2.8045638385374896, 2.7484624480161397, 2.6934849057803247, + 2.6396084071788155, 2.5868106626367977, 2.5350698767373405, + 2.484364729632282, 2.434674360391981, 2.3859783519730637, + 2.338256717540519, 2.29148988792738, 2.2456587000537227, + 2.2007443861582585, 2.1567285637216935, 2.1135932259822723, + 2.071320732961337, 2.0298938029310527, 1.989295504268164, + 1.9495092476473015, 1.9105187785352462, 1.8723081699540718, + 1.8348618154864051, 1.7981644225004372, 1.7622010055759116, + 1.7269568801152897, 1.692417656126719, 1.6585692321674472, + 1.6253977894379736, 1.592889786018592, 1.56103195124111, + 1.5298112801894574, 1.4992150283236676, 1.4692307062223533, + 1.4398460744393367, 1.411049138470523, 1.3828281438274905, + 1.3551715712145709, 1.3280681318064558, 1.3015067626235914, + 1.2754766220027967, 1.2499670851607159, 1.2249677398478411, + 1.2004683820909632, 1.1764590120220082, 1.1529298297913166 + ] + } + } + ], + "simulation_parameters": { + "cardiac_period": 1.0, + "number_of_cardiac_cycles": 6, + "number_of_time_pts_per_cardiac_cycle": 129, + "output_all_cycles": false, + "steady_initial": false, + "absolute_tolerance": 1e-9 + }, + "vessels": [ + { + "boundary_conditions": { + "inlet": "INFLOW", + "outlet": "OUT" + }, + "vessel_id": 0, + "vessel_length": 10.0, + "vessel_name": "branch0_seg0", + "zero_d_element_type": "BloodVessel", + "zero_d_element_values": { + "R_poiseuille": 75.0 + } + } + ] +} diff --git a/tests/cases/results/result_pulsatileFlow_R_impedance.json b/tests/cases/results/result_pulsatileFlow_R_impedance.json new file mode 100644 index 000000000..1423e62b6 --- /dev/null +++ b/tests/cases/results/result_pulsatileFlow_R_impedance.json @@ -0,0 +1,788 @@ +{ + "name": { + "0": "branch0_seg0", + "1": "branch0_seg0", + "2": "branch0_seg0", + "3": "branch0_seg0", + "4": "branch0_seg0", + "5": "branch0_seg0", + "6": "branch0_seg0", + "7": "branch0_seg0", + "8": "branch0_seg0", + "9": "branch0_seg0", + "10": "branch0_seg0", + "11": "branch0_seg0", + "12": "branch0_seg0", + "13": "branch0_seg0", + "14": "branch0_seg0", + "15": "branch0_seg0", + "16": "branch0_seg0", + "17": "branch0_seg0", + "18": "branch0_seg0", + "19": "branch0_seg0", + "20": "branch0_seg0", + "21": "branch0_seg0", + "22": "branch0_seg0", + "23": "branch0_seg0", + "24": "branch0_seg0", + "25": "branch0_seg0", + "26": "branch0_seg0", + "27": "branch0_seg0", + "28": "branch0_seg0", + "29": "branch0_seg0", + "30": "branch0_seg0", + "31": "branch0_seg0", + "32": "branch0_seg0", + "33": "branch0_seg0", + "34": "branch0_seg0", + "35": "branch0_seg0", + "36": "branch0_seg0", + "37": "branch0_seg0", + "38": "branch0_seg0", + "39": "branch0_seg0", + "40": "branch0_seg0", + "41": "branch0_seg0", + "42": "branch0_seg0", + "43": "branch0_seg0", + "44": "branch0_seg0", + "45": "branch0_seg0", + "46": "branch0_seg0", + "47": "branch0_seg0", + "48": "branch0_seg0", + "49": "branch0_seg0", + "50": "branch0_seg0", + "51": "branch0_seg0", + "52": "branch0_seg0", + "53": "branch0_seg0", + "54": "branch0_seg0", + "55": "branch0_seg0", + "56": "branch0_seg0", + "57": "branch0_seg0", + "58": "branch0_seg0", + "59": "branch0_seg0", + "60": "branch0_seg0", + "61": "branch0_seg0", + "62": "branch0_seg0", + "63": "branch0_seg0", + "64": "branch0_seg0", + "65": "branch0_seg0", + "66": "branch0_seg0", + "67": "branch0_seg0", + "68": "branch0_seg0", + "69": "branch0_seg0", + "70": "branch0_seg0", + "71": "branch0_seg0", + "72": "branch0_seg0", + "73": "branch0_seg0", + "74": "branch0_seg0", + "75": "branch0_seg0", + "76": "branch0_seg0", + "77": "branch0_seg0", + "78": "branch0_seg0", + "79": "branch0_seg0", + "80": "branch0_seg0", + "81": "branch0_seg0", + "82": "branch0_seg0", + "83": "branch0_seg0", + "84": "branch0_seg0", + "85": "branch0_seg0", + "86": "branch0_seg0", + "87": "branch0_seg0", + "88": "branch0_seg0", + "89": "branch0_seg0", + "90": "branch0_seg0", + "91": "branch0_seg0", + "92": "branch0_seg0", + "93": "branch0_seg0", + "94": "branch0_seg0", + "95": "branch0_seg0", + "96": "branch0_seg0", + "97": "branch0_seg0", + "98": "branch0_seg0", + "99": "branch0_seg0", + "100": "branch0_seg0", + "101": "branch0_seg0", + "102": "branch0_seg0", + "103": "branch0_seg0", + "104": "branch0_seg0", + "105": "branch0_seg0", + "106": "branch0_seg0", + "107": "branch0_seg0", + "108": "branch0_seg0", + "109": "branch0_seg0", + "110": "branch0_seg0", + "111": "branch0_seg0", + "112": "branch0_seg0", + "113": "branch0_seg0", + "114": "branch0_seg0", + "115": "branch0_seg0", + "116": "branch0_seg0", + "117": "branch0_seg0", + "118": "branch0_seg0", + "119": "branch0_seg0", + "120": "branch0_seg0", + "121": "branch0_seg0", + "122": "branch0_seg0", + "123": "branch0_seg0", + "124": "branch0_seg0", + "125": "branch0_seg0", + "126": "branch0_seg0", + "127": "branch0_seg0", + "128": "branch0_seg0" + }, + "time": { + "0": 0.0, + "1": 0.0078125, + "2": 0.015625, + "3": 0.0234375, + "4": 0.03125, + "5": 0.0390625, + "6": 0.046875, + "7": 0.0546875, + "8": 0.0625, + "9": 0.0703125, + "10": 0.078125, + "11": 0.0859375, + "12": 0.09375, + "13": 0.1015625, + "14": 0.109375, + "15": 0.1171875, + "16": 0.125, + "17": 0.1328125, + "18": 0.140625, + "19": 0.1484375, + "20": 0.15625, + "21": 0.1640625, + "22": 0.171875, + "23": 0.1796875, + "24": 0.1875, + "25": 0.1953125, + "26": 0.203125, + "27": 0.2109375, + "28": 0.21875, + "29": 0.2265625, + "30": 0.234375, + "31": 0.2421875, + "32": 0.25, + "33": 0.2578125, + "34": 0.265625, + "35": 0.2734375, + "36": 0.28125, + "37": 0.2890625, + "38": 0.296875, + "39": 0.3046875, + "40": 0.3125, + "41": 0.3203125, + "42": 0.328125, + "43": 0.3359375, + "44": 0.34375, + "45": 0.3515625, + "46": 0.359375, + "47": 0.3671875, + "48": 0.375, + "49": 0.3828125, + "50": 0.390625, + "51": 0.3984375, + "52": 0.40625, + "53": 0.4140625, + "54": 0.421875, + "55": 0.4296875, + "56": 0.4375, + "57": 0.4453125, + "58": 0.453125, + "59": 0.4609375, + "60": 0.46875, + "61": 0.4765625, + "62": 0.484375, + "63": 0.4921875, + "64": 0.5, + "65": 0.5078125, + "66": 0.515625, + "67": 0.5234375, + "68": 0.53125, + "69": 0.5390625, + "70": 0.546875, + "71": 0.5546875, + "72": 0.5625, + "73": 0.5703125, + "74": 0.578125, + "75": 0.5859375, + "76": 0.59375, + "77": 0.6015625, + "78": 0.609375, + "79": 0.6171875, + "80": 0.625, + "81": 0.6328125, + "82": 0.640625, + "83": 0.6484375, + "84": 0.65625, + "85": 0.6640625, + "86": 0.671875, + "87": 0.6796875, + "88": 0.6875, + "89": 0.6953125, + "90": 0.703125, + "91": 0.7109375, + "92": 0.71875, + "93": 0.7265625, + "94": 0.734375, + "95": 0.7421875, + "96": 0.75, + "97": 0.7578125, + "98": 0.765625, + "99": 0.7734375, + "100": 0.78125, + "101": 0.7890625, + "102": 0.796875, + "103": 0.8046875, + "104": 0.8125, + "105": 0.8203125, + "106": 0.828125, + "107": 0.8359375, + "108": 0.84375, + "109": 0.8515625, + "110": 0.859375, + "111": 0.8671875, + "112": 0.875, + "113": 0.8828125, + "114": 0.890625, + "115": 0.8984375, + "116": 0.90625, + "117": 0.9140625, + "118": 0.921875, + "119": 0.9296875, + "120": 0.9375, + "121": 0.9453125, + "122": 0.953125, + "123": 0.9609375, + "124": 0.96875, + "125": 0.9765625, + "126": 0.984375, + "127": 0.9921875, + "128": 1.0 + }, + "flow_in": { + "0": 2.2, + "1": 2.2736015115, + "2": 2.3470257105, + "3": 2.4200957117, + "4": 2.492635483, + "5": 2.5644702699, + "6": 2.6354270159, + "7": 2.7053347801, + "8": 2.7740251485, + "9": 2.8413326401, + "10": 2.9070951052, + "11": 2.9711541163, + "12": 3.0333553495, + "13": 3.0935489567, + "14": 3.1515899262, + "15": 3.2073384323, + "16": 3.2606601718, + "17": 3.311426688, + "18": 3.35951568, + "19": 3.4048112972, + "20": 3.4472044185, + "21": 3.486592915, + "22": 3.5228818965, + "23": 3.5559839397, + "24": 3.5858192988, + "25": 3.6123160978, + "26": 3.6354105036, + "27": 3.6550468798, + "28": 3.6711779206, + "29": 3.6837647649, + "30": 3.69277709, + "31": 3.6981931843, + "32": 3.7, + "33": 3.6981931843, + "34": 3.69277709, + "35": 3.6837647649, + "36": 3.6711779206, + "37": 3.6550468798, + "38": 3.6354105036, + "39": 3.6123160978, + "40": 3.5858192988, + "41": 3.5559839397, + "42": 3.5228818965, + "43": 3.486592915, + "44": 3.4472044185, + "45": 3.4048112972, + "46": 3.35951568, + "47": 3.311426688, + "48": 3.2606601718, + "49": 3.2073384323, + "50": 3.1515899262, + "51": 3.0935489567, + "52": 3.0333553495, + "53": 2.9711541163, + "54": 2.9070951052, + "55": 2.8413326401, + "56": 2.7740251485, + "57": 2.7053347801, + "58": 2.6354270159, + "59": 2.5644702699, + "60": 2.492635483, + "61": 2.4200957117, + "62": 2.3470257105, + "63": 2.2736015115, + "64": 2.2, + "65": 2.1263984885, + "66": 2.0529742895, + "67": 1.9799042883, + "68": 1.907364517, + "69": 1.8355297301, + "70": 1.7645729841, + "71": 1.6946652199, + "72": 1.6259748515, + "73": 1.5586673599, + "74": 1.4929048948, + "75": 1.4288458837, + "76": 1.3666446505, + "77": 1.3064510433, + "78": 1.2484100738, + "79": 1.1926615677, + "80": 1.1393398282, + "81": 1.088573312, + "82": 1.04048432, + "83": 0.9951887028, + "84": 0.9527955815, + "85": 0.913407085, + "86": 0.8771181035, + "87": 0.8440160603, + "88": 0.8141807012, + "89": 0.7876839022, + "90": 0.7645894964, + "91": 0.7449531202, + "92": 0.7288220794, + "93": 0.7162352351, + "94": 0.70722291, + "95": 0.7018068157, + "96": 0.7, + "97": 0.7018068157, + "98": 0.70722291, + "99": 0.7162352351, + "100": 0.7288220794, + "101": 0.7449531202, + "102": 0.7645894964, + "103": 0.7876839022, + "104": 0.8141807012, + "105": 0.8440160603, + "106": 0.8771181035, + "107": 0.913407085, + "108": 0.9527955815, + "109": 0.9951887028, + "110": 1.04048432, + "111": 1.088573312, + "112": 1.1393398282, + "113": 1.1926615677, + "114": 1.2484100738, + "115": 1.3064510433, + "116": 1.3666446505, + "117": 1.4288458837, + "118": 1.4929048948, + "119": 1.5586673599, + "120": 1.6259748515, + "121": 1.6946652199, + "122": 1.7645729841, + "123": 1.8355297301, + "124": 1.907364517, + "125": 1.9799042883, + "126": 2.0529742895, + "127": 2.1263984885, + "128": 2.2 + }, + "flow_out": { + "0": 2.2, + "1": 2.2736015115, + "2": 2.3470257105, + "3": 2.4200957117, + "4": 2.492635483, + "5": 2.5644702699, + "6": 2.6354270159, + "7": 2.7053347801, + "8": 2.7740251485, + "9": 2.8413326401, + "10": 2.9070951052, + "11": 2.9711541163, + "12": 3.0333553495, + "13": 3.0935489567, + "14": 3.1515899262, + "15": 3.2073384323, + "16": 3.2606601718, + "17": 3.311426688, + "18": 3.35951568, + "19": 3.4048112972, + "20": 3.4472044185, + "21": 3.486592915, + "22": 3.5228818965, + "23": 3.5559839397, + "24": 3.5858192988, + "25": 3.6123160978, + "26": 3.6354105036, + "27": 3.6550468798, + "28": 3.6711779206, + "29": 3.6837647649, + "30": 3.69277709, + "31": 3.6981931843, + "32": 3.7, + "33": 3.6981931843, + "34": 3.69277709, + "35": 3.6837647649, + "36": 3.6711779206, + "37": 3.6550468798, + "38": 3.6354105036, + "39": 3.6123160978, + "40": 3.5858192988, + "41": 3.5559839397, + "42": 3.5228818965, + "43": 3.486592915, + "44": 3.4472044185, + "45": 3.4048112972, + "46": 3.35951568, + "47": 3.311426688, + "48": 3.2606601718, + "49": 3.2073384323, + "50": 3.1515899262, + "51": 3.0935489567, + "52": 3.0333553495, + "53": 2.9711541163, + "54": 2.9070951052, + "55": 2.8413326401, + "56": 2.7740251485, + "57": 2.7053347801, + "58": 2.6354270159, + "59": 2.5644702699, + "60": 2.492635483, + "61": 2.4200957117, + "62": 2.3470257105, + "63": 2.2736015115, + "64": 2.2, + "65": 2.1263984885, + "66": 2.0529742895, + "67": 1.9799042883, + "68": 1.907364517, + "69": 1.8355297301, + "70": 1.7645729841, + "71": 1.6946652199, + "72": 1.6259748515, + "73": 1.5586673599, + "74": 1.4929048948, + "75": 1.4288458837, + "76": 1.3666446505, + "77": 1.3064510433, + "78": 1.2484100738, + "79": 1.1926615677, + "80": 1.1393398282, + "81": 1.088573312, + "82": 1.04048432, + "83": 0.9951887028, + "84": 0.9527955815, + "85": 0.913407085, + "86": 0.8771181035, + "87": 0.8440160603, + "88": 0.8141807012, + "89": 0.7876839022, + "90": 0.7645894964, + "91": 0.7449531202, + "92": 0.7288220794, + "93": 0.7162352351, + "94": 0.70722291, + "95": 0.7018068157, + "96": 0.7, + "97": 0.7018068157, + "98": 0.70722291, + "99": 0.7162352351, + "100": 0.7288220794, + "101": 0.7449531202, + "102": 0.7645894964, + "103": 0.7876839022, + "104": 0.8141807012, + "105": 0.8440160603, + "106": 0.8771181035, + "107": 0.913407085, + "108": 0.9527955815, + "109": 0.9951887028, + "110": 1.04048432, + "111": 1.088573312, + "112": 1.1393398282, + "113": 1.1926615677, + "114": 1.2484100738, + "115": 1.3064510433, + "116": 1.3666446505, + "117": 1.4288458837, + "118": 1.4929048948, + "119": 1.5586673599, + "120": 1.6259748515, + "121": 1.6946652199, + "122": 1.7645729841, + "123": 1.8355297301, + "124": 1.907364517, + "125": 1.9799042883, + "126": 2.0529742895, + "127": 2.1263984885, + "128": 2.2 + }, + "pressure_in": { + "0": 90844.4239895008, + "1": 91210.3607938189, + "2": 91579.8638990727, + "3": 91952.0431399168, + "4": 92326.0019039613, + "5": 92700.8392917884, + "6": 93075.652287299, + "7": 93449.5379331573, + "8": 93821.5955060938, + "9": 94190.928686827, + "10": 94556.6477193745, + "11": 94917.8715545541, + "12": 95273.7299725069, + "13": 95623.365679135, + "14": 95965.9363713963, + "15": 96300.6167664877, + "16": 96626.6005900228, + "17": 96943.1025184179, + "18": 97249.3600708053, + "19": 97544.6354459162, + "20": 97828.2172995089, + "21": 98099.4224580593, + "22": 98357.5975645856, + "23": 98602.1206526428, + "24": 98832.4026446941, + "25": 99047.8887712502, + "26": 99248.0599073579, + "27": 99432.4338232176, + "28": 99600.5663459164, + "29": 99752.0524294807, + "30": 99886.5271306666, + "31": 100003.6664881402, + "32": 100103.1883029292, + "33": 100184.8528182647, + "34": 100248.4632971763, + "35": 100293.8664964487, + "36": 100320.953035798, + "37": 100329.6576613784, + "38": 100319.9594029844, + "39": 100291.88162457, + "40": 100245.4919679627, + "41": 100180.9021899084, + "42": 100098.2678928401, + "43": 99997.7881500171, + "44": 99879.7050259409, + "45": 99744.3029932004, + "46": 99591.9082471521, + "47": 99422.8879200875, + "48": 99237.649196779, + "49": 99036.638333536, + "50": 98820.3395831345, + "51": 98589.2740282098, + "52": 98343.9983259226, + "53": 98085.1033669234, + "54": 97813.2128518449, + "55": 97528.9817887527, + "56": 97233.0949151735, + "57": 96926.2650485025, + "58": 96609.2313687635, + "59": 96282.7576378598, + "60": 95947.630359605, + "61": 95604.6568849662, + "62": 95254.6634670845, + "63": 94898.4932707594, + "64": 94537.0043411909, + "65": 94171.0675368727, + "66": 93801.5644316189, + "67": 93429.3851907748, + "68": 93055.4264267303, + "69": 92680.5890389032, + "70": 92305.7760433926, + "71": 91931.8903975343, + "72": 91559.8328245979, + "73": 91190.4996438647, + "74": 90824.7806113171, + "75": 90463.5567761376, + "76": 90107.6983581847, + "77": 89758.0626515567, + "78": 89415.4919592954, + "79": 89080.811564204, + "80": 88754.8277406689, + "81": 88438.3258122737, + "82": 88132.0682598864, + "83": 87836.7928847755, + "84": 87553.2110311828, + "85": 87282.0058726324, + "86": 87023.8307661061, + "87": 86779.3076780489, + "88": 86549.0256859976, + "89": 86333.5395594415, + "90": 86133.3684233338, + "91": 85948.9945074741, + "92": 85780.8619847753, + "93": 85629.375901211, + "94": 85494.9012000251, + "95": 85377.7618425515, + "96": 85278.2400277625, + "97": 85196.5755124269, + "98": 85132.9650335153, + "99": 85087.5618342429, + "100": 85060.4752948936, + "101": 85051.7706693133, + "102": 85061.4689277073, + "103": 85089.5467061217, + "104": 85135.936362729, + "105": 85200.5261407832, + "106": 85283.1604378516, + "107": 85383.6401806746, + "108": 85501.7233047507, + "109": 85637.1253374913, + "110": 85789.5200835396, + "111": 85958.5404106041, + "112": 86143.7791339126, + "113": 86344.7899971556, + "114": 86561.0887475571, + "115": 86792.1543024818, + "116": 87037.430004769, + "117": 87296.3249637682, + "118": 87568.2154788467, + "119": 87852.4465419389, + "120": 88148.3334155181, + "121": 88455.1632821892, + "122": 88772.1969619282, + "123": 89098.6706928318, + "124": 89433.7979710866, + "125": 89776.7714457255, + "126": 90126.7648636072, + "127": 90482.9350599322, + "128": 90844.4239895008 + }, + "pressure_out": { + "0": 90679.4239895009, + "1": 91039.8406804571, + "2": 91403.8369707856, + "3": 91770.5359615407, + "4": 92139.0542427345, + "5": 92508.5040215493, + "6": 92877.9952611079, + "7": 93246.6378246507, + "8": 93613.5436199527, + "9": 93977.8287388161, + "10": 94338.6155864816, + "11": 94695.0349958323, + "12": 95046.2283212922, + "13": 95391.3495073796, + "14": 95729.5671269279, + "15": 96060.0663840674, + "16": 96382.0510771393, + "17": 96694.7455168155, + "18": 96997.396394802, + "19": 97289.2745986246, + "20": 97569.6769681249, + "21": 97837.9279894342, + "22": 98093.3814223464, + "23": 98335.4218571664, + "24": 98563.4661972865, + "25": 98776.9650639171, + "26": 98975.404119588, + "27": 99158.3053072332, + "28": 99325.228001871, + "29": 99475.7700721097, + "30": 99609.5688489159, + "31": 99726.3019993171, + "32": 99825.6883029292, + "33": 99907.4883294416, + "34": 99971.5050154257, + "35": 100017.5841390777, + "36": 100045.6146917527, + "37": 100055.529145394, + "38": 100047.3036152146, + "39": 100020.9579172368, + "40": 99976.5555205552, + "41": 99914.203394432, + "42": 99834.0517506009, + "43": 99736.293681392, + "44": 99621.1646945569, + "45": 99488.9421459088, + "46": 99339.9445711488, + "47": 99174.5309184851, + "48": 98993.0996838956, + "49": 98796.0879511157, + "50": 98583.9703386661, + "51": 98357.2578564544, + "52": 98116.4966747079, + "53": 97862.2668082017, + "54": 97595.180718952, + "55": 97315.8818407418, + "56": 97025.0430290325, + "57": 96723.3649399959, + "58": 96411.5743425723, + "59": 96090.4223676207, + "60": 95760.6826983782, + "61": 95423.1497065899, + "62": 95078.6365387974, + "63": 94727.9731573976, + "64": 94372.0043411908, + "65": 94011.5876502346, + "66": 93647.591359906, + "67": 93280.892369151, + "68": 92912.3740879571, + "69": 92542.9243091424, + "70": 92173.4330695837, + "71": 91804.790506041, + "72": 91437.8847107389, + "73": 91073.5995918756, + "74": 90712.8127442101, + "75": 90356.3933348594, + "76": 90005.2000093995, + "77": 89660.0788233121, + "78": 89321.8612037638, + "79": 88991.3619466243, + "80": 88669.3772535523, + "81": 88356.6828138761, + "82": 88054.0319358897, + "83": 87762.153732067, + "84": 87481.7513625668, + "85": 87213.5003412574, + "86": 86958.0469083452, + "87": 86716.0064735253, + "88": 86487.9621334051, + "89": 86274.4632667746, + "90": 86076.0242111036, + "91": 85893.1230234585, + "92": 85726.2003288206, + "93": 85575.658258582, + "94": 85441.8594817757, + "95": 85325.1263313746, + "96": 85225.7400277625, + "97": 85143.94000125, + "98": 85079.923315266, + "99": 85033.844191614, + "100": 85005.813638939, + "101": 84995.8991852976, + "102": 85004.1247154771, + "103": 85030.4704134548, + "104": 85074.8728101365, + "105": 85137.2249362596, + "106": 85217.3765800908, + "107": 85315.1346492996, + "108": 85430.2636361348, + "109": 85562.4861847829, + "110": 85711.4837595429, + "111": 85876.8974122066, + "112": 86058.3286467961, + "113": 86255.340379576, + "114": 86467.4579920255, + "115": 86694.1704742373, + "116": 86934.9316559837, + "117": 87189.16152249, + "118": 87456.2476117397, + "119": 87735.5464899498, + "120": 88026.3853016592, + "121": 88328.0633906958, + "122": 88639.8539881193, + "123": 88961.0059630709, + "124": 89290.7456323134, + "125": 89628.2786241017, + "126": 89972.7917918943, + "127": 90323.4551732941, + "128": 90679.4239895009 + } +} diff --git a/tests/test_interface/CMakeLists.txt b/tests/test_interface/CMakeLists.txt index 58ae0f827..e5ff9d3bb 100644 --- a/tests/test_interface/CMakeLists.txt +++ b/tests/test_interface/CMakeLists.txt @@ -6,3 +6,4 @@ project(svZeroDSolverInterfaceTests) add_subdirectory("test_01/") add_subdirectory("test_02/") add_subdirectory("test_03/") +add_subdirectory("test_04/") diff --git a/tests/test_interface/test_04/CMakeLists.txt b/tests/test_interface/test_04/CMakeLists.txt new file mode 100644 index 000000000..8d7f81d37 --- /dev/null +++ b/tests/test_interface/test_04/CMakeLists.txt @@ -0,0 +1,2 @@ +add_executable(svZeroD_interface_test04 ../LPNSolverInterface/LPNSolverInterface.cpp main.cpp) +target_link_libraries(svZeroD_interface_test04 ${CMAKE_DL_LIBS}) diff --git a/tests/test_interface/test_04/main.cpp b/tests/test_interface/test_04/main.cpp new file mode 100644 index 000000000..4c4cedbda --- /dev/null +++ b/tests/test_interface/test_04/main.cpp @@ -0,0 +1,104 @@ +// Test coupling trial/accept semantics for IMPEDANCE BC persistent memory. + +#include +#include +#include +#include +#include +#include + +#include "../LPNSolverInterface/LPNSolverInterface.h" +namespace fs = std::filesystem; + +int main(int argc, char** argv) { + LPNSolverInterface interface; + + if (argc != 3) { + throw std::runtime_error( + "Usage: svZeroD_interface_test04 " + ""); + } + + fs::path build_dir = argv[1]; + fs::path iface_dir = build_dir / "src" / "interface"; + fs::path lib_so = iface_dir / "libsvzero_interface.so"; + fs::path lib_dylib = iface_dir / "libsvzero_interface.dylib"; + fs::path lib_dll = iface_dir / "libsvzero_interface.dll"; + + if (fs::exists(lib_so)) { + interface.load_library(lib_so.string()); + } else if (fs::exists(lib_dylib)) { + interface.load_library(lib_dylib.string()); + } else if (fs::exists(lib_dll)) { + interface.load_library(lib_dll.string()); + } else { + throw std::runtime_error("Could not find shared libraries " + + lib_so.string() + " or " + lib_dylib.string() + + " or " + lib_dll.string() + " !"); + } + + interface.initialize(std::string(argv[2])); + interface.set_external_step_size(0.001); + + std::vector y0(interface.system_size_, 0.0); + std::vector ydot0(interface.system_size_, 0.0); + interface.return_y(y0); + interface.return_ydot(ydot0); + + std::vector params = {2.0, 0.0, 0.001, 1.0, 2.0}; + interface.update_block_params("FLOW_COUPLING", params); + + std::vector t(interface.num_output_steps_, 0.0); + std::vector sol1(interface.num_output_steps_ * interface.system_size_, + 0.0); + std::vector sol2(interface.num_output_steps_ * interface.system_size_, + 0.0); + std::vector sol3(interface.num_output_steps_ * interface.system_size_, + 0.0); + + int error_code = 0; + + interface.update_state(y0, ydot0); + interface.run_simulation(0.0, t, sol1, error_code); + if (error_code != 0) { + throw std::runtime_error("First trial run failed."); + } + + interface.update_state(y0, ydot0); + interface.run_simulation(0.0, t, sol2, error_code); + if (error_code != 0) { + throw std::runtime_error("Second trial run failed."); + } + + double max_trial_diff = 0.0; + for (size_t i = 0; i < sol1.size(); i++) { + max_trial_diff = std::max(max_trial_diff, std::abs(sol1[i] - sol2[i])); + } + if (max_trial_diff > 1.0e-12) { + throw std::runtime_error( + "Repeated trial runs from same committed state are not deterministic " + "for IMPEDANCE BC."); + } + + // Mark accepted state (and persistent memory) using current interface state. + std::vector ydot_commit(interface.system_size_, 0.0); + interface.return_ydot(ydot_commit); + + interface.update_state(y0, ydot0); + interface.run_simulation(0.0, t, sol3, error_code); + if (error_code != 0) { + throw std::runtime_error("Post-commit run failed."); + } + + double max_commit_diff = 0.0; + for (size_t i = 0; i < sol1.size(); i++) { + max_commit_diff = std::max(max_commit_diff, std::abs(sol3[i] - sol1[i])); + } + if (max_commit_diff < 1.0e-8) { + throw std::runtime_error( + "Persistent state commit had no observable effect for IMPEDANCE BC."); + } + + std::cout << "test_04 passed" << std::endl; + return 0; +} diff --git a/tests/test_interface/test_04/svzerod_3Dcoupling_impedance.json b/tests/test_interface/test_04/svzerod_3Dcoupling_impedance.json new file mode 100644 index 000000000..c126993b7 --- /dev/null +++ b/tests/test_interface/test_04/svzerod_3Dcoupling_impedance.json @@ -0,0 +1,36 @@ +{ + "simulation_parameters": { + "cardiac_period": 0.004, + "coupled_simulation": true, + "number_of_time_pts": 2, + "external_step_size": 0.001, + "output_all_cycles": true, + "steady_initial": false, + "absolute_tolerance": 1e-10 + }, + "boundary_conditions": [ + { + "bc_name": "IMP", + "bc_type": "IMPEDANCE", + "bc_values": { + "Pd": 5.0, + "z": [400.0, 120.0, 40.0, 10.0] + } + } + ], + "external_solver_coupling_blocks": [ + { + "name": "FLOW_COUPLING", + "type": "FLOW", + "location": "inlet", + "connected_block": "IMP", + "periodic": false, + "values": { + "t": [0.0, 0.004], + "Q": [1.0, 1.0] + } + } + ], + "junctions": [], + "vessels": [] +} diff --git a/tests/test_solver.py b/tests/test_solver.py index 8a7443835..08898740c 100644 --- a/tests/test_solver.py +++ b/tests/test_solver.py @@ -1,8 +1,6 @@ import os import json import numpy as np -import os -import json import pandas as pd import pytest @@ -54,6 +52,7 @@ 'piecewise_Chamber_and_Valve.json', 'closed_loop_two_hill.json', 'pulsatileFlow_CRL.json', + 'pulsatileFlow_R_impedance.json', 'pulsatileFlow_R_coronary_varres.json', 'closedLoopHeart_singleVessel_decomposed.json' ])