diff --git a/CMakeLists.txt b/CMakeLists.txt index 7e0468eed..ad4f94367 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -133,6 +133,25 @@ FetchContent_Declare( FetchContent_MakeAvailable(pybind11_json) +# ----------------------------------------------------------------------------- +# Fetch ExprTk +# ----------------------------------------------------------------------------- +# ExprTk is a header-only mathematical expression parser used by Parameter. +# +FetchContent_Declare( + exprtk + GIT_REPOSITORY https://github.com/ArashPartow/exprtk.git + GIT_TAG master + GIT_SHALLOW TRUE + GIT_PROGRESS TRUE +) + +FetchContent_MakeAvailable(exprtk) + +add_library(exprtk INTERFACE) +add_library(exprtk::exprtk ALIAS exprtk) +target_include_directories(exprtk INTERFACE ${exprtk_SOURCE_DIR}) + # ----------------------------------------------------------------------------- # Set executables for each application. # ----------------------------------------------------------------------------- diff --git a/src/model/Block.h b/src/model/Block.h index 849aeb956..14c7e3d13 100644 --- a/src/model/Block.h +++ b/src/model/Block.h @@ -25,7 +25,7 @@ * to the global system. */ struct TripletsContributions { - TripletsContributions() {}; + TripletsContributions() = default; /** * @brief Set the number of triplets that the element contributes * to the global system. @@ -33,7 +33,7 @@ struct TripletsContributions { * @param E Contributions to E matrix * @param D Contributions to dC/dy matrix */ - TripletsContributions(int F, int E, int D) : F(F), E(E), D(D) {}; + TripletsContributions(int F, int E, int D) : F(F), E(E), D(D) {} /** * @brief Set the number of triplets that the element contributes * to the global system. diff --git a/src/model/CMakeLists.txt b/src/model/CMakeLists.txt index 88ba98d5f..16d794828 100644 --- a/src/model/CMakeLists.txt +++ b/src/model/CMakeLists.txt @@ -72,6 +72,11 @@ set(HDRS add_library(${lib} OBJECT ${CXXSRCS} ${HDRS}) +# exprtk.hpp generates too many COFF sections for the default MSVC object format +if(MSVC) + set_source_files_properties(Parameter.cpp PROPERTIES COMPILE_FLAGS "/bigobj") +endif() + target_include_directories( ${lib} PUBLIC ${CMAKE_SOURCE_DIR}/src/algebra ${CMAKE_SOURCE_DIR}/src/model @@ -80,4 +85,4 @@ target_include_directories( ${lib} PUBLIC target_link_libraries( ${lib} Eigen3::Eigen ) target_link_libraries( ${lib} nlohmann_json::nlohmann_json ) - +target_link_libraries( ${lib} exprtk::exprtk ) diff --git a/src/model/FlowReferenceBC.h b/src/model/FlowReferenceBC.h index f626c710e..ecd7bdae4 100644 --- a/src/model/FlowReferenceBC.h +++ b/src/model/FlowReferenceBC.h @@ -84,8 +84,9 @@ class FlowReferenceBC : public Block { */ FlowReferenceBC(int id, Model* model) : Block(id, model, BlockType::flow_bc, BlockClass::boundary_condition, - {{"t", InputParameter(false, true)}, - {"Q", InputParameter(false, true)}}) {} + {{"t", InputParameter(true, true)}, + {"Q", InputParameter(true, true)}, + {"fn", InputParameter(true, false, false, true)}}) {} /** * @brief Set up the degrees of freedom (DOF) of the block diff --git a/src/model/Model.cpp b/src/model/Model.cpp index 5b7fb5066..d338443b5 100644 --- a/src/model/Model.cpp +++ b/src/model/Model.cpp @@ -158,6 +158,13 @@ int Model::add_parameter(const std::vector& times, return parameter_count++; } +int Model::add_parameter(const std::string expression_string) { + auto param = Parameter(parameter_count, expression_string); + parameter_values.push_back(param.get(0.0)); + parameters.push_back(std::move(param)); + return parameter_count++; +} + Parameter* Model::get_parameter(int param_id) { return ¶meters[param_id]; } double Model::get_parameter_value(int param_id) const { @@ -181,6 +188,10 @@ void Model::finalize() { for (auto& block : blocks) { block->setup_model_dependent_params(); } + + if (cardiac_cycle_period < 0.0) { + cardiac_cycle_period = 1.0; + } } int Model::get_num_blocks(bool internal) const { diff --git a/src/model/Model.h b/src/model/Model.h index 9412e932b..45a630d90 100644 --- a/src/model/Model.h +++ b/src/model/Model.h @@ -186,6 +186,14 @@ class Model { int add_parameter(const std::vector& times, const std::vector& values, bool periodic = true); + /** + * @brief Add a function-based model parameter + * + * @param expression_string String with time-dependent function expression + * @return int Global ID of the parameter + */ + int add_parameter(const std::string expression_string); + /** * @brief Get a parameter by its global ID * diff --git a/src/model/Parameter.cpp b/src/model/Parameter.cpp index 3a43154fd..cfe5bc6f0 100644 --- a/src/model/Parameter.cpp +++ b/src/model/Parameter.cpp @@ -2,6 +2,36 @@ // University of California, and others. SPDX-License-Identifier: BSD-3-Clause #include "Parameter.h" +#include +#include + +// Suppress -Wshadow warnings from the exprtk third-party header +#if defined(__clang__) +#pragma clang diagnostic push +#pragma clang diagnostic ignored "-Wshadow" +#elif defined(__GNUC__) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wshadow" +#endif +#include +#if defined(__clang__) +#pragma clang diagnostic pop +#elif defined(__GNUC__) +#pragma GCC diagnostic pop +#endif + +/// PIMPL implementation holding the exprtk expression state +struct Parameter::ExprtkImpl { + exprtk::symbol_table + symbol_table; ///< Symbol table mapping 't' to current time + exprtk::expression expression; ///< Compiled expression object + double* time_value = nullptr; ///< Pointer into the symbol table for 't' +}; + +Parameter::~Parameter() = default; +Parameter::Parameter(Parameter&&) noexcept = default; +Parameter& Parameter::operator=(Parameter&&) noexcept = default; + Parameter::Parameter(int id, double value) { this->id = id; update(value); @@ -14,6 +44,12 @@ Parameter::Parameter(int id, const std::vector& times, update(times, values); } +Parameter::Parameter(int id, const std::string expression_string) { + this->id = id; + this->expression_string = expression_string; + update(expression_string); +} + void Parameter::update(double update_value) { is_constant = true; is_periodic = true; @@ -35,6 +71,48 @@ void Parameter::update(const std::vector& update_times, } } +void Parameter::update(const std::string update_string) { + is_function = true; + is_constant = false; + expression_string = update_string; + + exprtk_impl_ = std::make_unique(); + + if (!exprtk_impl_->symbol_table.create_variable("t")) { + throw std::runtime_error( + "Error failed to create time_value in symbol_table."); + } + + exprtk_impl_->symbol_table.add_constants(); + + exprtk_impl_->time_value = + &exprtk_impl_->symbol_table.get_variable("t")->ref(); + exprtk_impl_->expression.register_symbol_table(exprtk_impl_->symbol_table); + + exprtk::parser parser; + + if (!parser.compile(expression_string, exprtk_impl_->expression)) { + is_function = false; + exprtk_impl_.reset(); + + printf("Error: %s\tExpression: %s\n", parser.error().c_str(), + expression_string.c_str()); + + for (std::size_t i = 0; i < parser.error_count(); ++i) { + const auto error = parser.get_error(i); + + printf( + "Error: %02d Position: %02d Type: [%14s] Msg: %s\tExpression: %s\n", + static_cast(i), + static_cast(error.token.position), + exprtk::parser_error::to_str(error.mode).c_str(), + error.diagnostic.c_str(), expression_string.c_str()); + } + throw std::runtime_error( + "Error when compiling the function provided in 'fn'."); + } +} + double Parameter::get(double time) { // Return the constant value if parameter is constant if (is_constant) { @@ -51,6 +129,13 @@ double Parameter::get(double time) { rtime = time; } + if (is_function) { + assert(exprtk_impl_ != nullptr); + assert(exprtk_impl_->time_value != nullptr); + *exprtk_impl_->time_value = time; + return exprtk_impl_->expression.value(); + } + // Determine the lower and upper element for interpolation auto i = lower_bound(times.begin(), times.end(), rtime); int k = i - times.begin(); diff --git a/src/model/Parameter.h b/src/model/Parameter.h index 515862ca9..6542c0ab0 100644 --- a/src/model/Parameter.h +++ b/src/model/Parameter.h @@ -10,7 +10,9 @@ #include #include #include +#include #include +#include #include #include "DOFHandler.h" @@ -43,6 +45,30 @@ class Parameter { Parameter(int id, const std::vector& times, const std::vector& values, bool periodic = true); + /** + * @brief Construct a new Parameter object + * + * @param id Global ID of the parameter + * @param expression_string string with time-dependent function to get values + */ + Parameter(int id, const std::string expression_string); + + /** + * @brief Destructor + */ + ~Parameter(); + + /** + * @brief Move constructor + */ + Parameter(Parameter&&) noexcept; + + /** + * @brief Move assignment operator + * @return Reference to this object + */ + Parameter& operator=(Parameter&&) noexcept; + int id; ///< Global ID of the parameter std::vector times; ///< Time steps if parameter is time-dependent std::vector values; ///< Values if parameter is time-dependent @@ -53,6 +79,9 @@ class Parameter { bool is_constant; ///< Bool value indicating if the parameter is constant bool is_periodic; ///< Bool value indicating if the parameter is periodic ///< with the cardiac cycle + bool is_function = false; ///< Bool value indicating if the parameter is a + ///< function + std::string expression_string; ///< String with value function /** * @brief Update the parameter @@ -70,6 +99,13 @@ class Parameter { void update(const std::vector& times, const std::vector& values); + /** + * @brief Update the parameter + * + * @param update_string String with function + */ + void update(const std::string update_string); + /** * @brief Get the parameter value at the specified time. * @@ -92,6 +128,8 @@ class Parameter { private: bool steady_converted = false; + struct ExprtkImpl; ///< Forward declared PIMPL struct + std::unique_ptr exprtk_impl_; ///< PIMPL for exprtk types }; /** @@ -101,6 +139,7 @@ struct InputParameter { bool is_optional; ///< Is this parameter optional? bool is_array; ///< Is this parameter an array? bool is_number; ///< Is this parameter a number? + bool is_function; ///< Is this parameter a time-dependent function? double default_val; ///< Default value (if parameter is optional) /** @@ -109,13 +148,16 @@ struct InputParameter { * @param is_optional Is this parameter optional? * @param is_array Is this parameter an array? * @param is_number Is this parameter a number? + * @param is_function Is this parameter a time-dependent function? * @param default_val Default value (if parameter is optional) */ InputParameter(bool is_optional = false, bool is_array = false, - bool is_number = true, double default_val = 0.0) + bool is_number = true, bool is_function = false, + double default_val = 0.0) : is_optional(is_optional), is_array(is_array), is_number(is_number), + is_function(is_function), default_val(default_val) {} }; diff --git a/src/model/PressureReferenceBC.h b/src/model/PressureReferenceBC.h index 7a5bda0e6..73dfc9f94 100644 --- a/src/model/PressureReferenceBC.h +++ b/src/model/PressureReferenceBC.h @@ -85,8 +85,9 @@ class PressureReferenceBC : public Block { */ PressureReferenceBC(int id, Model* model) : Block(id, model, BlockType::pressure_bc, BlockClass::boundary_condition, - {{"t", InputParameter(false, true)}, - {"P", InputParameter(false, true)}}) {} + {{"t", InputParameter(true, true)}, + {"P", InputParameter(true, true)}, + {"fn", InputParameter(true, false, false, true)}}) {} /** * @brief Set up the degrees of freedom (DOF) of the block diff --git a/src/solve/SimulationParameters.cpp b/src/solve/SimulationParameters.cpp index 036609994..7d198aa5c 100644 --- a/src/solve/SimulationParameters.cpp +++ b/src/solve/SimulationParameters.cpp @@ -30,6 +30,20 @@ bool get_param_vector(const nlohmann::json& data, const std::string& name, return false; } +bool get_param_string(const nlohmann::json& data, const std::string& name, + const InputParameter& param, std::string& val) { + if (data.contains(name)) { + val = data[name].get(); + } else { + if (param.is_optional) { + val = ""; + } else { + return true; + } + } + return false; +} + bool has_parameter( const std::vector>& params, const std::string& name) { @@ -85,21 +99,32 @@ int generate_block(Model& model, const nlohmann::json& block_params_json, continue; } - // Skip reading parameters that are not a number - if (!block_param.second.is_number) { + // Skip reading parameters that are not a number and not a function + if (!block_param.second.is_number && !block_param.second.is_function) { continue; } - // Get vector parameter - if (block_param.second.is_array) { + if (block_param.second.is_function) { + std::string expression_string; + err = get_param_string(block_params_json, block_param.first, + block_param.second, expression_string); + if (expression_string.length() <= 1) { + continue; + } + new_id = model.add_parameter(expression_string); + } else if (block_param.second.is_array) { // Get parameter vector std::vector val; err = get_param_vector(block_params_json, block_param.first, block_param.second, val); if (err) { - throw std::runtime_error("Array parameter " + block_param.first + - " is mandatory in " + block_type + - " block " + static_cast(name)); + if (!block_param.second.is_optional) { + throw std::runtime_error( + "Array parameter " + block_param.first + " is mandatory in " + + block_type + " block " + static_cast(name)); + } else { + continue; + } } // Get time vector @@ -107,9 +132,13 @@ int generate_block(Model& model, const nlohmann::json& block_params_json, std::vector time; err = get_param_vector(block_params_json, "t", t_param, time); if (err) { - throw std::runtime_error("Array parameter t is mandatory in " + - block_type + " block " + - static_cast(name)); + if (!block_param.second.is_optional) { + throw std::runtime_error("Array parameter t is mandatory in " + + block_type + " block " + + static_cast(name)); + } else { + continue; + } } // Add parameters to model diff --git a/tests/cases/dirgraph-results/timeDep_Flow_directed_graph.dot b/tests/cases/dirgraph-results/timeDep_Flow_directed_graph.dot new file mode 100644 index 000000000..93378161d --- /dev/null +++ b/tests/cases/dirgraph-results/timeDep_Flow_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/dirgraph-results/timeDep_Pressure_directed_graph.dot b/tests/cases/dirgraph-results/timeDep_Pressure_directed_graph.dot new file mode 100644 index 000000000..93378161d --- /dev/null +++ b/tests/cases/dirgraph-results/timeDep_Pressure_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/timeDep_Flow.json b/tests/cases/timeDep_Flow.json new file mode 100644 index 000000000..26516cdf2 --- /dev/null +++ b/tests/cases/timeDep_Flow.json @@ -0,0 +1,52 @@ +{ + "description": { + "description of test case": "pulsatile flow -> -> RC -> R", + "analytical results": [ + "Boundary conditions:", + "inlet:", + "flow = 2*pi*cos(2*pi*t)", + "outlet:", + "Resistance: Rd = 1000, Pd = 0" + ] + }, + "boundary_conditions": [ + { + "bc_name": "INLET", + "bc_type": "FLOW", + "bc_values": { + "fn": "2.0 * (4*ATAN(1.)) * COS(2.0 * (4*ATAN(1.)) * t)" + } + }, + { + "bc_name": "OUTLET", + "bc_type": "RESISTANCE", + "bc_values": { + "Pd": 0.0, + "R": 1000.0 + } + } + ], + "junctions": [], + "simulation_parameters": { + "number_of_cardiac_cycles": 3, + "number_of_time_pts_per_cardiac_cycle": 100, + "output_all_cycles": true, + "output_variable_based": true + }, + "vessels": [ + { + "boundary_conditions": { + "inlet": "INLET", + "outlet": "OUTLET" + }, + "vessel_id": 0, + "vessel_length": 10.0, + "vessel_name": "branch0_seg0", + "zero_d_element_type": "BloodVessel", + "zero_d_element_values": { + "C": 0.0001, + "R_poiseuille": 100.0 + } + } + ] +} diff --git a/tests/cases/timeDep_Pressure.json b/tests/cases/timeDep_Pressure.json new file mode 100644 index 000000000..bd1de7af9 --- /dev/null +++ b/tests/cases/timeDep_Pressure.json @@ -0,0 +1,52 @@ +{ + "description": { + "description of test case": "pulsatile pressure -> Blood Vessel -> R", + "analytical results": [ + "Boundary conditions:", + "inlet:", + "Pressure = 2*pi*cos(2*pi*t)", + "outlet:", + "Resistance: R = 1000, Pd = 0" + ] + }, + "boundary_conditions": [ + { + "bc_name": "INLET", + "bc_type": "PRESSURE", + "bc_values": { + "fn": "2.0 * (4*ATAN(1.)) * COS(2.0 * (4*ATAN(1.)) * t)" + } + }, + { + "bc_name": "OUTLET", + "bc_type": "RESISTANCE", + "bc_values": { + "Pd": 0.0, + "R": 1000.0 + } + } + ], + "junctions": [], + "simulation_parameters": { + "number_of_cardiac_cycles": 3, + "number_of_time_pts_per_cardiac_cycle": 100, + "output_all_cycles": true, + "output_variable_based": true + }, + "vessels": [ + { + "boundary_conditions": { + "inlet": "INLET", + "outlet": "OUTLET" + }, + "vessel_id": 0, + "vessel_length": 10.0, + "vessel_name": "branch0_seg0", + "zero_d_element_type": "BloodVessel", + "zero_d_element_values": { + "C": 0.0001, + "R_poiseuille": 100.0 + } + } + ] +} diff --git a/tests/test_io.py b/tests/test_io.py index d40136259..9e33518f5 100644 --- a/tests/test_io.py +++ b/tests/test_io.py @@ -1,6 +1,10 @@ import numpy as np +import pytest + +import pysvzerod + +from .utils import run_test_case_by_name, RTOL_FLOW, RTOL_PRES, get_test_case_by_name -from .utils import run_test_case_by_name, RTOL_FLOW, RTOL_PRES # use coarse absolute tolerances for gradient calculation # we're comparing gradients from gen-alpha in svZeroDSolver with central differences in np.gradient @@ -121,3 +125,69 @@ def test_pulsatile_flow_r_rcr_mean_derivative_variable(): rt = np.mean(np.gradient(res_time[f + "_" + k], dt)) rm = res_mean["ydot"][res_mean["name"] == f + ":" + v] assert np.isclose(rt, rm, atol=ATOL_MEAN[f[0]]) + + +def test_time_dependent_block(): + # time-dependent results + res_flow = run_test_case_by_name("timeDep_Flow", output_variable_based=True) + res_press = run_test_case_by_name( + "timeDep_Pressure", output_variable_based=True + ) + + time = res_flow["time"][ + res_flow["name"] == "flow:INLET:branch0_seg0" + ].to_numpy() + flow = res_flow["y"][ + res_flow["name"] == "flow:INLET:branch0_seg0" + ].to_numpy() + pressure = res_press["y"][ + res_press["name"] == "pressure:INLET:branch0_seg0" + ].to_numpy() + + # compare time-dependent results to the calculated values + for i, t in enumerate(time): + # may be unsteady at the beginning + if i < 15: + continue + calc_val = 2.0 * (4 * np.arctan(1.0)) * np.cos( + 2.0 * (4 * np.arctan(1.0)) * t + ) + assert np.isclose(flow[i], calc_val, rtol=0.003) + assert np.isclose(pressure[i], calc_val, rtol=0.003) + + + +def test_invalid_fn_expression(): + """Invalid exprtk expression string triggers compile-error path in Parameter.""" + + config = get_test_case_by_name("timeDep_Flow") + + # input invalid expression + config["boundary_conditions"][0]["bc_values"] = {"fn": "sin(t"} + + with pytest.raises(RuntimeError): + pysvzerod.simulate(config) + + +def test_missing_mandatory_array_param(): + """Missing mandatory array parameter raises RuntimeError from SimulationParameters.""" + + config = get_test_case_by_name("pulsatileFlow_R_coronary") + + # remove required t vector for array parameter Pim + del config["boundary_conditions"][1]["bc_values"]["Pim"] + + with pytest.raises(RuntimeError): + pysvzerod.simulate(config) + + +def test_missing_t_vector_for_array_param(): + """Mandatory array param present but missing t vector raises RuntimeError.""" + + config = get_test_case_by_name("pulsatileFlow_R_coronary") + + # remove required t vector for array parameter Pim + del config["boundary_conditions"][1]["bc_values"]["t"] + + with pytest.raises(RuntimeError): + pysvzerod.simulate(config) diff --git a/tests/utils.py b/tests/utils.py index 41fc68243..093807432 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -189,3 +189,11 @@ def get_result(result_array, field, branch, time_step): """ "Get results at specific field, branch, branch_node and time step.""" # extract result return result_array[field][branch][time_step] + + +def get_test_case_by_name(test_case_name): + """Get test case configuration by name.""" + testfile = os.path.join(this_file_dir, "cases", test_case_name + ".json") + with open(testfile) as ff: + config = json.load(ff) + return config \ No newline at end of file