diff --git a/.github/codecov.yml b/.github/codecov.yml new file mode 100644 index 000000000..529e569a6 --- /dev/null +++ b/.github/codecov.yml @@ -0,0 +1,11 @@ +# configuration for code coverage testing with codecov +coverage: + status: + project: + default: + informational: true + patch: + default: + informational: true +codecov: + token: ${{ secrets.CODECOV_TOKEN }} \ No newline at end of file diff --git a/.github/workflows/codechecks.yml b/.github/workflows/codechecks.yml new file mode 100644 index 000000000..d797f9da0 --- /dev/null +++ b/.github/workflows/codechecks.yml @@ -0,0 +1,17 @@ +# This workflow checks for compliance with the Google C++ style guide. +name: Codechecks +on: [push, pull_request] +jobs: + clang-format: + runs-on: macos-latest + steps: + - uses: actions/checkout@v3 + - name: Install dependencies + run: | + brew install clang-format + - name: Run clang-format + run: | + mkdir Release + cd Release + cmake .. + make codecheck diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml new file mode 100644 index 000000000..1b87c90bf --- /dev/null +++ b/.github/workflows/documentation.yml @@ -0,0 +1,30 @@ +# This workflow builds and deploys the html documentation for svZeroDSolver. +name: Documentation +on: [push, pull_request] +permissions: + contents: write +jobs: + documentation: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - name: Make build directory + run: mkdir docs/build + - name: Build doxygen documentation + continue-on-error: false + uses: mattnotmitt/doxygen-action@edge + with: + working-directory: '.' + doxyfile-path: 'docs/Doxyfile' + enable-latex: true + - name: Save documentation + uses: actions/upload-artifact@v4 + with: + name: documentation + path: ./docs/build/html + - name: Deploy documentation + if: github.ref == 'refs/heads/master' + uses: peaceiris/actions-gh-pages@v3 + with: + github_token: ${{ secrets.GITHUB_TOKEN }} + publish_dir: ./docs/build/html diff --git a/.github/workflows/gui.yml b/.github/workflows/gui.yml new file mode 100644 index 000000000..064ba1d18 --- /dev/null +++ b/.github/workflows/gui.yml @@ -0,0 +1,49 @@ +# This workflow uses Cypress for end-to-end testing of the 0D model GUI. + +name: GUI-tests + +on: [push, pull_request] + +jobs: + cypress-run: + strategy: + matrix: + os: [ubuntu-latest, macos-latest] + fail-fast: false + runs-on: ${{ matrix.os }} + + steps: + - name: Checkout code + uses: actions/checkout@v3 + + - name: Set up Node.js + uses: actions/setup-node@v3 + with: + node-version: '20' + + - name: Install dependencies + working-directory: tests/cypress + run: npm install + + - name: Set up Python + uses: actions/setup-python@v4 + with: + python-version: '3.10' # Use the Python version compatible with your Flask app + + - name: Install Flask dependencies + run: | + python -m pip install --upgrade pip + python -m pip install flask + + - name: Start Flask Application + run: | + cd applications/svZeroDGUI + FLASK_APP=app.py flask run --host=0.0.0.0 --port=8902 & + env: + FLASK_ENV: development + + - name: Run Cypress tests + uses: cypress-io/github-action@v5 + with: + start: npm start + working-directory: tests/cypress diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml new file mode 100644 index 000000000..e52875d9a --- /dev/null +++ b/.github/workflows/test.yml @@ -0,0 +1,170 @@ +# This workflow builds and tests svZeroDSolver. It is built and tested on +# different versions of ubuntu and macOS. +name: Build and test +on: [push, pull_request] +jobs: + test: + runs-on: ${{ matrix.os }} + strategy: + matrix: + os: [ubuntu-22.04, ubuntu-latest, macos-13, macos-latest, windows-latest] + version: [13] # GCC version + fail-fast: false + env: + GCC_V: ${{ matrix.version }} + CONDA_ENV: zerod + PYTHONPATH: ${{ github.workspace }} + steps: + - uses: actions/checkout@v4 + - uses: conda-incubator/setup-miniconda@v3 + with: + auto-update-conda: true + activate-environment: ${{env.CONDA_ENV}} + python-version: "3.11.4" + - name: Install ubuntu dependencies + if: startsWith(matrix.os, 'ubuntu') + run: sudo apt update && sudo apt install build-essential cmake lcov + + - name: Install dependencies to get correct version numbers (Ubuntu) + if: startsWith(matrix.os, 'ubuntu') + run: conda install -c conda-forge libstdcxx-ng=${GCC_V} gcc=${GCC_V} + + - name: Install dependencies to get correct version numbers (MacOS) + if: startsWith(matrix.os, 'macos') + run: | + brew install gcc@${GCC_V} + ln -s /usr/local/bin/gcc-${GCC_V} /usr/local/bin/gcc + + - name: Install dependencies for windows + if: startsWith(matrix.os, 'windows') + shell: pwsh + run: | + choco install mingw --no-progress + conda install -y -c conda-forge cmake graphviz python-graphviz pydot + pip install --upgrade cmake-setuptools + + - name: Install POISX-like svZeroDSolver + if: ${{!startsWith(matrix.os, 'windows')}} + run: conda run pip install -e ".[dev]" + + - name: Install Windows svZeroDSolver + if: startsWith(matrix.os, 'windows') + shell: pwsh + run: | + $Env:CMAKE_GENERATOR = 'MinGW Makefiles' + Write-Host "→ Using CMAKE_GENERATOR = $Env:CMAKE_GENERATOR" + pip install --no-build-isolation -v .[dev] + pip show pysvzerod + + - name: Install Networkx + run: | + conda run pip install networkx + + - name: Test the build + run: | + cd tests + conda run pytest -v --durations=0 --ignore=test_dirgraph.py + + - name: Build using CMake for POISX-like Systems + if: ${{!startsWith(matrix.os, 'windows')}} + run: | + mkdir Release + cd Release + cmake -DCMAKE_BUILD_TYPE=Release -DENABLE_DISTRIBUTION=ON .. + make -j2 + + - name: Build using CMake for Windows Systems + if: startsWith(matrix.os, 'windows') + shell: pwsh + run: | + mkdir Release + cd Release + cmake -G "MinGW Makefiles" ` + -DCMAKE_BUILD_TYPE=Release ` + -DENABLE_DISTRIBUTION=ON ` + .. + cmake --build . --parallel 2 + + - name: Test interface POISX-like Systems + if: ${{!startsWith(matrix.os, 'windows')}} + run: | + cd tests/test_interface + mkdir build_tests + cd build_tests + cmake ../ + make -j2 + cd test_01 + ./svZeroD_interface_test01 ../../../../Release ../../test_01/svzerod_3Dcoupling.json + cd ../test_02 + ./svZeroD_interface_test02 ../../../../Release ../../test_02/svzerod_tuned.json + + - name: Test interface Windows Systems + if: startsWith(matrix.os, 'windows') + shell: pwsh + run: | + cd tests/test_interface + mkdir build_tests + cd build_tests + cmake -G "MinGW Makefiles" .. + cmake --build . --parallel 2 + cd test_01 + ./svZeroD_interface_test01.exe ` + ../../../../Release ` + ../../test_01/svzerod_3Dcoupling.json + + cd ../test_02 + ./svZeroD_interface_test02 ` + ../../../../Release ` + ../../test_02/svzerod_tuned.json + + - name: Generate code coverage + if: startsWith(matrix.os, 'ubuntu-22.04') + run: | + cd Release + cmake -DENABLE_COVERAGE=ON .. + make -j2 + cd ../tests + conda run pytest -v --durations=0 --coverage --ignore=test_dirgraph.py + cd ../Release + make coverage + + - name: Save coverage report + if: startsWith(matrix.os, 'ubuntu-22.04') + uses: actions/upload-artifact@v4 + with: + name: coverage_report + path: Release/coverage + + - name: Upload coverage reports to Codecov + if: startsWith(matrix.os, 'ubuntu-22.04') + uses: codecov/codecov-action@v4 + env: + CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} + + - name: Build installer POISX-like Systems + if: ${{!startsWith(matrix.os, 'windows')}} + run: | + cd Release + cpack + cp distribution/svZeroDSolver_* .. + + - name: Install NSIS + if: startsWith(matrix.os, 'windows') + run: choco install nsis -y + + - name: Build installer Windows Systems + if: startsWith(matrix.os, 'windows') + shell: pwsh + run: | + # NSIS will be at C:\Program Files (x86)\NSIS\Bin + echo "C:\Program Files (x86)\NSIS\Bin" | Out-File -FilePath $env:GITHUB_PATH -Encoding utf8 -Append + cd Release + cpack -C Release + Copy-Item distribution\svZeroDSolver_* -Destination ..\ + + - name: Upload installer + uses: actions/upload-artifact@v4 + with: + name: ${{ matrix.os }} installer + path: svZeroDSolver_* + if-no-files-found: error diff --git a/.github/workflows/test_visualization.yml b/.github/workflows/test_visualization.yml new file mode 100644 index 000000000..b895662ce --- /dev/null +++ b/.github/workflows/test_visualization.yml @@ -0,0 +1,40 @@ +name: Test Visualization Application + +on: [push, pull_request] + +jobs: + test: + strategy: + matrix: + os: [ubuntu-latest, macos-latest] + fail-fast: false + runs-on: ${{ matrix.os }} + + steps: + - name: Checkout Code + uses: actions/checkout@v3 + + - name: Set up Python + uses: actions/setup-python@v4 + with: + python-version: '3.11' + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install -r tests/requirements.txt + + - name: Install Graphviz + if: runner.os == 'Linux' + run: | + sudo apt-get update + sudo apt-get install -y graphviz + + - name: Install Graphviz (macOS) + if: runner.os == 'macOS' + run: | + brew install graphviz + + - name: Run tests + run: | + pytest tests/test_dirgraph.py diff --git a/CMakeLists.txt b/CMakeLists.txt index 1cf48da8b..f611be68e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -65,10 +65,13 @@ endif() # ----------------------------------------------------------------------------- # Eigen is a header-only C++ template library for linear algebra. # +# Tell FetchContent not to re-check for updates once Eigen is downloaded +set(FETCHCONTENT_UPDATES_DISCONNECTED TRUE) + FetchContent_Declare( Eigen GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git - GIT_TAG master + GIT_TAG 3.4.0 GIT_SHALLOW TRUE GIT_PROGRESS TRUE) diff --git a/applications/pysvzerod.cpp b/applications/pysvzerod.cpp index 15bc27c13..281d42428 100644 --- a/applications/pysvzerod.cpp +++ b/applications/pysvzerod.cpp @@ -1,5 +1,5 @@ -// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. -// SPDX-License-Identifier: BSD-3-Clause +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause /** * @file pysvzerod.cpp * @brief Python interface for svZeroDSolver @@ -65,9 +65,8 @@ PYBIND11_MODULE(pysvzerod, m) { py::module_ sys = py::module_::import("sys"); auto argv = sys.attr("argv").cast>(); if (argv.size() != 3) { - std::cout - << "Usage: svzerodsolver path/to/config.json path/to/output.csv" - << std::endl; + std::cout << "Usage: svzerodsolver path/to/config.json path/to/output.csv" + << std::endl; exit(1); } std::ifstream ifs(argv[1]); diff --git a/applications/svzerodcalibrator.cpp b/applications/svzerodcalibrator.cpp index f9dc0479a..f9505807c 100644 --- a/applications/svzerodcalibrator.cpp +++ b/applications/svzerodcalibrator.cpp @@ -1,5 +1,5 @@ -// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. -// SPDX-License-Identifier: BSD-3-Clause +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause /** * @file svzerodcalibrator.cpp * @brief Main routine for svZeroDCalibrator @@ -26,7 +26,8 @@ int main(int argc, char* argv[]) { std::ifstream input_file(input_file_name); if (!input_file.is_open()) { - std::cerr << "[svzerodcalibrator] Error: The input file '" << input_file_name << "' cannot be opened." << std::endl; + std::cerr << "[svzerodcalibrator] Error: The input file '" + << input_file_name << "' cannot be opened." << std::endl; return 1; } @@ -36,8 +37,10 @@ int main(int argc, char* argv[]) { try { output_config = calibrate(config); } catch (const nlohmann::json::parse_error& e) { - std::cerr << "[svzerodcalibrator] Error: The input file '" << input_file_name - << "' does not have the parameters needed by the calibrate program." << std::endl; + std::cerr + << "[svzerodcalibrator] Error: The input file '" << input_file_name + << "' does not have the parameters needed by the calibrate program." + << std::endl; return 1; } @@ -45,7 +48,8 @@ int main(int argc, char* argv[]) { std::ofstream out_file(output_file_name); if (!out_file.is_open()) { - std::cerr << "[svzerodcalibrator] Error: The output file '" << output_file_name << "' cannot be opened." << std::endl; + std::cerr << "[svzerodcalibrator] Error: The output file '" + << output_file_name << "' cannot be opened." << std::endl; return 1; } diff --git a/applications/svzerodsolver.cpp b/applications/svzerodsolver.cpp index b1d7d3448..873d6dc11 100644 --- a/applications/svzerodsolver.cpp +++ b/applications/svzerodsolver.cpp @@ -1,5 +1,5 @@ -// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. -// SPDX-License-Identifier: BSD-3-Clause +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause /** * @file svzerodsolver.cpp * @brief Main routine of svZeroDSolver @@ -30,7 +30,9 @@ int main(int argc, char* argv[]) { // Get input and output file name if (argc < 2 || argc > 3) { - throw std::runtime_error("Usage: svzerodsolver path/to/config.json [optional:path/to/output.csv]"); + throw std::runtime_error( + "Usage: svzerodsolver path/to/config.json " + "[optional:path/to/output.csv]"); } std::string input_file_name = argv[1]; @@ -56,23 +58,27 @@ int main(int argc, char* argv[]) { } output_file_name = output_file_path + "/output.csv"; - std::cout << "[svzerodsolver] Output will be written to '" << output_file_name << "'." << std::endl;; + std::cout << "[svzerodsolver] Output will be written to '" + << output_file_name << "'." << std::endl; + ; } std::ifstream input_file(input_file_name); if (!input_file.is_open()) { - std::cerr << "[svzerodsolver] Error: The input file '" << input_file_name << "' cannot be opened." << std::endl; + std::cerr << "[svzerodsolver] Error: The input file '" << input_file_name + << "' cannot be opened." << std::endl; return 1; } nlohmann::json config; - try { + try { config = nlohmann::json::parse(input_file); } catch (const nlohmann::json::parse_error& e) { - std::cout << "[svzerodsolver] Error: Parsing the input file '" << input_file_name << "' has failed." << std::endl; + std::cout << "[svzerodsolver] Error: Parsing the input file '" + << input_file_name << "' has failed." << std::endl; std::cout << "[svzerodsolver] Details of the parsing error: " << std::endl; std::cout << e.what() << std::endl; return 1; diff --git a/src/algebra/SparseSystem.cpp b/src/algebra/SparseSystem.cpp index c0e0b3156..689d449d5 100644 --- a/src/algebra/SparseSystem.cpp +++ b/src/algebra/SparseSystem.cpp @@ -27,7 +27,7 @@ void SparseSystem::clean() { // delete solver; } -void SparseSystem::reserve(Model *model) { +void SparseSystem::reserve(Model* model) { auto num_triplets = model->get_num_triplets(); F.reserve(num_triplets.F); E.reserve(num_triplets.E); @@ -56,8 +56,8 @@ void SparseSystem::reserve(Model *model) { } void SparseSystem::update_residual( - Eigen::Matrix &y, - Eigen::Matrix &ydot) { + Eigen::Matrix& y, + Eigen::Matrix& ydot) { residual.setZero(); residual -= C; residual.noalias() -= E * ydot; diff --git a/src/algebra/SparseSystem.h b/src/algebra/SparseSystem.h index c3efdd1a1..373a858f0 100644 --- a/src/algebra/SparseSystem.h +++ b/src/algebra/SparseSystem.h @@ -70,7 +70,7 @@ class SparseSystem { * * @param model The model to reserve space for in the system */ - void reserve(Model *model); + void reserve(Model* model); /** * @brief Update the residual of the system @@ -78,8 +78,8 @@ class SparseSystem { * @param y Vector of current solution quantities * @param ydot Derivate of y */ - void update_residual(Eigen::Matrix &y, - Eigen::Matrix &ydot); + void update_residual(Eigen::Matrix& y, + Eigen::Matrix& ydot); /** * @brief Update the jacobian of the system diff --git a/src/algebra/State.cpp b/src/algebra/State.cpp index dae11d364..03d498b61 100644 --- a/src/algebra/State.cpp +++ b/src/algebra/State.cpp @@ -12,7 +12,7 @@ State::State(int n) { State::~State() {} -State::State(const State &state) { +State::State(const State& state) { y = state.y; ydot = state.ydot; } diff --git a/src/algebra/State.h b/src/algebra/State.h index 8b58aa7ab..e9b286a8e 100644 --- a/src/algebra/State.h +++ b/src/algebra/State.h @@ -46,7 +46,7 @@ class State { * * @param state */ - State(const State &state); + State(const State& state); /** * @brief Construct a new State object and initilaize with all zeros. diff --git a/src/model/Block.cpp b/src/model/Block.cpp index c25fb24c5..7800fa10d 100644 --- a/src/model/Block.cpp +++ b/src/model/Block.cpp @@ -11,26 +11,26 @@ void Block::update_vessel_type(VesselType type) { vessel_type = type; } Block::~Block() {} -void Block::setup_params_(const std::vector ¶m_ids) { +void Block::setup_params_(const std::vector& param_ids) { this->global_param_ids = param_ids; } -void Block::setup_dofs_(DOFHandler &dofhandler, int num_equations, - const std::list &internal_var_names) { +void Block::setup_dofs_(DOFHandler& dofhandler, int num_equations, + const std::list& internal_var_names) { // Collect external DOFs from inlet nodes - for (auto &inlet_node : inlet_nodes) { + for (auto& inlet_node : inlet_nodes) { global_var_ids.push_back(inlet_node->pres_dof); global_var_ids.push_back(inlet_node->flow_dof); } // Collect external DOFs from outlet nodes - for (auto &outlet_node : outlet_nodes) { + for (auto& outlet_node : outlet_nodes) { global_var_ids.push_back(outlet_node->pres_dof); global_var_ids.push_back(outlet_node->flow_dof); } // Register internal variables of block - for (auto &int_name : internal_var_names) { + for (auto& int_name : internal_var_names) { global_var_ids.push_back( dofhandler.register_variable(int_name + ":" + this->get_name())); } @@ -41,30 +41,30 @@ void Block::setup_dofs_(DOFHandler &dofhandler, int num_equations, } } -void Block::setup_dofs(DOFHandler &dofhandler) {} +void Block::setup_dofs(DOFHandler& dofhandler) {} void Block::setup_model_dependent_params() {} void Block::setup_initial_state_dependent_params( - State initial_state, std::vector ¶meters) {} + State initial_state, std::vector& parameters) {} -void Block::update_constant(SparseSystem &system, - std::vector ¶meters) {} +void Block::update_constant(SparseSystem& system, + std::vector& parameters) {} -void Block::update_time(SparseSystem &system, std::vector ¶meters) { +void Block::update_time(SparseSystem& system, std::vector& parameters) { } void Block::update_solution( - SparseSystem &system, std::vector ¶meters, - const Eigen::Matrix &y, - const Eigen::Matrix &dy) {} + SparseSystem& system, std::vector& parameters, + const Eigen::Matrix& y, + const Eigen::Matrix& dy) {} -void Block::post_solve(Eigen::Matrix &y) {} +void Block::post_solve(Eigen::Matrix& y) {} -void Block::update_gradient(Eigen::SparseMatrix &jacobian, - Eigen::Matrix &residual, - Eigen::Matrix &alpha, - std::vector &y, std::vector &dy) { +void Block::update_gradient(Eigen::SparseMatrix& jacobian, + Eigen::Matrix& residual, + Eigen::Matrix& alpha, + std::vector& y, std::vector& dy) { throw std::runtime_error("Gradient calculation not implemented for block " + get_name()); } diff --git a/src/model/Block.h b/src/model/Block.h index e328bcea7..835ca837b 100644 --- a/src/model/Block.h +++ b/src/model/Block.h @@ -40,7 +40,7 @@ struct TripletsContributions { * number of contributions * @return The number of triplets */ - TripletsContributions operator+=(const TripletsContributions &other) { + TripletsContributions operator+=(const TripletsContributions& other) { F += other.F; E += other.E; D += other.D; @@ -75,15 +75,15 @@ class Model; class Block { public: const int id; ///< Global ID of the block - const Model *model; ///< The model to which the block belongs + const Model* model; ///< The model to which the block belongs const BlockType block_type; ///< Type of this block const BlockClass block_class; ///< Class of this block VesselType vessel_type = VesselType::neither; ///< Vessel type of this block const std::vector> input_params; ///< Map from name to input parameter - std::vector inlet_nodes; ///< Inlet nodes - std::vector outlet_nodes; ///< Outlet nodes + std::vector inlet_nodes; ///< Inlet nodes + std::vector outlet_nodes; ///< Outlet nodes bool steady = false; ///< Toggle steady behavior bool input_params_list = false; ///< Are input parameters given as a list? @@ -97,7 +97,7 @@ class Block { * @param block_class The class the block belongs to (e.g. vessel, junction) * @param input_params The parameters the block takes from the input file */ - Block(int id, Model *model, BlockType block_type, BlockClass block_class, + Block(int id, Model* model, BlockType block_type, BlockClass block_class, std::vector> input_params) : id(id), model(model), @@ -115,7 +115,7 @@ class Block { * @brief Copy the Block object * */ - Block(const Block &) = delete; + Block(const Block&) = delete; /** * @brief Global IDs for the block parameters. * @@ -166,7 +166,7 @@ class Block { * @brief Setup parameter IDs for the block * @param param_ids Global IDs of the block parameters */ - void setup_params_(const std::vector ¶m_ids); + void setup_params_(const std::vector& param_ids); /** * @brief Set up the degrees of freedom (DOF) of the block @@ -181,8 +181,8 @@ class Block { * @param internal_var_names Number of internal variables of the block */ - void setup_dofs_(DOFHandler &dofhandler, int num_equations, - const std::list &internal_var_names); + void setup_dofs_(DOFHandler& dofhandler, int num_equations, + const std::list& internal_var_names); /** * @brief Set up the degrees of freedom (DOF) of the block @@ -194,7 +194,7 @@ class Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - virtual void setup_dofs(DOFHandler &dofhandler); + virtual void setup_dofs(DOFHandler& dofhandler); /** * @brief Setup parameters that depend on the model @@ -209,7 +209,7 @@ class Block { * @param parameters The parameter values vector (at time 0) */ virtual void setup_initial_state_dependent_params( - State initial_state, std::vector ¶meters); + State initial_state, std::vector& parameters); /** * @brief Update the constant contributions of the element in a sparse system @@ -217,8 +217,8 @@ class Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - virtual void update_constant(SparseSystem &system, - std::vector ¶meters); + virtual void update_constant(SparseSystem& system, + std::vector& parameters); /** * @brief Update the time-dependent contributions of the element in a sparse * system @@ -226,8 +226,8 @@ class Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - virtual void update_time(SparseSystem &system, - std::vector ¶meters); + virtual void update_time(SparseSystem& system, + std::vector& parameters); /** * @brief Update the solution-dependent contributions of the element in a @@ -239,16 +239,16 @@ class Block { * @param dy Current derivate of the solution */ virtual void update_solution( - SparseSystem &system, std::vector ¶meters, - const Eigen::Matrix &y, - const Eigen::Matrix &dy); + SparseSystem& system, std::vector& parameters, + const Eigen::Matrix& y, + const Eigen::Matrix& dy); /** * @brief Modify the solution after solving it * * @param y Current solution */ - virtual void post_solve(Eigen::Matrix &y); + virtual void post_solve(Eigen::Matrix& y); /** * @brief Set the gradient of the block contributions with respect to the @@ -261,10 +261,10 @@ class Block { * @param dy Time-derivative of the current solution */ virtual void update_gradient( - Eigen::SparseMatrix &jacobian, - Eigen::Matrix &residual, - Eigen::Matrix &alpha, std::vector &y, - std::vector &dy); + Eigen::SparseMatrix& jacobian, + Eigen::Matrix& residual, + Eigen::Matrix& alpha, std::vector& y, + std::vector& dy); /** * @brief Number of triplets of element diff --git a/src/model/BloodVessel.cpp b/src/model/BloodVessel.cpp index aa290e0d3..faa6f25c1 100644 --- a/src/model/BloodVessel.cpp +++ b/src/model/BloodVessel.cpp @@ -3,12 +3,12 @@ #include "BloodVessel.h" -void BloodVessel::setup_dofs(DOFHandler &dofhandler) { +void BloodVessel::setup_dofs(DOFHandler& dofhandler) { Block::setup_dofs_(dofhandler, 2, {}); } -void BloodVessel::update_constant(SparseSystem &system, - std::vector ¶meters) { +void BloodVessel::update_constant(SparseSystem& system, + std::vector& parameters) { // Get parameters double capacitance = parameters[global_param_ids[ParamId::CAPACITANCE]]; double inductance = parameters[global_param_ids[ParamId::INDUCTANCE]]; @@ -31,9 +31,9 @@ void BloodVessel::update_constant(SparseSystem &system, } void BloodVessel::update_solution( - SparseSystem &system, std::vector ¶meters, - const Eigen::Matrix &y, - const Eigen::Matrix &dy) { + SparseSystem& system, std::vector& parameters, + const Eigen::Matrix& y, + const Eigen::Matrix& dy) { // Get parameters double capacitance = parameters[global_param_ids[ParamId::CAPACITANCE]]; double stenosis_coeff = @@ -57,10 +57,10 @@ void BloodVessel::update_solution( } void BloodVessel::update_gradient( - Eigen::SparseMatrix &jacobian, - Eigen::Matrix &residual, - Eigen::Matrix &alpha, std::vector &y, - std::vector &dy) { + Eigen::SparseMatrix& jacobian, + Eigen::Matrix& residual, + Eigen::Matrix& alpha, std::vector& y, + std::vector& dy) { auto y0 = y[global_var_ids[0]]; auto y1 = y[global_var_ids[1]]; auto y2 = y[global_var_ids[2]]; diff --git a/src/model/BloodVessel.h b/src/model/BloodVessel.h index 4182cc952..b2e029e43 100644 --- a/src/model/BloodVessel.h +++ b/src/model/BloodVessel.h @@ -136,7 +136,7 @@ class BloodVessel : public Block { * @param id Global ID of the block * @param model The model to which the block belongs */ - BloodVessel(int id, Model *model) + BloodVessel(int id, Model* model) : Block(id, model, BlockType::blood_vessel, BlockClass::vessel, {{"R_poiseuille", InputParameter()}, {"C", InputParameter(true)}, @@ -153,7 +153,7 @@ class BloodVessel : public Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); /** * @brief Update the constant contributions of the element in a sparse @@ -162,7 +162,7 @@ class BloodVessel : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_constant(SparseSystem &system, std::vector ¶meters); + void update_constant(SparseSystem& system, std::vector& parameters); /** * @brief Update the solution-dependent contributions of the element in a @@ -173,9 +173,9 @@ class BloodVessel : public Block { * @param y Current solution * @param dy Current derivate of the solution */ - void update_solution(SparseSystem &system, std::vector ¶meters, - const Eigen::Matrix &y, - const Eigen::Matrix &dy); + void update_solution(SparseSystem& system, std::vector& parameters, + const Eigen::Matrix& y, + const Eigen::Matrix& dy); /** * @brief Set the gradient of the block contributions with respect to the @@ -187,10 +187,10 @@ class BloodVessel : public Block { * @param y Current solution * @param dy Time-derivative of the current solution */ - void update_gradient(Eigen::SparseMatrix &jacobian, - Eigen::Matrix &residual, - Eigen::Matrix &alpha, - std::vector &y, std::vector &dy); + void update_gradient(Eigen::SparseMatrix& jacobian, + Eigen::Matrix& residual, + Eigen::Matrix& alpha, + std::vector& y, std::vector& dy); /** * @brief Number of triplets of element diff --git a/src/model/BloodVesselJunction.cpp b/src/model/BloodVesselJunction.cpp index 392cd9ae6..56469e31e 100644 --- a/src/model/BloodVesselJunction.cpp +++ b/src/model/BloodVesselJunction.cpp @@ -3,7 +3,7 @@ #include "BloodVesselJunction.h" -void BloodVesselJunction::setup_dofs(DOFHandler &dofhandler) { +void BloodVesselJunction::setup_dofs(DOFHandler& dofhandler) { if (inlet_nodes.size() != 1) { throw std::runtime_error( "Blood vessel junction does not support multiple inlets."); @@ -16,8 +16,8 @@ void BloodVesselJunction::setup_dofs(DOFHandler &dofhandler) { num_triplets.D = 2 * num_outlets; } -void BloodVesselJunction::update_constant(SparseSystem &system, - std::vector ¶meters) { +void BloodVesselJunction::update_constant(SparseSystem& system, + std::vector& parameters) { // Mass conservation system.F.coeffRef(global_eqn_ids[0], global_var_ids[1]) = 1.0; @@ -36,9 +36,9 @@ void BloodVesselJunction::update_constant(SparseSystem &system, } void BloodVesselJunction::update_solution( - SparseSystem &system, std::vector ¶meters, - const Eigen::Matrix &y, - const Eigen::Matrix &dy) { + SparseSystem& system, std::vector& parameters, + const Eigen::Matrix& y, + const Eigen::Matrix& dy) { for (size_t i = 0; i < num_outlets; i++) { // Get parameters auto stenosis_coeff = parameters[global_param_ids[2 * num_outlets + i]]; @@ -53,10 +53,10 @@ void BloodVesselJunction::update_solution( } void BloodVesselJunction::update_gradient( - Eigen::SparseMatrix &jacobian, - Eigen::Matrix &residual, - Eigen::Matrix &alpha, std::vector &y, - std::vector &dy) { + Eigen::SparseMatrix& jacobian, + Eigen::Matrix& residual, + Eigen::Matrix& alpha, std::vector& y, + std::vector& dy) { auto p_in = y[global_var_ids[0]]; auto q_in = y[global_var_ids[1]]; diff --git a/src/model/BloodVesselJunction.h b/src/model/BloodVesselJunction.h index 02ec7fefd..8e4bfa893 100644 --- a/src/model/BloodVesselJunction.h +++ b/src/model/BloodVesselJunction.h @@ -126,7 +126,7 @@ class BloodVesselJunction : public Block { * @param id Global ID of the block * @param model The model to which the block belongs */ - BloodVesselJunction(int id, Model *model) + BloodVesselJunction(int id, Model* model) : Block(id, model, BlockType::blood_vessel_junction, BlockClass::junction, {{"R_poiseuille", InputParameter()}, {"L", InputParameter()}, @@ -144,7 +144,7 @@ class BloodVesselJunction : public Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); /** * @brief Update the constant contributions of the element in a sparse system @@ -152,7 +152,7 @@ class BloodVesselJunction : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_constant(SparseSystem &system, std::vector ¶meters); + void update_constant(SparseSystem& system, std::vector& parameters); /** * @brief Update the solution-dependent contributions of the element in a @@ -164,9 +164,9 @@ class BloodVesselJunction : public Block { * @param dy Current derivate of the solution */ virtual void update_solution( - SparseSystem &system, std::vector ¶meters, - const Eigen::Matrix &y, - const Eigen::Matrix &dy); + SparseSystem& system, std::vector& parameters, + const Eigen::Matrix& y, + const Eigen::Matrix& dy); /** * @brief Set the gradient of the block contributions with respect to the @@ -178,10 +178,10 @@ class BloodVesselJunction : public Block { * @param y Current solution * @param dy Time-derivative of the current solution */ - void update_gradient(Eigen::SparseMatrix &jacobian, - Eigen::Matrix &residual, - Eigen::Matrix &alpha, - std::vector &y, std::vector &dy); + void update_gradient(Eigen::SparseMatrix& jacobian, + Eigen::Matrix& residual, + Eigen::Matrix& alpha, + std::vector& y, std::vector& dy); /** * @brief Number of triplets of element diff --git a/src/model/ChamberElastanceInductor.cpp b/src/model/ChamberElastanceInductor.cpp index fc1080347..440835e8a 100644 --- a/src/model/ChamberElastanceInductor.cpp +++ b/src/model/ChamberElastanceInductor.cpp @@ -2,13 +2,13 @@ // University of California, and others. SPDX-License-Identifier: BSD-3-Clause #include "ChamberElastanceInductor.h" -void ChamberElastanceInductor::setup_dofs(DOFHandler &dofhandler) { +void ChamberElastanceInductor::setup_dofs(DOFHandler& dofhandler) { // Internal variable is chamber volume Block::setup_dofs_(dofhandler, 3, {"Vc"}); } void ChamberElastanceInductor::update_constant( - SparseSystem &system, std::vector ¶meters) { + SparseSystem& system, std::vector& parameters) { double L = parameters[global_param_ids[ParamId::IMPEDANCE]]; // Eq 0: P_in - E(t)(Vc - Vrest) = 0 @@ -25,8 +25,8 @@ void ChamberElastanceInductor::update_constant( system.E.coeffRef(global_eqn_ids[2], global_var_ids[4]) = -1.0; } -void ChamberElastanceInductor::update_time(SparseSystem &system, - std::vector ¶meters) { +void ChamberElastanceInductor::update_time(SparseSystem& system, + std::vector& parameters) { get_elastance_values(parameters); // Eq 0: P_in - E(t)(Vc - Vrest) = P_in - E(t)*Vc + E(t)*Vrest = 0 @@ -35,7 +35,7 @@ void ChamberElastanceInductor::update_time(SparseSystem &system, } void ChamberElastanceInductor::get_elastance_values( - std::vector ¶meters) { + std::vector& parameters) { double Emax = parameters[global_param_ids[ParamId::EMAX]]; double Emin = parameters[global_param_ids[ParamId::EMIN]]; double Vrd = parameters[global_param_ids[ParamId::VRD]]; diff --git a/src/model/ChamberElastanceInductor.h b/src/model/ChamberElastanceInductor.h index aa5956603..7e398281c 100644 --- a/src/model/ChamberElastanceInductor.h +++ b/src/model/ChamberElastanceInductor.h @@ -124,7 +124,7 @@ class ChamberElastanceInductor : public Block { * @param id Global ID of the block * @param model The model to which the block belongs */ - ChamberElastanceInductor(int id, Model *model) + ChamberElastanceInductor(int id, Model* model) : Block(id, model, BlockType::chamber_elastance_inductor, BlockClass::chamber, {{"Emax", InputParameter()}, @@ -159,7 +159,7 @@ class ChamberElastanceInductor : public Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); /** * @brief Update the constant contributions of the element in a sparse @@ -168,7 +168,7 @@ class ChamberElastanceInductor : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_constant(SparseSystem &system, std::vector ¶meters); + void update_constant(SparseSystem& system, std::vector& parameters); /** * @brief Update the time-dependent contributions of the element in a sparse @@ -177,7 +177,7 @@ class ChamberElastanceInductor : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_time(SparseSystem &system, std::vector ¶meters); + void update_time(SparseSystem& system, std::vector& parameters); /** * @brief Number of triplets of element @@ -196,7 +196,7 @@ class ChamberElastanceInductor : public Block { * * @param parameters Parameters of the model */ - void get_elastance_values(std::vector ¶meters); + void get_elastance_values(std::vector& parameters); }; #endif // SVZERODSOLVER_MODEL_CHAMBERELASTANCEINDUCTOR_HPP_ diff --git a/src/model/ChamberSphere.cpp b/src/model/ChamberSphere.cpp index d6c25cfbc..db940432d 100644 --- a/src/model/ChamberSphere.cpp +++ b/src/model/ChamberSphere.cpp @@ -5,13 +5,13 @@ #include "Model.h" -void ChamberSphere::setup_dofs(DOFHandler &dofhandler) { +void ChamberSphere::setup_dofs(DOFHandler& dofhandler) { Block::setup_dofs_(dofhandler, 7, {"radius", "velo", "stress", "tau", "volume"}); } -void ChamberSphere::update_constant(SparseSystem &system, - std::vector ¶meters) { +void ChamberSphere::update_constant(SparseSystem& system, + std::vector& parameters) { const double thick0 = parameters[global_param_ids[ParamId::thick0]]; const double rho = parameters[global_param_ids[ParamId::rho]]; @@ -42,17 +42,17 @@ void ChamberSphere::update_constant(SparseSystem &system, system.F.coeffRef(global_eqn_ids[6], global_var_ids[2]) = -1; } -void ChamberSphere::update_time(SparseSystem &system, - std::vector ¶meters) { +void ChamberSphere::update_time(SparseSystem& system, + std::vector& parameters) { // active stress get_elastance_values(parameters); system.F.coeffRef(global_eqn_ids[3], global_var_ids[7]) = act; } void ChamberSphere::update_solution( - SparseSystem &system, std::vector ¶meters, - const Eigen::Matrix &y, - const Eigen::Matrix &dy) { + SparseSystem& system, std::vector& parameters, + const Eigen::Matrix& y, + const Eigen::Matrix& dy) { const double W1 = parameters[global_param_ids[ParamId::W1]]; const double W2 = parameters[global_param_ids[ParamId::W2]]; const double eta = parameters[global_param_ids[ParamId::eta]]; @@ -107,7 +107,7 @@ void ChamberSphere::update_solution( system.C.coeffRef(global_eqn_ids[3]) = -act_plus * sigma_max; } -void ChamberSphere::get_elastance_values(std::vector ¶meters) { +void ChamberSphere::get_elastance_values(std::vector& parameters) { const double alpha_max = parameters[global_param_ids[ParamId::alpha_max]]; const double alpha_min = parameters[global_param_ids[ParamId::alpha_min]]; const double tsys = parameters[global_param_ids[ParamId::tsys]]; diff --git a/src/model/ChamberSphere.h b/src/model/ChamberSphere.h index 7ef94f71d..734220d94 100644 --- a/src/model/ChamberSphere.h +++ b/src/model/ChamberSphere.h @@ -131,7 +131,7 @@ class ChamberSphere : public Block { * @param id Global ID of the block * @param model The model to which the block belongs */ - ChamberSphere(int id, Model *model) + ChamberSphere(int id, Model* model) : Block(id, model, BlockType::chamber_sphere, BlockClass::vessel, {{"rho", InputParameter()}, {"thick0", InputParameter()}, @@ -156,7 +156,7 @@ class ChamberSphere : public Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); /** * @brief Update the constant contributions of the element in a sparse @@ -165,7 +165,7 @@ class ChamberSphere : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_constant(SparseSystem &system, std::vector ¶meters); + void update_constant(SparseSystem& system, std::vector& parameters); /** * @brief Update the time-dependent contributions of the element in a sparse @@ -174,7 +174,7 @@ class ChamberSphere : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_time(SparseSystem &system, std::vector ¶meters); + void update_time(SparseSystem& system, std::vector& parameters); /** * @brief Update the solution-dependent contributions of the element in a @@ -185,16 +185,16 @@ class ChamberSphere : public Block { * @param y Current solution * @param dy Current derivate of the solution */ - void update_solution(SparseSystem &system, std::vector ¶meters, - const Eigen::Matrix &y, - const Eigen::Matrix &dy); + void update_solution(SparseSystem& system, std::vector& parameters, + const Eigen::Matrix& y, + const Eigen::Matrix& dy); /** * @brief Update the elastance functions which depend on time * * @param parameters Parameters of the model */ - void get_elastance_values(std::vector ¶meters); + void get_elastance_values(std::vector& parameters); private: double act = 0.0; // activation function diff --git a/src/model/ClosedLoopCoronaryBC.cpp b/src/model/ClosedLoopCoronaryBC.cpp index 9283dfa56..ef1f419ed 100644 --- a/src/model/ClosedLoopCoronaryBC.cpp +++ b/src/model/ClosedLoopCoronaryBC.cpp @@ -4,12 +4,12 @@ #include "Model.h" -void ClosedLoopCoronaryBC::setup_dofs(DOFHandler &dofhandler) { +void ClosedLoopCoronaryBC::setup_dofs(DOFHandler& dofhandler) { Block::setup_dofs_(dofhandler, 3, {"volume_im"}); } -void ClosedLoopCoronaryBC::update_constant(SparseSystem &system, - std::vector ¶meters) { +void ClosedLoopCoronaryBC::update_constant(SparseSystem& system, + std::vector& parameters) { auto ra = parameters[global_param_ids[ParamId::RA]]; auto ram = parameters[global_param_ids[ParamId::RAM]]; auto rv = parameters[global_param_ids[ParamId::RV]]; @@ -34,9 +34,9 @@ void ClosedLoopCoronaryBC::update_constant(SparseSystem &system, } void ClosedLoopCoronaryBC::update_solution( - SparseSystem &system, std::vector ¶meters, - const Eigen::Matrix &y, - const Eigen::Matrix &dy) { + SparseSystem& system, std::vector& parameters, + const Eigen::Matrix& y, + const Eigen::Matrix& dy) { auto cim = parameters[global_param_ids[ParamId::CIM]]; auto im = parameters[im_param_id]; auto pim = im * y[ventricle_var_id]; diff --git a/src/model/ClosedLoopCoronaryBC.h b/src/model/ClosedLoopCoronaryBC.h index b1c2cec14..41790042a 100644 --- a/src/model/ClosedLoopCoronaryBC.h +++ b/src/model/ClosedLoopCoronaryBC.h @@ -103,7 +103,7 @@ class ClosedLoopCoronaryBC : public Block { * @param model The model to which the block belongs * @param block_type The specific type of block (left or right) */ - ClosedLoopCoronaryBC(int id, Model *model, BlockType block_type) + ClosedLoopCoronaryBC(int id, Model* model, BlockType block_type) : Block(id, model, block_type, BlockClass::closed_loop, {{"Ra", InputParameter()}, {"Ram", InputParameter()}, @@ -127,7 +127,7 @@ class ClosedLoopCoronaryBC : public Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); /** * @brief Setup parameters that depend on the model @@ -141,7 +141,7 @@ class ClosedLoopCoronaryBC : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_constant(SparseSystem &system, std::vector ¶meters); + void update_constant(SparseSystem& system, std::vector& parameters); /** * @brief Update the solution-dependent contributions of the element in a @@ -152,9 +152,9 @@ class ClosedLoopCoronaryBC : public Block { * @param y Current solution * @param dy Current derivate of the solution */ - void update_solution(SparseSystem &system, std::vector ¶meters, - const Eigen::Matrix &y, - const Eigen::Matrix &dy); + void update_solution(SparseSystem& system, std::vector& parameters, + const Eigen::Matrix& y, + const Eigen::Matrix& dy); /** * @brief Number of triplets of element diff --git a/src/model/ClosedLoopCoronaryLeftBC.h b/src/model/ClosedLoopCoronaryLeftBC.h index 9fe083eb9..6830d0c4f 100644 --- a/src/model/ClosedLoopCoronaryLeftBC.h +++ b/src/model/ClosedLoopCoronaryLeftBC.h @@ -21,7 +21,7 @@ class ClosedLoopCoronaryLeftBC : public ClosedLoopCoronaryBC { * @param id Global ID of the block * @param model The model to which the block belongs */ - ClosedLoopCoronaryLeftBC(int id, Model *model) + ClosedLoopCoronaryLeftBC(int id, Model* model) : ClosedLoopCoronaryBC(id, model, BlockType::closed_loop_coronary_left_bc) {} diff --git a/src/model/ClosedLoopCoronaryRightBC.h b/src/model/ClosedLoopCoronaryRightBC.h index fff6c4908..f49457f73 100644 --- a/src/model/ClosedLoopCoronaryRightBC.h +++ b/src/model/ClosedLoopCoronaryRightBC.h @@ -21,7 +21,7 @@ class ClosedLoopCoronaryRightBC : public ClosedLoopCoronaryBC { * @param id Global ID of the block * @param model The model to which the block belongs */ - ClosedLoopCoronaryRightBC(int id, Model *model) + ClosedLoopCoronaryRightBC(int id, Model* model) : ClosedLoopCoronaryBC(id, model, BlockType::closed_loop_coronary_right_bc) {} diff --git a/src/model/ClosedLoopHeartPulmonary.cpp b/src/model/ClosedLoopHeartPulmonary.cpp index 2a317d567..fd684f632 100644 --- a/src/model/ClosedLoopHeartPulmonary.cpp +++ b/src/model/ClosedLoopHeartPulmonary.cpp @@ -4,14 +4,14 @@ #include "Model.h" -void ClosedLoopHeartPulmonary::setup_dofs(DOFHandler &dofhandler) { +void ClosedLoopHeartPulmonary::setup_dofs(DOFHandler& dofhandler) { Block::setup_dofs_(dofhandler, 14, {"V_RA", "Q_RA", "P_RV", "V_RV", "Q_RV", "P_pul", "P_LA", "V_LA", "Q_LA", "P_LV", "V_LV", "Q_LV"}); } void ClosedLoopHeartPulmonary::update_constant( - SparseSystem &system, std::vector ¶meters) { + SparseSystem& system, std::vector& parameters) { // DOF 0, Eq 0: Right atrium pressure system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = 1.0; @@ -75,8 +75,8 @@ void ClosedLoopHeartPulmonary::update_constant( parameters[global_param_ids[ParamId::LLV_A]]; } -void ClosedLoopHeartPulmonary::update_time(SparseSystem &system, - std::vector ¶meters) { +void ClosedLoopHeartPulmonary::update_time(SparseSystem& system, + std::vector& parameters) { get_activation_and_elastance_functions(parameters); // DOF 0, Eq 0: Right atrium pressure @@ -99,9 +99,9 @@ void ClosedLoopHeartPulmonary::update_time(SparseSystem &system, } void ClosedLoopHeartPulmonary::update_solution( - SparseSystem &system, std::vector ¶meters, - const Eigen::Matrix &y, - const Eigen::Matrix &dy) { + SparseSystem& system, std::vector& parameters, + const Eigen::Matrix& y, + const Eigen::Matrix& dy) { get_psi_ra_la(parameters, y); get_valve_positions(y); @@ -169,7 +169,7 @@ void ClosedLoopHeartPulmonary::update_solution( } void ClosedLoopHeartPulmonary::get_activation_and_elastance_functions( - std::vector ¶meters) { + std::vector& parameters) { auto T_cardiac = model->cardiac_cycle_period; auto Tsa = T_cardiac * parameters[global_param_ids[ParamId::TSA]]; auto tpwave = T_cardiac / parameters[global_param_ids[ParamId::TPWAVE]]; @@ -219,8 +219,8 @@ void ClosedLoopHeartPulmonary::get_activation_and_elastance_functions( } void ClosedLoopHeartPulmonary::get_psi_ra_la( - std::vector ¶meters, - const Eigen::Matrix &y) { + std::vector& parameters, + const Eigen::Matrix& y) { auto RA_volume = y[global_var_ids[4]]; auto LA_volume = y[global_var_ids[11]]; psi_ra = parameters[global_param_ids[ParamId::KXP_RA]] * @@ -245,7 +245,7 @@ void ClosedLoopHeartPulmonary::get_psi_ra_la( } void ClosedLoopHeartPulmonary::get_valve_positions( - const Eigen::Matrix &y) { + const Eigen::Matrix& y) { std::fill(valves, valves + 16, 1.0); // RA to RV @@ -280,7 +280,7 @@ void ClosedLoopHeartPulmonary::get_valve_positions( } void ClosedLoopHeartPulmonary::post_solve( - Eigen::Matrix &y) { + Eigen::Matrix& y) { for (size_t i = 0; i < 16; i++) if (valves[i] < 0.5) y[global_var_ids[i]] = 0.0; } diff --git a/src/model/ClosedLoopHeartPulmonary.h b/src/model/ClosedLoopHeartPulmonary.h index 81753cb5d..a775c6301 100644 --- a/src/model/ClosedLoopHeartPulmonary.h +++ b/src/model/ClosedLoopHeartPulmonary.h @@ -85,7 +85,7 @@ class ClosedLoopHeartPulmonary : public Block { * @param id Global ID of the block * @param model The model to which the block belongs */ - ClosedLoopHeartPulmonary(int id, Model *model) + ClosedLoopHeartPulmonary(int id, Model* model) : Block(id, model, BlockType::closed_loop_heart_pulmonary, BlockClass::closed_loop, {{"Tsa", InputParameter()}, {"tpwave", InputParameter()}, @@ -147,7 +147,7 @@ class ClosedLoopHeartPulmonary : public Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); /** * @brief Update the constant contributions of the element in a sparse @@ -156,7 +156,7 @@ class ClosedLoopHeartPulmonary : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_constant(SparseSystem &system, std::vector ¶meters); + void update_constant(SparseSystem& system, std::vector& parameters); /** * @brief Update the time-dependent contributions of the element in a sparse @@ -165,7 +165,7 @@ class ClosedLoopHeartPulmonary : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_time(SparseSystem &system, std::vector ¶meters); + void update_time(SparseSystem& system, std::vector& parameters); /** * @brief Update the solution-dependent contributions of the element in a @@ -176,16 +176,16 @@ class ClosedLoopHeartPulmonary : public Block { * @param y Current solution * @param dy Current derivate of the solution */ - void update_solution(SparseSystem &system, std::vector ¶meters, - const Eigen::Matrix &y, - const Eigen::Matrix &dy); + void update_solution(SparseSystem& system, std::vector& parameters, + const Eigen::Matrix& y, + const Eigen::Matrix& dy); /** * @brief Modify the solution after solving it * * @param y Current solution */ - void post_solve(Eigen::Matrix &y); + void post_solve(Eigen::Matrix& y); /** * @brief Number of triplets of element @@ -211,7 +211,7 @@ class ClosedLoopHeartPulmonary : public Block { * * @param parameters Parameters of the model */ - void get_activation_and_elastance_functions(std::vector ¶meters); + void get_activation_and_elastance_functions(std::vector& parameters); /** * @brief Compute sub-expressions that are part of atrial elastance and @@ -220,15 +220,15 @@ class ClosedLoopHeartPulmonary : public Block { * @param parameters Parameters of the model * @param y Current solution */ - void get_psi_ra_la(std::vector ¶meters, - const Eigen::Matrix &y); + void get_psi_ra_la(std::vector& parameters, + const Eigen::Matrix& y); /** * @brief Valve positions for each heart chamber * * @param y Current solution */ - void get_valve_positions(const Eigen::Matrix &y); + void get_valve_positions(const Eigen::Matrix& y); }; #endif // SVZERODSOLVER_MODEL_CLOSEDLOOPHEARTPULMONARY_HPP_ diff --git a/src/model/ClosedLoopRCRBC.cpp b/src/model/ClosedLoopRCRBC.cpp index c8f2d9677..b7fccdedd 100644 --- a/src/model/ClosedLoopRCRBC.cpp +++ b/src/model/ClosedLoopRCRBC.cpp @@ -2,12 +2,12 @@ // University of California, and others. SPDX-License-Identifier: BSD-3-Clause #include "ClosedLoopRCRBC.h" -void ClosedLoopRCRBC::setup_dofs(DOFHandler &dofhandler) { +void ClosedLoopRCRBC::setup_dofs(DOFHandler& dofhandler) { Block::setup_dofs_(dofhandler, 3, {"P_c"}); } -void ClosedLoopRCRBC::update_constant(SparseSystem &system, - std::vector ¶meters) { +void ClosedLoopRCRBC::update_constant(SparseSystem& system, + std::vector& parameters) { system.F.coeffRef(global_eqn_ids[0], global_var_ids[1]) = -1.0; system.F.coeffRef(global_eqn_ids[0], global_var_ids[3]) = 1.0; system.F.coeffRef(global_eqn_ids[1], global_var_ids[0]) = 1.0; diff --git a/src/model/ClosedLoopRCRBC.h b/src/model/ClosedLoopRCRBC.h index 6b895cc0b..5d46866f5 100644 --- a/src/model/ClosedLoopRCRBC.h +++ b/src/model/ClosedLoopRCRBC.h @@ -97,7 +97,7 @@ class ClosedLoopRCRBC : public Block { * @param id Global ID of the block * @param model The model to which the block belongs */ - ClosedLoopRCRBC(int id, Model *model) + ClosedLoopRCRBC(int id, Model* model) : Block(id, model, BlockType::closed_loop_rcr_bc, BlockClass::boundary_condition, {{"Rp", InputParameter()}, @@ -125,7 +125,7 @@ class ClosedLoopRCRBC : public Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); /** * @brief Update the constant contributions of the element in a sparse @@ -134,7 +134,7 @@ class ClosedLoopRCRBC : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_constant(SparseSystem &system, std::vector ¶meters); + void update_constant(SparseSystem& system, std::vector& parameters); /** * @brief Number of triplets of element diff --git a/src/model/FlowReferenceBC.cpp b/src/model/FlowReferenceBC.cpp index 9208c29a7..f57bf08e9 100644 --- a/src/model/FlowReferenceBC.cpp +++ b/src/model/FlowReferenceBC.cpp @@ -2,16 +2,16 @@ // University of California, and others. SPDX-License-Identifier: BSD-3-Clause #include "FlowReferenceBC.h" -void FlowReferenceBC::setup_dofs(DOFHandler &dofhandler) { +void FlowReferenceBC::setup_dofs(DOFHandler& dofhandler) { Block::setup_dofs_(dofhandler, 1, {}); } -void FlowReferenceBC::update_constant(SparseSystem &system, - std::vector ¶meters) { +void FlowReferenceBC::update_constant(SparseSystem& system, + std::vector& parameters) { system.F.coeffRef(global_eqn_ids[0], global_var_ids[1]) = 1.0; } -void FlowReferenceBC::update_time(SparseSystem &system, - std::vector ¶meters) { +void FlowReferenceBC::update_time(SparseSystem& system, + std::vector& parameters) { system.C(global_eqn_ids[0]) = -parameters[global_param_ids[0]]; } diff --git a/src/model/FlowReferenceBC.h b/src/model/FlowReferenceBC.h index 502f45bb7..e5c702228 100644 --- a/src/model/FlowReferenceBC.h +++ b/src/model/FlowReferenceBC.h @@ -63,7 +63,7 @@ class FlowReferenceBC : public Block { * @param id Global ID of the block * @param model The model to which the block belongs */ - FlowReferenceBC(int id, Model *model) + FlowReferenceBC(int id, Model* model) : Block(id, model, BlockType::flow_bc, BlockClass::boundary_condition, {{"t", InputParameter(false, true)}, {"Q", InputParameter(false, true)}}) {} @@ -78,7 +78,7 @@ class FlowReferenceBC : public Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); /** * @brief Update the constant contributions of the element in a sparse system @@ -86,7 +86,7 @@ class FlowReferenceBC : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_constant(SparseSystem &system, std::vector ¶meters); + void update_constant(SparseSystem& system, std::vector& parameters); /** * @brief Update the time-dependent contributions of the element in a sparse @@ -95,7 +95,7 @@ class FlowReferenceBC : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_time(SparseSystem &system, std::vector ¶meters); + void update_time(SparseSystem& system, std::vector& parameters); /** * @brief Number of triplets of element diff --git a/src/model/Junction.cpp b/src/model/Junction.cpp index 80a2a25ee..4f773df09 100644 --- a/src/model/Junction.cpp +++ b/src/model/Junction.cpp @@ -2,7 +2,7 @@ // University of California, and others. SPDX-License-Identifier: BSD-3-Clause #include "Junction.h" -void Junction::setup_dofs(DOFHandler &dofhandler) { +void Junction::setup_dofs(DOFHandler& dofhandler) { // Set number of equations of a junction block based on number of // inlets/outlets. Must be set before calling parent constructor num_inlets = inlet_nodes.size(); @@ -12,8 +12,8 @@ void Junction::setup_dofs(DOFHandler &dofhandler) { (num_inlets + num_outlets - 1) * 2 + num_inlets + num_outlets; } -void Junction::update_constant(SparseSystem &system, - std::vector ¶meters) { +void Junction::update_constant(SparseSystem& system, + std::vector& parameters) { // Pressure conservation for (size_t i = 0; i < (num_inlets + num_outlets - 1); i++) { system.F.coeffRef(global_eqn_ids[i], global_var_ids[0]) = 1.0; @@ -33,10 +33,10 @@ void Junction::update_constant(SparseSystem &system, } void Junction::update_gradient( - Eigen::SparseMatrix &jacobian, - Eigen::Matrix &residual, - Eigen::Matrix &alpha, std::vector &y, - std::vector &dy) { + Eigen::SparseMatrix& jacobian, + Eigen::Matrix& residual, + Eigen::Matrix& alpha, std::vector& y, + std::vector& dy) { // Pressure conservation residual(global_eqn_ids[0]) = y[global_var_ids[0]] - y[global_var_ids[2]]; diff --git a/src/model/Junction.h b/src/model/Junction.h index d365616a1..d8b7b0deb 100644 --- a/src/model/Junction.h +++ b/src/model/Junction.h @@ -81,7 +81,7 @@ class Junction : public Block { * @param id Global ID of the block * @param model The model to which the block belongs */ - Junction(int id, Model *model) + Junction(int id, Model* model) : Block(id, model, BlockType::junction, BlockClass::junction, {}) {} /** * @brief Set up the degrees of freedom (DOF) of the block @@ -93,7 +93,7 @@ class Junction : public Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); /** * @brief Update the constant contributions of the element in a sparse system @@ -101,7 +101,7 @@ class Junction : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_constant(SparseSystem &system, std::vector ¶meters); + void update_constant(SparseSystem& system, std::vector& parameters); /** * @brief Set the gradient of the block contributions with respect to the @@ -113,10 +113,10 @@ class Junction : public Block { * @param y Current solution * @param dy Time-derivative of the current solution */ - void update_gradient(Eigen::SparseMatrix &jacobian, - Eigen::Matrix &residual, - Eigen::Matrix &alpha, - std::vector &y, std::vector &dy); + void update_gradient(Eigen::SparseMatrix& jacobian, + Eigen::Matrix& residual, + Eigen::Matrix& alpha, + std::vector& y, std::vector& dy); /** * @brief Number of triplets of element diff --git a/src/model/Model.cpp b/src/model/Model.cpp index eae74f6d0..340c54d1f 100644 --- a/src/model/Model.cpp +++ b/src/model/Model.cpp @@ -4,7 +4,7 @@ template BlockFactoryFunc block_factory() { - return [](int count, Model *model) -> Block * { + return [](int count, Model* model) -> Block* { return new block_type(count, model); }; } @@ -33,18 +33,18 @@ Model::Model() { Model::~Model() {} -Block *Model::create_block(const std::string &block_type) { +Block* Model::create_block(const std::string& block_type) { // Get block from factory auto it = block_factory_map.find(block_type); if (it == block_factory_map.end()) { throw std::runtime_error("Invalid block type " + block_type); } - Block *block = it->second(block_count, this); + Block* block = it->second(block_count, this); return block; } -int Model::add_block(Block *block, const std::string_view &name, - const std::vector &block_param_ids, bool internal) { +int Model::add_block(Block* block, const std::string_view& name, + const std::vector& block_param_ids, bool internal) { // Set global parameter IDs block->setup_params_(block_param_ids); @@ -63,9 +63,9 @@ int Model::add_block(Block *block, const std::string_view &name, return block_count++; } -int Model::add_block(const std::string &block_name, - const std::vector &block_param_ids, - const std::string_view &name, bool internal) { +int Model::add_block(const std::string& block_name, + const std::vector& block_param_ids, + const std::string_view& name, bool internal) { // Generate block from factory auto block = this->create_block(block_name); @@ -73,7 +73,7 @@ int Model::add_block(const std::string &block_name, return this->add_block(block, name, block_param_ids, internal); } -bool Model::has_block(const std::string &name) const { +bool Model::has_block(const std::string& name) const { if (block_index_map.find(name) == block_index_map.end()) { return false; } else { @@ -81,7 +81,7 @@ bool Model::has_block(const std::string &name) const { } } -Block *Model::get_block(const std::string_view &name) const { +Block* Model::get_block(const std::string_view& name) const { auto name_string = static_cast(name); if (!has_block(name_string)) { @@ -91,7 +91,7 @@ Block *Model::get_block(const std::string_view &name) const { return blocks[block_index_map.at(name_string)].get(); } -Block *Model::get_block(int block_id) const { +Block* Model::get_block(int block_id) const { if (block_id >= blocks.size()) { return hidden_blocks[block_id - blocks.size()].get(); } @@ -99,7 +99,7 @@ Block *Model::get_block(int block_id) const { return blocks[block_id].get(); } -BlockType Model::get_block_type(const std::string_view &name) const { +BlockType Model::get_block_type(const std::string_view& name) const { auto name_string = static_cast(name); if (block_index_map.find(name_string) == block_index_map.end()) { @@ -113,9 +113,9 @@ std::string Model::get_block_name(int block_id) const { return block_names[block_id]; } -int Model::add_node(const std::vector &inlet_eles, - const std::vector &outlet_eles, - const std::string_view &name) { +int Model::add_node(const std::vector& inlet_eles, + const std::vector& outlet_eles, + const std::string_view& name) { // DEBUG_MSG("Adding node " << name); auto node = std::shared_ptr( new Node(node_count, inlet_eles, outlet_eles, this)); @@ -135,8 +135,8 @@ int Model::add_parameter(double value) { return parameter_count++; } -int Model::add_parameter(const std::vector ×, - const std::vector &values, bool periodic) { +int Model::add_parameter(const std::vector& times, + const std::vector& values, bool periodic) { auto param = Parameter(parameter_count, times, values, periodic); if (periodic && (param.is_constant == false)) { if ((this->cardiac_cycle_period > 0.0) && @@ -151,7 +151,7 @@ int Model::add_parameter(const std::vector ×, return parameter_count++; } -Parameter *Model::get_parameter(int param_id) { return ¶meters[param_id]; } +Parameter* Model::get_parameter(int param_id) { return ¶meters[param_id]; } double Model::get_parameter_value(int param_id) const { return parameter_values[param_id]; @@ -163,15 +163,15 @@ void Model::update_parameter_value(int param_id, double param_value) { void Model::finalize() { DEBUG_MSG("Setup degrees-of-freedom of nodes"); - for (auto &node : nodes) { + for (auto& node : nodes) { node->setup_dofs(dofhandler); } DEBUG_MSG("Setup degrees-of-freedom of blocks"); - for (auto &block : blocks) { + for (auto& block : blocks) { block->setup_dofs(dofhandler); } DEBUG_MSG("Setup model-dependent parameters"); - for (auto &block : blocks) { + for (auto& block : blocks) { block->setup_model_dependent_params(); } } @@ -186,16 +186,16 @@ int Model::get_num_blocks(bool internal) const { return num_blocks; } -void Model::update_constant(SparseSystem &system) { +void Model::update_constant(SparseSystem& system) { for (auto block : blocks) { block->update_constant(system, parameter_values); } } -void Model::update_time(SparseSystem &system, double time) { +void Model::update_time(SparseSystem& system, double time) { this->time = time; - for (auto ¶m : parameters) { + for (auto& param : parameters) { parameter_values[param.id] = param.get(time); } @@ -204,22 +204,22 @@ void Model::update_time(SparseSystem &system, double time) { } } -void Model::update_solution(SparseSystem &system, - Eigen::Matrix &y, - Eigen::Matrix &dy) { +void Model::update_solution(SparseSystem& system, + Eigen::Matrix& y, + Eigen::Matrix& dy) { for (auto block : blocks) { block->update_solution(system, parameter_values, y, dy); } } -void Model::post_solve(Eigen::Matrix &y) { +void Model::post_solve(Eigen::Matrix& y) { for (auto block : blocks) { block->post_solve(y); } } void Model::to_steady() { - for (auto ¶m : parameters) { + for (auto& param : parameters) { param.to_steady(); } @@ -237,10 +237,10 @@ void Model::to_steady() { } void Model::to_unsteady() { - for (auto ¶m : parameters) { + for (auto& param : parameters) { param.to_unsteady(); } - for (auto &[param_id_capacitance, value] : param_value_cache) { + for (auto& [param_id_capacitance, value] : param_value_cache) { // DEBUG_MSG("Setting Windkessel capacitance back to " << value); parameters[param_id_capacitance].update(value); } @@ -252,7 +252,7 @@ void Model::to_unsteady() { TripletsContributions Model::get_num_triplets() const { TripletsContributions triplets_sum; - for (auto &elem : blocks) { + for (auto& elem : blocks) { triplets_sum += elem->get_num_triplets(); } @@ -261,7 +261,7 @@ TripletsContributions Model::get_num_triplets() const { void Model::setup_initial_state_dependent_parameters(State initial_state) { DEBUG_MSG("Setup initial state dependent parameters"); - for (auto &block : blocks) { + for (auto& block : blocks) { block->setup_initial_state_dependent_params(initial_state, parameter_values); } diff --git a/src/model/Model.h b/src/model/Model.h index 5ab0662c5..2a2dcc52c 100644 --- a/src/model/Model.h +++ b/src/model/Model.h @@ -74,7 +74,7 @@ class Model { * @param block_name The block name (defined in block_factory_map) * @return int Global ID of the block */ - Block *create_block(const std::string &block_name); + Block* create_block(const std::string& block_name); /** * @brief Add a block to the model (without parameters) @@ -85,8 +85,8 @@ class Model { * @param internal Toggle whether block is internal * @return int Global ID of the block */ - int add_block(Block *block, const std::string_view &name, - const std::vector &block_param_ids, bool internal = false); + int add_block(Block* block, const std::string_view& name, + const std::vector& block_param_ids, bool internal = false); /** * @brief Add a block to the model (with parameters) @@ -97,9 +97,9 @@ class Model { * @param internal Toggle whether block is internal * @return int Global ID of the block */ - int add_block(const std::string &block_name, - const std::vector &block_param_ids, - const std::string_view &name, bool internal = false); + int add_block(const std::string& block_name, + const std::vector& block_param_ids, + const std::string_view& name, bool internal = false); /** * @brief Check if a block with given name exists @@ -107,7 +107,7 @@ class Model { * @param name Name of the Block * @return bool whether block exists */ - bool has_block(const std::string &name) const; + bool has_block(const std::string& name) const; /** * @brief Get a block by its name @@ -115,7 +115,7 @@ class Model { * @param name Name of the Block * @return Block* The block */ - Block *get_block(const std::string_view &name) const; + Block* get_block(const std::string_view& name) const; /** * @brief Get a block by its global ID @@ -123,7 +123,7 @@ class Model { * @param block_id Global ID of the Block * @return Block* The block */ - Block *get_block(int block_id) const; + Block* get_block(int block_id) const; /** * @brief Get a block type by its name @@ -131,7 +131,7 @@ class Model { * @param name The name of the block * @return BlockType The block type */ - BlockType get_block_type(const std::string_view &name) const; + BlockType get_block_type(const std::string_view& name) const; /** * @brief Get the name of a block by it's ID @@ -149,9 +149,9 @@ class Model { * @param name Name of node * @return int Global ID of the node */ - int add_node(const std::vector &inlet_eles, - const std::vector &outlet_eles, - const std::string_view &name); + int add_node(const std::vector& inlet_eles, + const std::vector& outlet_eles, + const std::string_view& name); /** * @brief Get the name of a node by it's ID @@ -177,8 +177,8 @@ class Model { * @param periodic Toggle whether parameter is periodic * @return int Global ID of the parameter */ - int add_parameter(const std::vector ×, - const std::vector &values, bool periodic = true); + int add_parameter(const std::vector& times, + const std::vector& values, bool periodic = true); /** * @brief Get a parameter by its global ID @@ -186,7 +186,7 @@ class Model { * @param param_id Global ID of the parameter * @return Parameter* The parameter */ - Parameter *get_parameter(int param_id); + Parameter* get_parameter(int param_id); /** * @brief Get the current value of a parameter @@ -218,7 +218,7 @@ class Model { * * @param system System to update contributions at */ - void update_constant(SparseSystem &system); + void update_constant(SparseSystem& system); /** * @brief Update the time-dependent contributions of all elements in a sparse @@ -227,7 +227,7 @@ class Model { * @param system System to update contributions at * @param time Current time */ - void update_time(SparseSystem &system, double time); + void update_time(SparseSystem& system, double time); /** * @brief Update the solution-dependent contributions of all elements in a @@ -237,16 +237,16 @@ class Model { * @param y Current solution * @param dy Current derivate of the solution */ - void update_solution(SparseSystem &system, - Eigen::Matrix &y, - Eigen::Matrix &dy); + void update_solution(SparseSystem& system, + Eigen::Matrix& y, + Eigen::Matrix& dy); /** * @brief Modify the solution after solving it * * @param y Current solution */ - void post_solve(Eigen::Matrix &y); + void post_solve(Eigen::Matrix& y); /** * @brief Convert the blocks to a steady behavior diff --git a/src/model/Node.cpp b/src/model/Node.cpp index f0bd67f06..5eddcbf8e 100644 --- a/src/model/Node.cpp +++ b/src/model/Node.cpp @@ -5,25 +5,25 @@ #include "Block.h" #include "Model.h" -Node::Node(int id, const std::vector &inlet_eles, - const std::vector &outlet_eles, Model *model) { +Node::Node(int id, const std::vector& inlet_eles, + const std::vector& outlet_eles, Model* model) { this->id = id; this->inlet_eles = inlet_eles; this->outlet_eles = outlet_eles; this->model = model; - for (auto &inlet_ele : inlet_eles) { + for (auto& inlet_ele : inlet_eles) { inlet_ele->outlet_nodes.push_back(this); } - for (auto &outlet_ele : outlet_eles) { + for (auto& outlet_ele : outlet_eles) { outlet_ele->inlet_nodes.push_back(this); } } std::string Node::get_name() { return this->model->get_node_name(this->id); } -void Node::setup_dofs(DOFHandler &dofhandler) { +void Node::setup_dofs(DOFHandler& dofhandler) { flow_dof = dofhandler.register_variable("flow:" + get_name()); pres_dof = dofhandler.register_variable("pressure:" + get_name()); } diff --git a/src/model/Node.h b/src/model/Node.h index fd0e43343..3a1cad65c 100644 --- a/src/model/Node.h +++ b/src/model/Node.h @@ -33,13 +33,13 @@ class Node { * @param outlet_eles Outlet element of the node * @param model The model to which the node belongs */ - Node(int id, const std::vector &inlet_eles, - const std::vector &outlet_eles, Model *model); + Node(int id, const std::vector& inlet_eles, + const std::vector& outlet_eles, Model* model); - int id; ///< Global ID of the block - std::vector inlet_eles; ///< Inlet element of the node - std::vector outlet_eles; ///< Outlet element of the node - Model *model{nullptr}; ///< The model to which the node belongs + int id; ///< Global ID of the block + std::vector inlet_eles; ///< Inlet element of the node + std::vector outlet_eles; ///< Outlet element of the node + Model* model{nullptr}; ///< The model to which the node belongs int flow_dof{0}; ///< Global flow degree-of-freedom of the node int pres_dof{0}; ///< Global pressure degree-of-freedom of the node @@ -61,7 +61,7 @@ class Node { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); }; #endif // SVZERODSOLVER_MODEL_NODE_HPP_ diff --git a/src/model/OpenLoopCoronaryBC.cpp b/src/model/OpenLoopCoronaryBC.cpp index 7cb295765..dc3f924ba 100644 --- a/src/model/OpenLoopCoronaryBC.cpp +++ b/src/model/OpenLoopCoronaryBC.cpp @@ -2,12 +2,12 @@ // University of California, and others. SPDX-License-Identifier: BSD-3-Clause #include "OpenLoopCoronaryBC.h" -void OpenLoopCoronaryBC::setup_dofs(DOFHandler &dofhandler) { +void OpenLoopCoronaryBC::setup_dofs(DOFHandler& dofhandler) { Block::setup_dofs_(dofhandler, 2, {"volume_im"}); } -void OpenLoopCoronaryBC::update_constant(SparseSystem &system, - std::vector ¶meters) { +void OpenLoopCoronaryBC::update_constant(SparseSystem& system, + std::vector& parameters) { auto Ra = parameters[global_param_ids[0]]; auto Ram = parameters[global_param_ids[1]]; auto Rv = parameters[global_param_ids[2]]; @@ -37,8 +37,8 @@ void OpenLoopCoronaryBC::update_constant(SparseSystem &system, } } -void OpenLoopCoronaryBC::update_time(SparseSystem &system, - std::vector ¶meters) { +void OpenLoopCoronaryBC::update_time(SparseSystem& system, + std::vector& parameters) { auto Ram = parameters[global_param_ids[1]]; auto Rv = parameters[global_param_ids[2]]; auto Cim = parameters[global_param_ids[4]]; @@ -57,7 +57,7 @@ void OpenLoopCoronaryBC::update_time(SparseSystem &system, } void OpenLoopCoronaryBC::setup_initial_state_dependent_params( - State initial_state, std::vector ¶meters) { + State initial_state, std::vector& parameters) { auto P_in = initial_state.y[global_var_ids[0]]; auto Q_in = initial_state.y[global_var_ids[1]]; auto P_in_dot = initial_state.ydot[global_var_ids[0]]; diff --git a/src/model/OpenLoopCoronaryBC.h b/src/model/OpenLoopCoronaryBC.h index bbbc4ccca..ba710dca3 100644 --- a/src/model/OpenLoopCoronaryBC.h +++ b/src/model/OpenLoopCoronaryBC.h @@ -95,7 +95,7 @@ class OpenLoopCoronaryBC : public Block { * @param id Global ID of the block * @param model The model to which the block belongs */ - OpenLoopCoronaryBC(int id, Model *model) + OpenLoopCoronaryBC(int id, Model* model) : Block(id, model, BlockType::open_loop_coronary_bc, BlockClass::boundary_condition, {{"Ra1", InputParameter()}, @@ -118,7 +118,7 @@ class OpenLoopCoronaryBC : public Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); /** * @brief Setup parameters that depend on the initial state @@ -127,7 +127,7 @@ class OpenLoopCoronaryBC : public Block { * @param parameters The parameter values vector (at time 0) */ void setup_initial_state_dependent_params(State initial_state, - std::vector ¶meters); + std::vector& parameters); /** * @brief Update the constant contributions of the element in a sparse system @@ -135,7 +135,7 @@ class OpenLoopCoronaryBC : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_constant(SparseSystem &system, std::vector ¶meters); + void update_constant(SparseSystem& system, std::vector& parameters); /** * @brief Update the time-dependent contributions of the element in a sparse @@ -144,7 +144,7 @@ class OpenLoopCoronaryBC : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_time(SparseSystem &system, std::vector ¶meters); + void update_time(SparseSystem& system, std::vector& parameters); /** * @brief Number of triplets of element diff --git a/src/model/Parameter.cpp b/src/model/Parameter.cpp index f7f28c2b8..3a43154fd 100644 --- a/src/model/Parameter.cpp +++ b/src/model/Parameter.cpp @@ -7,8 +7,8 @@ Parameter::Parameter(int id, double value) { update(value); } -Parameter::Parameter(int id, const std::vector ×, - const std::vector &values, bool periodic) { +Parameter::Parameter(int id, const std::vector& times, + const std::vector& values, bool periodic) { this->id = id; this->is_periodic = periodic; update(times, values); @@ -20,8 +20,8 @@ void Parameter::update(double update_value) { value = update_value; } -void Parameter::update(const std::vector &update_times, - const std::vector &update_values) { +void Parameter::update(const std::vector& update_times, + const std::vector& update_values) { this->size = update_values.size(); if (size == 1) { diff --git a/src/model/PressureReferenceBC.cpp b/src/model/PressureReferenceBC.cpp index 8f28e62a7..2d589be35 100644 --- a/src/model/PressureReferenceBC.cpp +++ b/src/model/PressureReferenceBC.cpp @@ -2,16 +2,16 @@ // University of California, and others. SPDX-License-Identifier: BSD-3-Clause #include "PressureReferenceBC.h" -void PressureReferenceBC::setup_dofs(DOFHandler &dofhandler) { +void PressureReferenceBC::setup_dofs(DOFHandler& dofhandler) { Block::setup_dofs_(dofhandler, 1, {}); } -void PressureReferenceBC::update_constant(SparseSystem &system, - std::vector ¶meters) { +void PressureReferenceBC::update_constant(SparseSystem& system, + std::vector& parameters) { system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = 1.0; } -void PressureReferenceBC::update_time(SparseSystem &system, - std::vector ¶meters) { +void PressureReferenceBC::update_time(SparseSystem& system, + std::vector& parameters) { system.C(global_eqn_ids[0]) = -parameters[global_param_ids[0]]; } diff --git a/src/model/PressureReferenceBC.h b/src/model/PressureReferenceBC.h index 8bd38d2d9..19bbe5373 100644 --- a/src/model/PressureReferenceBC.h +++ b/src/model/PressureReferenceBC.h @@ -64,7 +64,7 @@ class PressureReferenceBC : public Block { * @param id Global ID of the block * @param model The model to which the block belongs */ - PressureReferenceBC(int id, Model *model) + PressureReferenceBC(int id, Model* model) : Block(id, model, BlockType::pressure_bc, BlockClass::boundary_condition, {{"t", InputParameter(false, true)}, {"P", InputParameter(false, true)}}) {} @@ -79,7 +79,7 @@ class PressureReferenceBC : public Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); /** * @brief Update the constant contributions of the element in a sparse system @@ -87,7 +87,7 @@ class PressureReferenceBC : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_constant(SparseSystem &system, std::vector ¶meters); + void update_constant(SparseSystem& system, std::vector& parameters); /** * @brief Update the time-dependent contributions of the element in a sparse @@ -96,7 +96,7 @@ class PressureReferenceBC : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_time(SparseSystem &system, std::vector ¶meters); + void update_time(SparseSystem& system, std::vector& parameters); /** * @brief Number of triplets of element diff --git a/src/model/ResistanceBC.cpp b/src/model/ResistanceBC.cpp index 5d706e97b..44c0b99b8 100644 --- a/src/model/ResistanceBC.cpp +++ b/src/model/ResistanceBC.cpp @@ -2,17 +2,17 @@ // University of California, and others. SPDX-License-Identifier: BSD-3-Clause #include "ResistanceBC.h" -void ResistanceBC::setup_dofs(DOFHandler &dofhandler) { +void ResistanceBC::setup_dofs(DOFHandler& dofhandler) { Block::setup_dofs_(dofhandler, 1, {}); } -void ResistanceBC::update_constant(SparseSystem &system, - std::vector ¶meters) { +void ResistanceBC::update_constant(SparseSystem& system, + std::vector& parameters) { system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = 1.0; } -void ResistanceBC::update_time(SparseSystem &system, - std::vector ¶meters) { +void ResistanceBC::update_time(SparseSystem& system, + std::vector& parameters) { system.F.coeffRef(global_eqn_ids[0], global_var_ids[1]) = -parameters[global_param_ids[0]]; system.C(global_eqn_ids[0]) = -parameters[global_param_ids[1]]; diff --git a/src/model/ResistanceBC.h b/src/model/ResistanceBC.h index ea320d210..f8b56d6cc 100644 --- a/src/model/ResistanceBC.h +++ b/src/model/ResistanceBC.h @@ -63,7 +63,7 @@ class ResistanceBC : public Block { * @param id Global ID of the block * @param model The model to which the block belongs */ - ResistanceBC(int id, Model *model) + ResistanceBC(int id, Model* model) : Block(id, model, BlockType::resistance_bc, BlockClass::boundary_condition, {{"R", InputParameter()}, {"Pd", InputParameter()}}) {} @@ -78,7 +78,7 @@ class ResistanceBC : public Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); /** * @brief Update the constant contributions of the element in a sparse system @@ -86,7 +86,7 @@ class ResistanceBC : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_constant(SparseSystem &system, std::vector ¶meters); + void update_constant(SparseSystem& system, std::vector& parameters); /** * @brief Update the time-dependent contributions of the element in a sparse @@ -95,7 +95,7 @@ class ResistanceBC : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_time(SparseSystem &system, std::vector ¶meters); + void update_time(SparseSystem& system, std::vector& parameters); /** * @brief Number of triplets of element diff --git a/src/model/ResistiveJunction.cpp b/src/model/ResistiveJunction.cpp index 80504d320..9cab54eac 100644 --- a/src/model/ResistiveJunction.cpp +++ b/src/model/ResistiveJunction.cpp @@ -2,7 +2,7 @@ // University of California, and others. SPDX-License-Identifier: BSD-3-Clause #include "ResistiveJunction.h" -void ResistiveJunction::setup_dofs(DOFHandler &dofhandler) { +void ResistiveJunction::setup_dofs(DOFHandler& dofhandler) { // Set number of equations of a junction block based on number of // inlets/outlets. Must be set before calling parent constructor num_inlets = inlet_nodes.size(); @@ -11,8 +11,8 @@ void ResistiveJunction::setup_dofs(DOFHandler &dofhandler) { num_triplets.F = (num_inlets + num_outlets) * 4; } -void ResistiveJunction::update_constant(SparseSystem &system, - std::vector ¶meters) { +void ResistiveJunction::update_constant(SparseSystem& system, + std::vector& parameters) { for (size_t i = 0; i < num_inlets; i++) { system.F.coeffRef(global_eqn_ids[i], global_var_ids[i * 2]) = 1.0; system.F.coeffRef(global_eqn_ids[i], global_var_ids[i * 2 + 1]) = diff --git a/src/model/ResistiveJunction.h b/src/model/ResistiveJunction.h index 10f7b21a1..bebda9f7c 100644 --- a/src/model/ResistiveJunction.h +++ b/src/model/ResistiveJunction.h @@ -89,7 +89,7 @@ class ResistiveJunction : public Block { * @param id Global ID of the block * @param model The model to which the block belongs */ - ResistiveJunction(int id, Model *model) + ResistiveJunction(int id, Model* model) : Block(id, model, BlockType::resistive_junction, BlockClass::junction, {{"R", InputParameter()}}) {} @@ -103,7 +103,7 @@ class ResistiveJunction : public Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); /** * @brief Update the constant contributions of the element in a sparse system @@ -111,7 +111,7 @@ class ResistiveJunction : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_constant(SparseSystem &system, std::vector ¶meters); + void update_constant(SparseSystem& system, std::vector& parameters); /** * @brief Number of triplets of element diff --git a/src/model/ValveTanh.cpp b/src/model/ValveTanh.cpp index a18c316cd..aaef5029d 100644 --- a/src/model/ValveTanh.cpp +++ b/src/model/ValveTanh.cpp @@ -2,7 +2,7 @@ // University of California, and others. SPDX-License-Identifier: BSD-3-Clause #include "ValveTanh.h" -void ValveTanh::setup_dofs(DOFHandler &dofhandler) { +void ValveTanh::setup_dofs(DOFHandler& dofhandler) { // set_up_dofs args: dofhandler (passed in), num equations, list of internal // variable names (strings) 2 eqns, one for Pressure, one for Flow Block::setup_dofs_(dofhandler, 2, {}); @@ -10,8 +10,8 @@ void ValveTanh::setup_dofs(DOFHandler &dofhandler) { // update_constant updates matrices E and F from E(y,t)*y_dot + F(y,t)*y + // c(y,t) = 0 with terms that DO NOT DEPEND ON THE SOLUTION -void ValveTanh::update_constant(SparseSystem &system, - std::vector ¶meters) { +void ValveTanh::update_constant(SparseSystem& system, + std::vector& parameters) { // Set element contributions // coeffRef args are the indices (i,j) of the matrix // global_eqn_ids: number of rows in the matrix, set in setup_dofs @@ -31,9 +31,9 @@ void ValveTanh::update_constant(SparseSystem &system, // c(y,t) = 0 with terms that DO DEPEND ON THE SOLUTION (will change with each // time step) void ValveTanh::update_solution( - SparseSystem &system, std::vector ¶meters, - const Eigen::Matrix &y, - const Eigen::Matrix &dy) { + SparseSystem& system, std::vector& parameters, + const Eigen::Matrix& y, + const Eigen::Matrix& dy) { // Get states double p_in = y[global_var_ids[0]]; double p_out = y[global_var_ids[2]]; diff --git a/src/model/ValveTanh.h b/src/model/ValveTanh.h index 8e173bcd9..92a963920 100644 --- a/src/model/ValveTanh.h +++ b/src/model/ValveTanh.h @@ -125,7 +125,7 @@ class ValveTanh : public Block { * @param id Global ID of the block * @param model The model to which the block belongs */ - ValveTanh(int id, Model *model) + ValveTanh(int id, Model* model) : Block(id, model, BlockType::valve_tanh, BlockClass::valve, {{"Rmax", InputParameter()}, {"Rmin", InputParameter()}, @@ -143,7 +143,7 @@ class ValveTanh : public Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); /** * @brief Update the constant contributions of the element in a sparse @@ -152,7 +152,7 @@ class ValveTanh : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_constant(SparseSystem &system, std::vector ¶meters); + void update_constant(SparseSystem& system, std::vector& parameters); /** * @brief Update the solution-dependent contributions of the element in a @@ -163,9 +163,9 @@ class ValveTanh : public Block { * @param y Current solution * @param dy Current derivate of the solution */ - void update_solution(SparseSystem &system, std::vector ¶meters, - const Eigen::Matrix &y, - const Eigen::Matrix &dy); + void update_solution(SparseSystem& system, std::vector& parameters, + const Eigen::Matrix& y, + const Eigen::Matrix& dy); /** * @brief Number of triplets of element diff --git a/src/model/WindkesselBC.cpp b/src/model/WindkesselBC.cpp index 7a2432cb3..90a0b9950 100644 --- a/src/model/WindkesselBC.cpp +++ b/src/model/WindkesselBC.cpp @@ -2,21 +2,21 @@ // University of California, and others. SPDX-License-Identifier: BSD-3-Clause #include "WindkesselBC.h" -void WindkesselBC::setup_dofs(DOFHandler &dofhandler) { +void WindkesselBC::setup_dofs(DOFHandler& dofhandler) { Block::setup_dofs_(dofhandler, 2, {"pressure_c"}); } -void WindkesselBC::update_constant(SparseSystem &system, +void WindkesselBC::update_constant(SparseSystem& system, - std::vector ¶meters) { + std::vector& parameters) { system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = 1.0; system.F.coeffRef(global_eqn_ids[0], global_var_ids[2]) = -1.0; system.F.coeffRef(global_eqn_ids[1], global_var_ids[2]) = -1.0; } -void WindkesselBC::update_time(SparseSystem &system, +void WindkesselBC::update_time(SparseSystem& system, - std::vector ¶meters) { + std::vector& parameters) { system.E.coeffRef(global_eqn_ids[1], global_var_ids[2]) = -parameters[global_param_ids[2]] * parameters[global_param_ids[1]]; system.F.coeffRef(global_eqn_ids[0], global_var_ids[1]) = diff --git a/src/model/WindkesselBC.h b/src/model/WindkesselBC.h index cde8e3e4b..435450a56 100644 --- a/src/model/WindkesselBC.h +++ b/src/model/WindkesselBC.h @@ -91,7 +91,7 @@ class WindkesselBC : public Block { * @param id Global ID of the block * @param model The model to which the block belongs */ - WindkesselBC(int id, Model *model) + WindkesselBC(int id, Model* model) : Block(id, model, BlockType::windkessel_bc, BlockClass::boundary_condition, {{"Rp", InputParameter()}, @@ -109,7 +109,7 @@ class WindkesselBC : public Block { * @param dofhandler Degree-of-freedom handler to register variables and * equations at */ - void setup_dofs(DOFHandler &dofhandler); + void setup_dofs(DOFHandler& dofhandler); /** * @brief Update the constant contributions of the element in a sparse @@ -118,7 +118,7 @@ class WindkesselBC : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_constant(SparseSystem &system, std::vector ¶meters); + void update_constant(SparseSystem& system, std::vector& parameters); /** * @brief Update the time-dependent contributions of the element in a sparse @@ -127,7 +127,7 @@ class WindkesselBC : public Block { * @param system System to update contributions at * @param parameters Parameters of the model */ - void update_time(SparseSystem &system, std::vector ¶meters); + void update_time(SparseSystem& system, std::vector& parameters); /** * @brief Number of triplets of element diff --git a/src/optimize/calibrate.cpp b/src/optimize/calibrate.cpp index 7471e5469..f5ff60604 100644 --- a/src/optimize/calibrate.cpp +++ b/src/optimize/calibrate.cpp @@ -5,12 +5,12 @@ #include "LevenbergMarquardtOptimizer.h" #include "SimulationParameters.h" -nlohmann::json calibrate(const nlohmann::json &config) { +nlohmann::json calibrate(const nlohmann::json& config) { auto output_config = nlohmann::json(config); // Read calibration parameters DEBUG_MSG("Parse calibration parameters"); - auto const &calibration_parameters = config["calibration_parameters"]; + auto const& calibration_parameters = config["calibration_parameters"]; double gradient_tol = calibration_parameters.value("tolerance_gradient", 1e-5); double increment_tol = @@ -37,7 +37,7 @@ nlohmann::json calibrate(const nlohmann::json &config) { DEBUG_MSG("Load vessels"); std::map vessel_id_map; int param_counter = 0; - for (auto const &vessel_config : config["vessels"]) { + for (auto const& vessel_config : config["vessels"]) { std::string vessel_name = vessel_config["vessel_name"]; // Create parameter IDs @@ -50,7 +50,7 @@ nlohmann::json calibrate(const nlohmann::json &config) { // Read connected boundary conditions if (vessel_config.contains("boundary_conditions")) { - auto const &vessel_bc_config = vessel_config["boundary_conditions"]; + auto const& vessel_bc_config = vessel_config["boundary_conditions"]; if (vessel_bc_config.contains("inlet")) { inlet_connections.push_back({vessel_bc_config["inlet"], vessel_name}); } @@ -61,9 +61,9 @@ nlohmann::json calibrate(const nlohmann::json &config) { } // Create junctions - for (auto const &junction_config : config["junctions"]) { + for (auto const& junction_config : config["junctions"]) { std::string junction_name = junction_config["junction_name"]; - auto const &outlet_vessels = junction_config["outlet_vessels"]; + auto const& outlet_vessels = junction_config["outlet_vessels"]; int num_outlets = outlet_vessels.size(); if (num_outlets == 1) { @@ -90,16 +90,16 @@ nlohmann::json calibrate(const nlohmann::json &config) { // Create Connections DEBUG_MSG("Created connection"); - for (auto &connection : connections) { + for (auto& connection : connections) { auto ele1 = model.get_block(std::get<0>(connection)); auto ele2 = model.get_block(std::get<1>(connection)); model.add_node({ele1}, {ele2}, ele1->get_name() + ":" + ele2->get_name()); } - for (auto &connection : inlet_connections) { + for (auto& connection : inlet_connections) { auto ele = model.get_block(std::get<1>(connection)); model.add_node({}, {ele}, std::get<0>(connection) + ":" + ele->get_name()); } - for (auto &connection : outlet_connections) { + for (auto& connection : outlet_connections) { auto ele = model.get_block(std::get<0>(connection)); model.add_node({ele}, {}, ele->get_name() + ":" + std::get<1>(connection)); } @@ -147,7 +147,7 @@ nlohmann::json calibrate(const nlohmann::json &config) { Eigen::Matrix alpha = Eigen::Matrix::Zero(param_counter); DEBUG_MSG("Reading initial alpha"); - for (auto &vessel_config : output_config["vessels"]) { + for (auto& vessel_config : output_config["vessels"]) { std::string vessel_name = vessel_config["vessel_name"]; DEBUG_MSG("Reading initial alpha for " << vessel_name); auto block = model.get_block(vessel_name); @@ -163,7 +163,7 @@ nlohmann::json calibrate(const nlohmann::json &config) { 0.0); } } - for (auto &junction_config : output_config["junctions"]) { + for (auto& junction_config : output_config["junctions"]) { std::string junction_name = junction_config["junction_name"]; DEBUG_MSG("Reading initial alpha for " << junction_name); auto block = model.get_block(junction_name); @@ -208,7 +208,7 @@ nlohmann::json calibrate(const nlohmann::json &config) { alpha = lm_alg.run(alpha, y_all, dy_all); // Write optimized simulation config file - for (auto &vessel_config : output_config["vessels"]) { + for (auto& vessel_config : output_config["vessels"]) { std::string vessel_name = vessel_config["vessel_name"]; auto block = model.get_block(vessel_name); double stenosis_coeff = 0.0; @@ -225,7 +225,7 @@ nlohmann::json calibrate(const nlohmann::json &config) { {"L", std::max(alpha[block->global_param_ids[2]], 0.0)}, {"stenosis_coefficient", stenosis_coeff}}; } - for (auto &junction_config : output_config["junctions"]) { + for (auto& junction_config : output_config["junctions"]) { std::string junction_name = junction_config["junction_name"]; auto block = model.get_block(junction_name); int num_outlets = block->outlet_nodes.size(); diff --git a/src/optimize/calibrate.h b/src/optimize/calibrate.h index 128388efe..a8b24b857 100644 --- a/src/optimize/calibrate.h +++ b/src/optimize/calibrate.h @@ -21,6 +21,6 @@ * @param config JSON configuration for 0D model * @return Calibrated JSON configuration for the 0D model */ -nlohmann::json calibrate(const nlohmann::json &config); +nlohmann::json calibrate(const nlohmann::json& config); #endif // SVZERODSOLVER_OPTIMIZE_CALIBRATOR_HPP_ diff --git a/src/solve/csv_writer.cpp b/src/solve/csv_writer.cpp index deca81057..fa64bb16d 100644 --- a/src/solve/csv_writer.cpp +++ b/src/solve/csv_writer.cpp @@ -15,8 +15,8 @@ * @param derivative Toggle whether to output time-derivatives * @return CSV encoded output string */ -std::string to_vessel_csv(const std::vector ×, - const std::vector &states, const Model &model, +std::string to_vessel_csv(const std::vector& times, + const std::vector& states, const Model& model, bool mean, bool derivative) { // Create string stream to buffer output std::stringstream out; @@ -43,8 +43,8 @@ std::string to_vessel_csv(const std::vector ×, auto block = model.get_block(i); // Extract global solution indices of the block - if (dynamic_cast(block) == nullptr && - dynamic_cast(block) == nullptr) { + if (dynamic_cast(block) == nullptr && + dynamic_cast(block) == nullptr) { continue; } @@ -145,9 +145,9 @@ std::string to_vessel_csv(const std::vector ×, * @param derivative Toggle whether to output time-derivatives * @return CSV encoded output string */ -std::string to_variable_csv(const std::vector ×, - const std::vector &states, - const Model &model, bool mean, bool derivative) { +std::string to_variable_csv(const std::vector& times, + const std::vector& states, + const Model& model, bool mean, bool derivative) { // Create string stream to buffer output std::stringstream out; diff --git a/src/solve/csv_writer.h b/src/solve/csv_writer.h index 0d31ffdae..c4cec3e59 100644 --- a/src/solve/csv_writer.h +++ b/src/solve/csv_writer.h @@ -14,13 +14,13 @@ #include "Model.h" #include "State.h" -std::string to_variable_csv(const std::vector ×, - const std::vector &states, - const Model &model, bool mean = false, +std::string to_variable_csv(const std::vector& times, + const std::vector& states, + const Model& model, bool mean = false, bool derivative = false); -std::string to_vessel_csv(const std::vector ×, - const std::vector &states, const Model &model, +std::string to_vessel_csv(const std::vector& times, + const std::vector& states, const Model& model, bool mean = false, bool derivative = false); #endif // SVZERODSOLVER_IO_CSVWRITER_HPP_ diff --git a/tests/test_interface/LPNSolverInterface/LPNSolverInterface.cpp b/tests/test_interface/LPNSolverInterface/LPNSolverInterface.cpp index 829e818a6..a69294204 100644 --- a/tests/test_interface/LPNSolverInterface/LPNSolverInterface.cpp +++ b/tests/test_interface/LPNSolverInterface/LPNSolverInterface.cpp @@ -1,4 +1,5 @@ #include "LPNSolverInterface.h" + #include #include @@ -6,9 +7,8 @@ // LPNSolverInterface //-------------------- // -LPNSolverInterface::LPNSolverInterface() -{ -//// Set the default names of the LPN interface functions. +LPNSolverInterface::LPNSolverInterface() { + //// Set the default names of the LPN interface functions. lpn_initialize_name_ = "initialize"; lpn_increment_time_name_ = "increment_time"; lpn_run_simulation_name_ = "run_simulation"; @@ -21,29 +21,27 @@ LPNSolverInterface::LPNSolverInterface() lpn_set_external_step_size_name_ = "set_external_step_size"; } -LPNSolverInterface::~LPNSolverInterface() -{ - dlclose(library_handle_); -} +LPNSolverInterface::~LPNSolverInterface() { dlclose(library_handle_); } //-------------- // load_library //-------------- // Load the LPN shared library and get pointers to its interface functions. // -void LPNSolverInterface::load_library(const std::string& interface_lib) -{ +void LPNSolverInterface::load_library(const std::string& interface_lib) { library_handle_ = dlopen(interface_lib.c_str(), RTLD_LAZY); if (!library_handle_) { - std::cerr << "Error loading shared library '" << interface_lib << "' with error: " << dlerror() << std::endl; + std::cerr << "Error loading shared library '" << interface_lib + << "' with error: " << dlerror() << std::endl; return; } // Get a pointer to the svzero 'initialize' function. *(void**)(&lpn_initialize_) = dlsym(library_handle_, "initialize"); if (!lpn_initialize_) { - std::cerr << "Error loading function 'initialize' with error: " << dlerror() << std::endl; + std::cerr << "Error loading function 'initialize' with error: " << dlerror() + << std::endl; dlclose(library_handle_); return; } @@ -51,71 +49,85 @@ void LPNSolverInterface::load_library(const std::string& interface_lib) // Get a pointer to the svzero 'increment_time' function. *(void**)(&lpn_increment_time_) = dlsym(library_handle_, "increment_time"); if (!lpn_increment_time_) { - std::cerr << "Error loading function 'increment_time' with error: " << dlerror() << std::endl; + std::cerr << "Error loading function 'increment_time' with error: " + << dlerror() << std::endl; dlclose(library_handle_); return; } - + // Get a pointer to the svzero 'run_simulation' function. *(void**)(&lpn_run_simulation_) = dlsym(library_handle_, "run_simulation"); if (!lpn_run_simulation_) { - std::cerr << "Error loading function 'run_simulation' with error: " << dlerror() << std::endl; + std::cerr << "Error loading function 'run_simulation' with error: " + << dlerror() << std::endl; dlclose(library_handle_); return; } // Get a pointer to the svzero 'update_block_params' function. - *(void**)(&lpn_update_block_params_) = dlsym(library_handle_, "update_block_params"); + *(void**)(&lpn_update_block_params_) = + dlsym(library_handle_, "update_block_params"); if (!lpn_update_block_params_) { - std::cerr << "Error loading function 'update_block_params' with error: " << dlerror() << std::endl; + std::cerr << "Error loading function 'update_block_params' with error: " + << dlerror() << std::endl; dlclose(library_handle_); return; } // Get a pointer to the svzero 'read_block_params' function. - *(void**)(&lpn_read_block_params_) = dlsym(library_handle_, "read_block_params"); + *(void**)(&lpn_read_block_params_) = + dlsym(library_handle_, "read_block_params"); if (!lpn_read_block_params_) { - std::cerr << "Error loading function 'read_block_params' with error: " << dlerror() << std::endl; + std::cerr << "Error loading function 'read_block_params' with error: " + << dlerror() << std::endl; dlclose(library_handle_); return; } - + // Get a pointer to the svzero 'get_block_node_IDs' function. - *(void**)(&lpn_get_block_node_IDs_) = dlsym(library_handle_, "get_block_node_IDs"); + *(void**)(&lpn_get_block_node_IDs_) = + dlsym(library_handle_, "get_block_node_IDs"); if (!lpn_get_block_node_IDs_) { - std::cerr << "Error loading function 'lpn_get_block_node_IDs' with error: " << dlerror() << std::endl; + std::cerr << "Error loading function 'lpn_get_block_node_IDs' with error: " + << dlerror() << std::endl; dlclose(library_handle_); return; } - + // Get a pointer to the svzero 'update_state' function. *(void**)(&lpn_update_state_) = dlsym(library_handle_, "update_state"); if (!lpn_update_state_) { - std::cerr << "Error loading function 'lpn_update_state' with error: " << dlerror() << std::endl; + std::cerr << "Error loading function 'lpn_update_state' with error: " + << dlerror() << std::endl; dlclose(library_handle_); return; } - + // Get a pointer to the svzero 'return_y' function. *(void**)(&lpn_return_y_) = dlsym(library_handle_, "return_y"); if (!lpn_return_y_) { - std::cerr << "Error loading function 'lpn_return_y' with error: " << dlerror() << std::endl; + std::cerr << "Error loading function 'lpn_return_y' with error: " + << dlerror() << std::endl; dlclose(library_handle_); return; } - + // Get a pointer to the svzero 'return_ydot' function. *(void**)(&lpn_return_ydot_) = dlsym(library_handle_, "return_ydot"); if (!lpn_return_ydot_) { - std::cerr << "Error loading function 'lpn_return_ydot' with error: " << dlerror() << std::endl; + std::cerr << "Error loading function 'lpn_return_ydot' with error: " + << dlerror() << std::endl; dlclose(library_handle_); return; } - + // Get a pointer to the svzero 'set_external_step_size' function. - *(void**)(&lpn_set_external_step_size_) = dlsym(library_handle_, "set_external_step_size"); + *(void**)(&lpn_set_external_step_size_) = + dlsym(library_handle_, "set_external_step_size"); if (!lpn_set_external_step_size_) { - std::cerr << "Error loading function 'lpn_set_external_step_size' with error: " << dlerror() << std::endl; + std::cerr + << "Error loading function 'lpn_set_external_step_size' with error: " + << dlerror() << std::endl; dlclose(library_handle_); return; } @@ -127,22 +139,23 @@ void LPNSolverInterface::load_library(const std::string& interface_lib) // // file_name: The name of the LPN configuration file (JSON). // -void LPNSolverInterface::initialize(std::string file_name) -{ - lpn_initialize_(file_name, problem_id_, pts_per_cycle_, num_cycles_, num_output_steps_, block_names_, variable_names_); - std::cout << "[LPNSolverInterface::initialize] Problem ID: " << problem_id_ << std::endl; +void LPNSolverInterface::initialize(std::string file_name) { + lpn_initialize_(file_name, problem_id_, pts_per_cycle_, num_cycles_, + num_output_steps_, block_names_, variable_names_); + std::cout << "[LPNSolverInterface::initialize] Problem ID: " << problem_id_ + << std::endl; system_size_ = variable_names_.size(); - std::cout << "[LPNSolverInterface::initialize] System size: " << system_size_ << std::endl; + std::cout << "[LPNSolverInterface::initialize] System size: " << system_size_ + << std::endl; } // Set the external time step variable in the svZeroD interface. // // Parameters: // -// step_size: The time step in the 3D (external) solver. +// step_size: The time step in the 3D (external) solver. // -void LPNSolverInterface::set_external_step_size(double step_size) -{ +void LPNSolverInterface::set_external_step_size(double step_size) { lpn_set_external_step_size_(problem_id_, step_size); } @@ -154,8 +167,8 @@ void LPNSolverInterface::set_external_step_size(double step_size) // // solution: The returned LPN solution. // -void LPNSolverInterface::increment_time(const double time, std::vector& solution) -{ +void LPNSolverInterface::increment_time(const double time, + std::vector& solution) { lpn_increment_time_(problem_id_, time, solution); } @@ -171,9 +184,12 @@ void LPNSolverInterface::increment_time(const double time, std::vector& // // error_code: Either 0 or 1 depending on whether the 0D simulation diverged. // -void LPNSolverInterface::run_simulation(const double time, std::vector& output_times, std::vector& output_solutions, int& error_code) -{ - lpn_run_simulation_(problem_id_, time, output_times, output_solutions, error_code); +void LPNSolverInterface::run_simulation(const double time, + std::vector& output_times, + std::vector& output_solutions, + int& error_code) { + lpn_run_simulation_(problem_id_, time, output_times, output_solutions, + error_code); } // Update the parameters of a particular 0D block @@ -184,8 +200,8 @@ void LPNSolverInterface::run_simulation(const double time, std::vector& // // new_params: The new parameters for the 0D block. // -void LPNSolverInterface::update_block_params(std::string block_name, std::vector& new_params) -{ +void LPNSolverInterface::update_block_params(std::string block_name, + std::vector& new_params) { lpn_update_block_params_(problem_id_, block_name, new_params); } @@ -197,12 +213,13 @@ void LPNSolverInterface::update_block_params(std::string block_name, std::vector // // new_params: The parameters for the 0D block. // -void LPNSolverInterface::read_block_params(std::string block_name, std::vector& block_params) -{ +void LPNSolverInterface::read_block_params(std::string block_name, + std::vector& block_params) { lpn_read_block_params_(problem_id_, block_name, block_params); } -// Get the IDs of the inlet/outlet variables of a given block in the state vector +// Get the IDs of the inlet/outlet variables of a given block in the state +// vector // // Parameters: // @@ -210,8 +227,8 @@ void LPNSolverInterface::read_block_params(std::string block_name, std::vector& IDs) -{ +void LPNSolverInterface::get_block_node_IDs(std::string block_name, + std::vector& IDs) { lpn_get_block_node_IDs_(problem_id_, block_name, IDs); } @@ -222,8 +239,8 @@ void LPNSolverInterface::get_block_node_IDs(std::string block_name, std::vector< // state_y: The y state vector // // state_ydot: The ydot state vector -void LPNSolverInterface::update_state(std::vector state_y, std::vector state_ydot) -{ +void LPNSolverInterface::update_state(std::vector state_y, + std::vector state_ydot) { lpn_update_state_(problem_id_, state_y, state_ydot); } @@ -233,8 +250,7 @@ void LPNSolverInterface::update_state(std::vector state_y, std::vector& y) -{ +void LPNSolverInterface::return_y(std::vector& y) { lpn_return_y_(problem_id_, y); } @@ -244,7 +260,6 @@ void LPNSolverInterface::return_y(std::vector& y) // // ydot: The y state vector // -void LPNSolverInterface::return_ydot(std::vector& ydot) -{ +void LPNSolverInterface::return_ydot(std::vector& ydot) { lpn_return_ydot_(problem_id_, ydot); } diff --git a/tests/test_interface/LPNSolverInterface/LPNSolverInterface.h b/tests/test_interface/LPNSolverInterface/LPNSolverInterface.h index a3a568204..deb9fa7da 100644 --- a/tests/test_interface/LPNSolverInterface/LPNSolverInterface.h +++ b/tests/test_interface/LPNSolverInterface/LPNSolverInterface.h @@ -3,48 +3,41 @@ /* ---------- Windows implementation ---------- */ #include #ifdef interface - #undef interface +#undef interface #endif using dl_handle_t = HMODULE; // Define windows flags to emulate #ifndef RTLD_LAZY - #define RTLD_LAZY 0 - #define RTLD_NOW 0 - #define RTLD_GLOBAL 0 - #define RTLD_LOCAL 0 +#define RTLD_LAZY 0 +#define RTLD_NOW 0 +#define RTLD_GLOBAL 0 +#define RTLD_LOCAL 0 #endif -inline dl_handle_t dlopen(const char* file, int /*flags*/) -{ +inline dl_handle_t dlopen(const char* file, int /*flags*/) { /* LoadLibraryA allows UTF-8 compatible narrow strings under MSVC ≥ 2015 */ return ::LoadLibraryA(file); } -inline void* dlsym(dl_handle_t handle, const char* symbol) -{ +inline void* dlsym(dl_handle_t handle, const char* symbol) { return reinterpret_cast(::GetProcAddress(handle, symbol)); } -inline int dlclose(dl_handle_t handle) -{ +inline int dlclose(dl_handle_t handle) { return ::FreeLibrary(handle) ? 0 : 1; // 0 = success, POSIX-style } /* Store the last error message in a local static buffer and return a C-string * (roughly mimicking the POSIX API). */ -inline const char* dlerror() -{ - static char buf[256]; - DWORD code = ::GetLastError(); +inline const char* dlerror() { + static char buf[256]; + DWORD code = ::GetLastError(); if (code == 0) return nullptr; - ::FormatMessageA(FORMAT_MESSAGE_FROM_SYSTEM | - FORMAT_MESSAGE_IGNORE_INSERTS, - nullptr, - code, - MAKELANGID(LANG_NEUTRAL, SUBLANG_DEFAULT), + ::FormatMessageA(FORMAT_MESSAGE_FROM_SYSTEM | FORMAT_MESSAGE_IGNORE_INSERTS, + nullptr, code, MAKELANGID(LANG_NEUTRAL, SUBLANG_DEFAULT), buf, sizeof(buf), nullptr); return buf; } @@ -56,10 +49,11 @@ using dl_handle_t = void*; #include #include + +#include #include #include -#include - + #ifndef LPNSolverInterface_h #define LPNSolverInterface_h @@ -67,64 +61,75 @@ using dl_handle_t = void*; // LPNSolverInterface //-------------------- // -class LPNSolverInterface -{ - public: - LPNSolverInterface(); - ~LPNSolverInterface(); - - void load_library(const std::string& interface_lib); - void initialize(std::string file_name); - void increment_time(const double time, std::vector& solution); - void run_simulation(const double time, std::vector& output_times, std::vector& output_solutions, int& error_code); - void update_block_params(std::string block_name, std::vector& new_params); - void read_block_params(std::string block_name, std::vector& block_params); - void get_block_node_IDs(std::string block_name, std::vector& IDs); - void update_state(std::vector state_y, std::vector state_ydot); - void return_y(std::vector& y); - void return_ydot(std::vector& ydot); - void set_external_step_size(double step_size); - - // Interface functions. - std::string lpn_initialize_name_; - void (*lpn_initialize_)(std::string, int&, int&, int&, int&, std::vector&, std::vector&); - - std::string lpn_increment_time_name_; - void (*lpn_increment_time_)(const int, const double, std::vector& solution); - - std::string lpn_run_simulation_name_; - void (*lpn_run_simulation_)(const int, const double, std::vector& output_times, std::vector& output_solutions, int& error_code); - - std::string lpn_update_block_params_name_; - void (*lpn_update_block_params_)(const int, std::string, std::vector& new_params); - - std::string lpn_read_block_params_name_; - void (*lpn_read_block_params_)(const int, std::string, std::vector& block_params); - - std::string lpn_get_block_node_IDs_name_; - void (*lpn_get_block_node_IDs_)(const int, std::string, std::vector& block_params); - - std::string lpn_update_state_name_; - void (*lpn_update_state_)(const int, std::vector, std::vector); - - std::string lpn_return_y_name_; - void (*lpn_return_y_)(const int, std::vector&); - - std::string lpn_return_ydot_name_; - void (*lpn_return_ydot_)(const int, std::vector&); - - std::string lpn_set_external_step_size_name_; - void (*lpn_set_external_step_size_)(const int, double); - - dl_handle_t library_handle_ = nullptr; - int problem_id_ = 0; - int system_size_ = 0; - int num_cycles_ = 0; - int pts_per_cycle_ = 0; - int num_output_steps_ = 0; - std::vector block_names_; - std::vector variable_names_; +class LPNSolverInterface { + public: + LPNSolverInterface(); + ~LPNSolverInterface(); + + void load_library(const std::string& interface_lib); + void initialize(std::string file_name); + void increment_time(const double time, std::vector& solution); + void run_simulation(const double time, std::vector& output_times, + std::vector& output_solutions, int& error_code); + void update_block_params(std::string block_name, + std::vector& new_params); + void read_block_params(std::string block_name, + std::vector& block_params); + void get_block_node_IDs(std::string block_name, std::vector& IDs); + void update_state(std::vector state_y, + std::vector state_ydot); + void return_y(std::vector& y); + void return_ydot(std::vector& ydot); + void set_external_step_size(double step_size); + + // Interface functions. + std::string lpn_initialize_name_; + void (*lpn_initialize_)(std::string, int&, int&, int&, int&, + std::vector&, std::vector&); + + std::string lpn_increment_time_name_; + void (*lpn_increment_time_)(const int, const double, + std::vector& solution); + + std::string lpn_run_simulation_name_; + void (*lpn_run_simulation_)(const int, const double, + std::vector& output_times, + std::vector& output_solutions, + int& error_code); + + std::string lpn_update_block_params_name_; + void (*lpn_update_block_params_)(const int, std::string, + std::vector& new_params); + + std::string lpn_read_block_params_name_; + void (*lpn_read_block_params_)(const int, std::string, + std::vector& block_params); + + std::string lpn_get_block_node_IDs_name_; + void (*lpn_get_block_node_IDs_)(const int, std::string, + std::vector& block_params); + + std::string lpn_update_state_name_; + void (*lpn_update_state_)(const int, std::vector, + std::vector); + + std::string lpn_return_y_name_; + void (*lpn_return_y_)(const int, std::vector&); + + std::string lpn_return_ydot_name_; + void (*lpn_return_ydot_)(const int, std::vector&); + + std::string lpn_set_external_step_size_name_; + void (*lpn_set_external_step_size_)(const int, double); + + dl_handle_t library_handle_ = nullptr; + int problem_id_ = 0; + int system_size_ = 0; + int num_cycles_ = 0; + int pts_per_cycle_ = 0; + int num_output_steps_ = 0; + std::vector block_names_; + std::vector variable_names_; }; #endif - diff --git a/tests/test_interface/test_01/main.cpp b/tests/test_interface/test_01/main.cpp index 3df83a1c3..0a1e1f8d3 100644 --- a/tests/test_interface/test_01/main.cpp +++ b/tests/test_interface/test_01/main.cpp @@ -1,31 +1,33 @@ -// Test interfacing to svZeroSolver. -// This test mimics an external 3D solver (svSolver/svFSI) interfacing with svZeroDSolver -// The model consists of a closed-loop heart model with coronary BCs -// It is run for one time step of the external solver +// Test interfacing to svZeroSolver. +// This test mimics an external 3D solver (svSolver/svFSI) interfacing with +// svZeroDSolver The model consists of a closed-loop heart model with coronary +// BCs It is run for one time step of the external solver -#include "../LPNSolverInterface/LPNSolverInterface.h" +#include +#include #include #include -#include -#include + +#include "../LPNSolverInterface/LPNSolverInterface.h" namespace fs = std::filesystem; //------ // main //------ // -int main(int argc, char** argv) -{ +int main(int argc, char** argv) { LPNSolverInterface interface; if (argc != 3) { - std::runtime_error("Usage: svZeroD_interface_test01 "); + std::runtime_error( + "Usage: svZeroD_interface_test01 " + ""); } - + // Load shared library and get interface functions. // File extension of the shared library depends on the system - fs::path build_dir = argv[1]; - fs::path iface_dir = build_dir / "src" / "interface"; + 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"; @@ -36,13 +38,15 @@ int main(int argc, char** argv) } 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() + " !"); + throw std::runtime_error("Could not find shared libraries " + + lib_so.string() + " or " + lib_dylib.string() + + " or " + lib_dll.string() + " !"); } // Set up the svZeroD model std::string file_name = std::string(argv[2]); interface.initialize(file_name); - + // Check number of variables and blocks if (interface.system_size_ != 133) { throw std::runtime_error("interface.system_size_ != 133"); @@ -54,38 +58,89 @@ int main(int argc, char** argv) // Set external time step size (flow solver step size) double external_step_size = 0.000923; interface.set_external_step_size(external_step_size); - + // Save the initial condition - std::vector init_state_y = {220.655, 113.454, 0.146379, 107.558, 0.0840239, 108.31, 0.0917877, 108.966, 0.0539358, 109.893, 0.0997981, 107.152, 0.168397, 109.693, 0.0478851, 108.683, 0.0969178, 106.587, 0.0745793, 111.186, 0.117854, 109.86, 0.063784, 108.403, 0.131471, 110.377, 0.326023, 101.013, 0.127284, 101.488, 0.27798, 110.105, 0.148945, 103.229, 0.14454, 103.893, 0.221119, 104.849, 0.127339, 101.74, 0.156511, 102.527, 0.162979, 103.859, 0.172369, 103.141, 57.563, 1.64141, 54.3487, 1.64141, 0.223534, 1.64141, 0.124233, 1.64141, 0.135591, 1.64141, 0.0763416, 1.64141, 0.151687, 1.64141, 0.253774, 1.64141, 0.0683957, 1.64141, 0.148502, 1.64141, 0.105813, 1.64141, 0.174386, 1.64141, 0.0934222, 1.64141, 0.193053, 1.64141, 0.268681, 1.64141, 0.0993405, 1.64141, 0.211272, 1.64141, 0.11724, 1.64141, 0.112843, 1.64141, 0.175487, 1.64141, 0.0993353, 1.64141, 0.121866, 1.64141, 0.125395, 1.64141, 0.134067, 1.64141, 223.7, 113.546, 223.7, 113.546, 81.4203, -0.00625658, -0.00343448, -0.00367393, -0.00204192, -0.00426901, -0.00667232, -0.00188009, -0.00422262, -0.00273856, -0.00460985, -0.00256085, -0.00507891, -6.67398e-05, 8.96751e-05, 0.0014474, 0.00033538, 0.000384206, 0.000664401, 0.000112356, 0.000156257, 0.000271718, 0.000217284, 35.0055, 1.72547e-12, 45.0271, 68.2839, 555.623, 25.0539, 4.60447, 60.4161, 2.74931e-10, 123.048, 74.8772, 405.595}; - - std::vector init_state_ydot = {-407.383, 603.025, -0.12541, 586.776, -0.143589, 579.533, -0.143206, 573.381, -0.140919, 563.241, -0.117593, 583.876, -0.149131, 559.217, -0.125706, 567.295, -0.111758, 583.808, -0.152064, 563.677, -0.16802, 563.084, -0.123088, 571.513, -0.201867, 564.857, 1.0169, 633.091, 0.273686, 633.771, 0.040692, 533.643, 0.191667, 575.913, 0.19006, 577.538, 0.223164, 553.187, 0.252782, 625.121, 0.35634, 639.502, 0.327164, 633.898, 0.355656, 632.605, 466.061, -19.3723, 464.294, -19.3723, 0.309113, -19.3723, 0.191591, -19.3723, 0.208691, -19.3723, 0.128617, -19.3723, 0.222605, -19.3723, 0.358526, -19.3723, 0.115032, -19.3723, 0.216634, -19.3723, 0.174236, -19.3723, 0.262547, -19.3723, 0.150376, -19.3723, 0.288608, -19.3723, -0.19198, -19.3723, -0.0722613, -19.3723, -0.0502622, -19.3723, -0.0729683, -19.3723, -0.0668817, -19.3723, -0.0913291, -19.3723, -0.070668, -19.3723, -0.0819415, -19.3723, -0.0760433, -19.3723, -0.085503, -19.3723, -404.441, 515.61, -404.441, 515.61, 662.168, -0.139623, -0.0804254, -0.0861598, -0.0501321, -0.0982909, -0.151062, -0.0461619, -0.0972004, -0.0662369, -0.106878, -0.0616797, -0.116491, -0.0378684, -0.0180084, -0.00754156, -0.0171239, -0.0157096, -0.0189088, -0.0175604, -0.01926, -0.0174746, -0.0195021, 57.5594, -115695, -23.5015, -555.659, -4642.24, 257.474, 24.2288, 555.659, -315843, 345.199, -405.655, -7759.96}; - - // Interface blocks flow boundary conditions (neumann boundary conditions for the 3D flow solver) - std::map> interface_block_params = { - {"inlet_BC_RCR1", {220.655, 220.143}}, - {"inlet_BC_lca1", {0.146379, 0.146236}}, - {"inlet_BC_lca10", {0.0840239, 0.0838506}}, - {"inlet_BC_lca11", {0.0917877, 0.0916064}}, - {"inlet_BC_lca12", {0.0539358, 0.0537849}}, - {"inlet_BC_lca2", {0.0997981, 0.0996647}}, - {"inlet_BC_lca3", {0.168397, 0.168252}}, - {"inlet_BC_lca4", {0.0478851, 0.0477658}}, - {"inlet_BC_lca5", {0.0969178, 0.0967786}}, - {"inlet_BC_lca6", {0.0745793, 0.0743957}}, - {"inlet_BC_lca7", {0.117854, 0.117657}}, - {"inlet_BC_lca8", {0.063784, 0.0636423}}, - {"inlet_BC_lca9", {0.131471, 0.131264}}, - {"inlet_BC_rca1", {0.326023, 0.326807}}, - {"inlet_BC_rca10", {0.127284, 0.12748}}, - {"inlet_BC_rca2", {0.27798, 0.278109}}, - {"inlet_BC_rca3", {0.148945, 0.149096}}, - {"inlet_BC_rca4", {0.14454, 0.144655}}, - {"inlet_BC_rca5", {0.221119, 0.221271}}, - {"inlet_BC_rca6", {0.127339, 0.127523}}, - {"inlet_BC_rca7", {0.156511, 0.15675}}, - {"inlet_BC_rca8", {0.162979, 0.163184}}, - {"inlet_BC_rca9", {0.172369, 0.172628}}, - {"outlet_aorta", {223.7, 223.19}}}; + std::vector init_state_y = { + 220.655, 113.454, 0.146379, 107.558, 0.0840239, + 108.31, 0.0917877, 108.966, 0.0539358, 109.893, + 0.0997981, 107.152, 0.168397, 109.693, 0.0478851, + 108.683, 0.0969178, 106.587, 0.0745793, 111.186, + 0.117854, 109.86, 0.063784, 108.403, 0.131471, + 110.377, 0.326023, 101.013, 0.127284, 101.488, + 0.27798, 110.105, 0.148945, 103.229, 0.14454, + 103.893, 0.221119, 104.849, 0.127339, 101.74, + 0.156511, 102.527, 0.162979, 103.859, 0.172369, + 103.141, 57.563, 1.64141, 54.3487, 1.64141, + 0.223534, 1.64141, 0.124233, 1.64141, 0.135591, + 1.64141, 0.0763416, 1.64141, 0.151687, 1.64141, + 0.253774, 1.64141, 0.0683957, 1.64141, 0.148502, + 1.64141, 0.105813, 1.64141, 0.174386, 1.64141, + 0.0934222, 1.64141, 0.193053, 1.64141, 0.268681, + 1.64141, 0.0993405, 1.64141, 0.211272, 1.64141, + 0.11724, 1.64141, 0.112843, 1.64141, 0.175487, + 1.64141, 0.0993353, 1.64141, 0.121866, 1.64141, + 0.125395, 1.64141, 0.134067, 1.64141, 223.7, + 113.546, 223.7, 113.546, 81.4203, -0.00625658, + -0.00343448, -0.00367393, -0.00204192, -0.00426901, -0.00667232, + -0.00188009, -0.00422262, -0.00273856, -0.00460985, -0.00256085, + -0.00507891, -6.67398e-05, 8.96751e-05, 0.0014474, 0.00033538, + 0.000384206, 0.000664401, 0.000112356, 0.000156257, 0.000271718, + 0.000217284, 35.0055, 1.72547e-12, 45.0271, 68.2839, + 555.623, 25.0539, 4.60447, 60.4161, 2.74931e-10, + 123.048, 74.8772, 405.595}; + + std::vector init_state_ydot = { + -407.383, 603.025, -0.12541, 586.776, -0.143589, 579.533, + -0.143206, 573.381, -0.140919, 563.241, -0.117593, 583.876, + -0.149131, 559.217, -0.125706, 567.295, -0.111758, 583.808, + -0.152064, 563.677, -0.16802, 563.084, -0.123088, 571.513, + -0.201867, 564.857, 1.0169, 633.091, 0.273686, 633.771, + 0.040692, 533.643, 0.191667, 575.913, 0.19006, 577.538, + 0.223164, 553.187, 0.252782, 625.121, 0.35634, 639.502, + 0.327164, 633.898, 0.355656, 632.605, 466.061, -19.3723, + 464.294, -19.3723, 0.309113, -19.3723, 0.191591, -19.3723, + 0.208691, -19.3723, 0.128617, -19.3723, 0.222605, -19.3723, + 0.358526, -19.3723, 0.115032, -19.3723, 0.216634, -19.3723, + 0.174236, -19.3723, 0.262547, -19.3723, 0.150376, -19.3723, + 0.288608, -19.3723, -0.19198, -19.3723, -0.0722613, -19.3723, + -0.0502622, -19.3723, -0.0729683, -19.3723, -0.0668817, -19.3723, + -0.0913291, -19.3723, -0.070668, -19.3723, -0.0819415, -19.3723, + -0.0760433, -19.3723, -0.085503, -19.3723, -404.441, 515.61, + -404.441, 515.61, 662.168, -0.139623, -0.0804254, -0.0861598, + -0.0501321, -0.0982909, -0.151062, -0.0461619, -0.0972004, -0.0662369, + -0.106878, -0.0616797, -0.116491, -0.0378684, -0.0180084, -0.00754156, + -0.0171239, -0.0157096, -0.0189088, -0.0175604, -0.01926, -0.0174746, + -0.0195021, 57.5594, -115695, -23.5015, -555.659, -4642.24, + 257.474, 24.2288, 555.659, -315843, 345.199, -405.655, + -7759.96}; + + // Interface blocks flow boundary conditions (neumann boundary conditions for + // the 3D flow solver) + std::map> interface_block_params = { + {"inlet_BC_RCR1", {220.655, 220.143}}, + {"inlet_BC_lca1", {0.146379, 0.146236}}, + {"inlet_BC_lca10", {0.0840239, 0.0838506}}, + {"inlet_BC_lca11", {0.0917877, 0.0916064}}, + {"inlet_BC_lca12", {0.0539358, 0.0537849}}, + {"inlet_BC_lca2", {0.0997981, 0.0996647}}, + {"inlet_BC_lca3", {0.168397, 0.168252}}, + {"inlet_BC_lca4", {0.0478851, 0.0477658}}, + {"inlet_BC_lca5", {0.0969178, 0.0967786}}, + {"inlet_BC_lca6", {0.0745793, 0.0743957}}, + {"inlet_BC_lca7", {0.117854, 0.117657}}, + {"inlet_BC_lca8", {0.063784, 0.0636423}}, + {"inlet_BC_lca9", {0.131471, 0.131264}}, + {"inlet_BC_rca1", {0.326023, 0.326807}}, + {"inlet_BC_rca10", {0.127284, 0.12748}}, + {"inlet_BC_rca2", {0.27798, 0.278109}}, + {"inlet_BC_rca3", {0.148945, 0.149096}}, + {"inlet_BC_rca4", {0.14454, 0.144655}}, + {"inlet_BC_rca5", {0.221119, 0.221271}}, + {"inlet_BC_rca6", {0.127339, 0.127523}}, + {"inlet_BC_rca7", {0.156511, 0.15675}}, + {"inlet_BC_rca8", {0.162979, 0.163184}}, + {"inlet_BC_rca9", {0.172369, 0.172628}}, + {"outlet_aorta", {223.7, 223.19}}}; std::vector interface_times = {0.082147, 0.08307}; // Get variable IDs @@ -96,65 +151,70 @@ int main(int argc, char** argv) double inlet_or_outlet; interface.get_block_node_IDs(block_name, IDs); // IDs in the above function stores info in the following format: - // {num inlet nodes, inlet flow[0], inlet pressure[0],..., num outlet nodes, outlet flow[0], outlet pressure[0],...} + // {num inlet nodes, inlet flow[0], inlet pressure[0],..., num outlet nodes, + // outlet flow[0], outlet pressure[0],...} int num_inlet_nodes = IDs[0]; - int num_outlet_nodes = IDs[1+num_inlet_nodes*2]; + int num_outlet_nodes = IDs[1 + num_inlet_nodes * 2]; if (block_name == "outlet_aorta") { if ((num_inlet_nodes != 1) || (num_outlet_nodes != 0)) { - throw std::runtime_error("Wrong number of inlets/outlets for outlet_aorta"); + throw std::runtime_error( + "Wrong number of inlets/outlets for outlet_aorta"); } } else { if ((num_inlet_nodes != 0) || (num_outlet_nodes != 1)) { - throw std::runtime_error("Wrong number of inlets/outlets for " + block_name); + throw std::runtime_error("Wrong number of inlets/outlets for " + + block_name); } } } // --- For outlet from heart block std::vector IDs; - std::string block_name = "J_heart_outlet"; + std::string block_name = "J_heart_outlet"; interface.get_block_node_IDs(block_name, IDs); int num_inlet_nodes = IDs[0]; - int num_outlet_nodes = IDs[1+num_inlet_nodes*2]; + int num_outlet_nodes = IDs[1 + num_inlet_nodes * 2]; if ((num_inlet_nodes != 1) && (num_outlet_nodes != 1)) { - throw std::runtime_error("Wrong number of inlets/outlets for J_heart_outlet"); + throw std::runtime_error( + "Wrong number of inlets/outlets for J_heart_outlet"); } int aortic_inlet_flow_id = IDs[1]; int aortic_inlet_pressure_id = IDs[2]; // --- For outlet from coronary - block_name = "BC_lca1"; + block_name = "BC_lca1"; interface.get_block_node_IDs(block_name, IDs); num_inlet_nodes = IDs[0]; - num_outlet_nodes = IDs[1+num_inlet_nodes*2]; + num_outlet_nodes = IDs[1 + num_inlet_nodes * 2]; if ((num_inlet_nodes != 1) && (num_outlet_nodes != 1)) { throw std::runtime_error("Wrong number of inlets/outlets for BC_lca1"); } int bc_lca1_outlet_flow_id = IDs[4]; int bc_lca1_outlet_pressure_id = IDs[5]; - + // Update block parameters with current flow from 3D solver for (const auto block_params : interface_block_params) { std::vector new_params(5); std::vector params = block_params.second; - // Format of new_params for flow/pressure blocks: + // Format of new_params for flow/pressure blocks: // [N, time_1, time_2, ..., time_N, value1, value2, ..., value_N] // where N is number of time points and value* is flow/pressure new_params[0] = 2.0; for (int i = 0; i < 2; i++) { - new_params[1+i] = interface_times[i]; - new_params[3+i] = params[i]; + new_params[1 + i] = interface_times[i]; + new_params[3 + i] = params[i]; } - interface.update_block_params(block_params.first, new_params); + interface.update_block_params(block_params.first, new_params); } - + // Set up vectors to run svZeroD simulation - std::vector solutions(interface.system_size_*interface.num_output_steps_); + std::vector solutions(interface.system_size_ * + interface.num_output_steps_); std::vector times(interface.num_output_steps_); int error_code = 0; - + // Run svZeroD simulation - interface.update_state(init_state_y, init_state_ydot); + interface.update_state(init_state_y, init_state_ydot); interface.run_simulation(0.0, times, solutions, error_code); - + // Parse output and calculate mean flow/pressure in aorta and coronary int sol_idx = 0; double mean_aortic_flow = 0.0; @@ -163,7 +223,7 @@ int main(int argc, char** argv) double mean_bc_lca1_outlet_pressure = 0.0; for (int tstep = 0; tstep < interface.num_output_steps_; tstep++) { for (int state = 0; state < interface.system_size_; state++) { - sol_idx = interface.system_size_*tstep + state; + sol_idx = interface.system_size_ * tstep + state; if (state == aortic_inlet_flow_id) { mean_aortic_flow += solutions[sol_idx]; } else if (state == aortic_inlet_pressure_id) { @@ -180,14 +240,17 @@ int main(int argc, char** argv) mean_aortic_pressure /= (double)interface.num_output_steps_; mean_bc_lca1_outlet_flow /= (double)interface.num_output_steps_; mean_bc_lca1_outlet_pressure /= (double)interface.num_output_steps_; - - std::cout <<"Simulation output: " << std::endl; - std::cout <<"Mean aortic flow = " << mean_aortic_flow << std::endl; - std::cout <<"Mean aortic pressure = " << mean_aortic_pressure << std::endl; - std::cout <<"Mean BC_lca1 outlet flow = " << mean_bc_lca1_outlet_flow << std::endl; - std::cout <<"Mean BC_lca1 outlet pressure = " << mean_bc_lca1_outlet_pressure << std::endl; - - // Compare mean flow/pressure in aorta and coronary with pre-computed ("correct") values + + std::cout << "Simulation output: " << std::endl; + std::cout << "Mean aortic flow = " << mean_aortic_flow << std::endl; + std::cout << "Mean aortic pressure = " << mean_aortic_pressure << std::endl; + std::cout << "Mean BC_lca1 outlet flow = " << mean_bc_lca1_outlet_flow + << std::endl; + std::cout << "Mean BC_lca1 outlet pressure = " << mean_bc_lca1_outlet_pressure + << std::endl; + + // Compare mean flow/pressure in aorta and coronary with pre-computed + // ("correct") values double error_limit = 0.05; std::vector wrong_quantities; bool is_wrong = false; @@ -207,14 +270,14 @@ int main(int argc, char** argv) is_wrong = true; wrong_quantities.push_back("Mean BC_lca1 outlet pressure"); } - + if (is_wrong) { std::string error_msg = " "; for (int i = 0; i < wrong_quantities.size(); i++) { - error_msg = error_msg + wrong_quantities[i] + "; "; + error_msg = error_msg + wrong_quantities[i] + "; "; } throw std::runtime_error("Error in the following quantities:" + error_msg); } - + return 0; } diff --git a/tests/test_interface/test_02/main.cpp b/tests/test_interface/test_02/main.cpp index 57bd803d7..a53330761 100644 --- a/tests/test_interface/test_02/main.cpp +++ b/tests/test_interface/test_02/main.cpp @@ -1,33 +1,49 @@ // Test interfacing to svZeroDSolver. -// This test mimics the coupling of svZeroDSolver with an external parameter estimation code (eg. Tulip) -// A coronary model is used and parameters of the BCs are updated and then compared to a reference solution +// This test mimics the coupling of svZeroDSolver with an external parameter +// estimation code (eg. Tulip) A coronary model is used and parameters of the +// BCs are updated and then compared to a reference solution -#include "../LPNSolverInterface/LPNSolverInterface.h" +#include +#include #include #include -#include -#include + +#include "../LPNSolverInterface/LPNSolverInterface.h" namespace fs = std::filesystem; //--------------------------------------------------------------------------------------- -// Compare mean flow/pressure in aorta and coronary with pre-computed ("correct") values +// Compare mean flow/pressure in aorta and coronary with pre-computed +// ("correct") values //--------------------------------------------------------------------------------------- -std::string check_simulation_results(double mean_aortic_flow, double mean_aortic_pressure, double mean_bc_lca1_outlet_flow, double mean_bc_lca1_outlet_pressure, double correct_mean_aortic_flow, double correct_mean_aortic_pressure, double correct_mean_bc_lca1_outlet_flow, double correct_mean_bc_lca1_outlet_pressure, bool &is_wrong) { +std::string check_simulation_results( + double mean_aortic_flow, double mean_aortic_pressure, + double mean_bc_lca1_outlet_flow, double mean_bc_lca1_outlet_pressure, + double correct_mean_aortic_flow, double correct_mean_aortic_pressure, + double correct_mean_bc_lca1_outlet_flow, + double correct_mean_bc_lca1_outlet_pressure, bool& is_wrong) { double error_limit = 0.01; std::vector wrong_quantities; - if (abs(mean_aortic_flow - correct_mean_aortic_flow)/correct_mean_aortic_flow > error_limit) { + if (abs(mean_aortic_flow - correct_mean_aortic_flow) / + correct_mean_aortic_flow > + error_limit) { is_wrong = true; wrong_quantities.push_back("Mean aortic flow"); } - if (abs(mean_aortic_pressure - correct_mean_aortic_pressure)/correct_mean_aortic_pressure > error_limit) { + if (abs(mean_aortic_pressure - correct_mean_aortic_pressure) / + correct_mean_aortic_pressure > + error_limit) { is_wrong = true; wrong_quantities.push_back("Mean aortic pressure"); } - if (abs(mean_bc_lca1_outlet_flow - correct_mean_bc_lca1_outlet_flow)/correct_mean_bc_lca1_outlet_flow > error_limit) { + if (abs(mean_bc_lca1_outlet_flow - correct_mean_bc_lca1_outlet_flow) / + correct_mean_bc_lca1_outlet_flow > + error_limit) { is_wrong = true; wrong_quantities.push_back("Mean BC_lca1 outlet flow"); } - if (abs(mean_bc_lca1_outlet_pressure - correct_mean_bc_lca1_outlet_pressure)/correct_mean_bc_lca1_outlet_pressure > error_limit) { + if (abs(mean_bc_lca1_outlet_pressure - correct_mean_bc_lca1_outlet_pressure) / + correct_mean_bc_lca1_outlet_pressure > + error_limit) { is_wrong = true; wrong_quantities.push_back("Mean BC_lca1 outlet pressure"); } @@ -40,22 +56,22 @@ std::string check_simulation_results(double mean_aortic_flow, double mean_aortic return error_msg; } - //------ // main //------ -int main(int argc, char** argv) -{ +int main(int argc, char** argv) { LPNSolverInterface interface; - + if (argc != 3) { - std::runtime_error("Usage: svZeroD_interface_test01 "); + std::runtime_error( + "Usage: svZeroD_interface_test01 " + ""); } // Load shared library and get interface functions. // File extension of the shared library depends on the system - fs::path build_dir = argv[1]; - fs::path iface_dir = build_dir / "src" / "interface"; + 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"; @@ -66,7 +82,9 @@ int main(int argc, char** argv) } 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() + " !"); + throw std::runtime_error("Could not find shared libraries " + + lib_so.string() + " or " + lib_dylib.string() + + " or " + lib_dll.string() + " !"); } // Set up the svZeroD model @@ -84,42 +102,45 @@ int main(int argc, char** argv) // Save important variable IDs // --- For outlet from heart block std::vector IDs; - std::string block_name = "J_heart_outlet"; + std::string block_name = "J_heart_outlet"; interface.get_block_node_IDs(block_name, IDs); // IDs in the above function stores info in the following format: - // {num inlet nodes, inlet flow[0], inlet pressure[0],..., num outlet nodes, outlet flow[0], outlet pressure[0],...} + // {num inlet nodes, inlet flow[0], inlet pressure[0],..., num outlet nodes, + // outlet flow[0], outlet pressure[0],...} int num_inlet_nodes = IDs[0]; - int num_outlet_nodes = IDs[1+num_inlet_nodes*2]; + int num_outlet_nodes = IDs[1 + num_inlet_nodes * 2]; if ((num_inlet_nodes != 1) && (num_outlet_nodes != 1)) { - throw std::runtime_error("Wrong number of inlets/outlets for J_heart_outlet"); + throw std::runtime_error( + "Wrong number of inlets/outlets for J_heart_outlet"); } int aortic_inlet_flow_id = IDs[1]; int aortic_inlet_pressure_id = IDs[2]; // --- For outlet from coronary - block_name = "BC_lca1"; + block_name = "BC_lca1"; interface.get_block_node_IDs(block_name, IDs); num_inlet_nodes = IDs[0]; - num_outlet_nodes = IDs[1+num_inlet_nodes*2]; + num_outlet_nodes = IDs[1 + num_inlet_nodes * 2]; if ((num_inlet_nodes != 1) && (num_outlet_nodes != 1)) { throw std::runtime_error("Wrong number of inlets/outlets for BC_lca1"); } int bc_lca1_outlet_flow_id = IDs[4]; int bc_lca1_outlet_pressure_id = IDs[5]; - + // Save the initial condition - std::vector init_state_y, init_state_ydot; + std::vector init_state_y, init_state_ydot; init_state_y.resize(interface.system_size_); init_state_ydot.resize(interface.system_size_); interface.return_y(init_state_y); interface.return_ydot(init_state_ydot); // Set up vectors to run svZeroD simulation - std::vector solutions(interface.system_size_*interface.num_output_steps_); + std::vector solutions(interface.system_size_ * + interface.num_output_steps_); std::vector times(interface.num_output_steps_); int error_code = 0; - + // Run svZeroD simulation - interface.update_state(init_state_y, init_state_ydot); + interface.update_state(init_state_y, init_state_ydot); interface.run_simulation(0.0, times, solutions, error_code); // Parse and print output @@ -136,7 +157,7 @@ int main(int argc, char** argv) double mean_bc_lca1_outlet_pressure = 0.0; for (int tstep = 0; tstep < interface.num_output_steps_; tstep++) { for (int state = 0; state < interface.system_size_; state++) { - sol_idx = interface.system_size_*tstep + state; + sol_idx = interface.system_size_ * tstep + state; if (state == aortic_inlet_flow_id) { aortic_flow[tstep] = solutions[sol_idx]; mean_aortic_flow += solutions[sol_idx]; @@ -156,45 +177,52 @@ int main(int argc, char** argv) mean_aortic_pressure /= (double)interface.num_output_steps_; mean_bc_lca1_outlet_flow /= (double)interface.num_output_steps_; mean_bc_lca1_outlet_pressure /= (double)interface.num_output_steps_; - + std::cout << "First simulation: " << std::endl; std::cout << "Mean aortic flow = " << mean_aortic_flow << std::endl; std::cout << "Mean aortic pressure = " << mean_aortic_pressure << std::endl; - std::cout << "Mean BC_lca1 outlet flow = " << mean_bc_lca1_outlet_flow << std::endl; - std::cout << "Mean BC_lca1 outlet pressure = " << mean_bc_lca1_outlet_pressure << std::endl; + std::cout << "Mean BC_lca1 outlet flow = " << mean_bc_lca1_outlet_flow + << std::endl; + std::cout << "Mean BC_lca1 outlet pressure = " << mean_bc_lca1_outlet_pressure + << std::endl; // Check if outputs are correct bool is_wrong = false; std::string error_msg; - error_msg = check_simulation_results(mean_aortic_flow, mean_aortic_pressure, mean_bc_lca1_outlet_flow, mean_bc_lca1_outlet_pressure, 63.3137, 101.139, 0.135942, 3.17525, is_wrong); + error_msg = check_simulation_results(mean_aortic_flow, mean_aortic_pressure, + mean_bc_lca1_outlet_flow, + mean_bc_lca1_outlet_pressure, 63.3137, + 101.139, 0.135942, 3.17525, is_wrong); if (is_wrong) { - throw std::runtime_error("After initial simulation, error in the following quantities: "+error_msg); + throw std::runtime_error( + "After initial simulation, error in the following quantities: " + + error_msg); } // Save the initial coronary params - std::map> initial_coronary_params = { - {"BC_lca1", {145.272, 236.068, 290.065, 0.00010326, 0.00105039}}, - {"BC_lca10", {264.502, 429.816, 528.13, 6.513e-05, 0.00066247}}, - {"BC_lca11", {242.343, 393.807, 483.885, 6.966e-05, 0.00070852}}, - {"BC_lca12", {434.805, 706.557, 868.172, 4.44e-05, 0.00045194}}, - {"BC_lca2", {215.568, 350.299, 430.424, 7.617e-05, 0.00077533}}, - {"BC_lca3", {128.321, 208.522, 256.219, 0.00011358, 0.00115556}}, - {"BC_lca4", {485.259, 788.546, 968.915, 4.083e-05, 0.00041538}}, - {"BC_lca5", {220.097, 357.658, 439.467, 7.498e-05, 0.00076305}}, - {"BC_lca6", {312.826, 508.343, 624.619, 5.727e-05, 0.00058227}}, - {"BC_lca7", {187.955, 305.427, 375.289, 8.467e-05, 0.00086153}}, - {"BC_lca8", {353.619, 574.63, 706.069, 5.21e-05, 0.00052984}}, - {"BC_lca9", {169.536, 275.495, 338.511, 9.166e-05, 0.00093274}}, - {"BC_rca1", {80.2918, 130.474, 160.318, 0.00017264, 0.00141371}}, - {"BC_rca10", {218.089, 354.394, 435.457, 8.004e-05, 0.00065544}}, - {"BC_rca2", {105.209, 170.965, 210.07, 0.00014026, 0.00114837}}, - {"BC_rca3", {186.143, 302.482, 371.671, 9.038e-05, 0.00074037}}, - {"BC_rca4", {193.761, 314.861, 386.881, 8.767e-05, 0.00071786}}, - {"BC_rca5", {124.951, 203.045, 249.489, 0.00012286, 0.0010061}}, - {"BC_rca6", {218.275, 354.697, 435.829, 7.994e-05, 0.00065505}}, - {"BC_rca7", {178.053, 289.337, 355.519, 9.357e-05, 0.0007661}}, - {"BC_rca8", {173.617, 282.127, 346.66, 9.54e-05, 0.00078117}}, - {"BC_rca9", {162.073, 263.368, 323.61, 0.00010053, 0.00082363}}}; + std::map> initial_coronary_params = { + {"BC_lca1", {145.272, 236.068, 290.065, 0.00010326, 0.00105039}}, + {"BC_lca10", {264.502, 429.816, 528.13, 6.513e-05, 0.00066247}}, + {"BC_lca11", {242.343, 393.807, 483.885, 6.966e-05, 0.00070852}}, + {"BC_lca12", {434.805, 706.557, 868.172, 4.44e-05, 0.00045194}}, + {"BC_lca2", {215.568, 350.299, 430.424, 7.617e-05, 0.00077533}}, + {"BC_lca3", {128.321, 208.522, 256.219, 0.00011358, 0.00115556}}, + {"BC_lca4", {485.259, 788.546, 968.915, 4.083e-05, 0.00041538}}, + {"BC_lca5", {220.097, 357.658, 439.467, 7.498e-05, 0.00076305}}, + {"BC_lca6", {312.826, 508.343, 624.619, 5.727e-05, 0.00058227}}, + {"BC_lca7", {187.955, 305.427, 375.289, 8.467e-05, 0.00086153}}, + {"BC_lca8", {353.619, 574.63, 706.069, 5.21e-05, 0.00052984}}, + {"BC_lca9", {169.536, 275.495, 338.511, 9.166e-05, 0.00093274}}, + {"BC_rca1", {80.2918, 130.474, 160.318, 0.00017264, 0.00141371}}, + {"BC_rca10", {218.089, 354.394, 435.457, 8.004e-05, 0.00065544}}, + {"BC_rca2", {105.209, 170.965, 210.07, 0.00014026, 0.00114837}}, + {"BC_rca3", {186.143, 302.482, 371.671, 9.038e-05, 0.00074037}}, + {"BC_rca4", {193.761, 314.861, 386.881, 8.767e-05, 0.00071786}}, + {"BC_rca5", {124.951, 203.045, 249.489, 0.00012286, 0.0010061}}, + {"BC_rca6", {218.275, 354.697, 435.829, 7.994e-05, 0.00065505}}, + {"BC_rca7", {178.053, 289.337, 355.519, 9.357e-05, 0.0007661}}, + {"BC_rca8", {173.617, 282.127, 346.66, 9.54e-05, 0.00078117}}, + {"BC_rca9", {162.073, 263.368, 323.61, 0.00010053, 0.00082363}}}; // Check if correct parameters are read and then update block parameters double param_update_factor = 2.0; @@ -203,33 +231,38 @@ int main(int argc, char** argv) coronary_params.resize(num_coronary_params); for (int i = 0; i < interface.block_names_.size(); i++) { std::string block_name = interface.block_names_[i]; - if ((block_name.substr(0,6) == "BC_lca") || (block_name.substr(0,6) == "BC_rca")) { + if ((block_name.substr(0, 6) == "BC_lca") || + (block_name.substr(0, 6) == "BC_rca")) { // Read current block parameters and compare with data above interface.read_block_params(block_name, coronary_params); auto correct_params = initial_coronary_params[block_name]; for (int j = 0; j < num_coronary_params; j++) { - if(abs(coronary_params[j]-correct_params[j])/correct_params[j] > 0.01) { - throw std::runtime_error("Wrong parameters read from block " + block_name); + if (abs(coronary_params[j] - correct_params[j]) / correct_params[j] > + 0.01) { + throw std::runtime_error("Wrong parameters read from block " + + block_name); } } // Update parameters by multiplying Ra, Ram and Rv by param_update_factor - coronary_params[0] *= param_update_factor; - coronary_params[1] *= param_update_factor; - coronary_params[2] *= param_update_factor; + coronary_params[0] *= param_update_factor; + coronary_params[1] *= param_update_factor; + coronary_params[2] *= param_update_factor; interface.update_block_params(block_name, coronary_params); // Verify that parameters were correctly updated interface.read_block_params(block_name, coronary_params); for (int j = 0; j < 3; j++) { correct_params[j] *= param_update_factor; - if(abs(coronary_params[j]-correct_params[j])/correct_params[j] > 0.01) { - throw std::runtime_error("Wrong parameters after update for block " + block_name); + if (abs(coronary_params[j] - correct_params[j]) / correct_params[j] > + 0.01) { + throw std::runtime_error("Wrong parameters after update for block " + + block_name); } } } } - + // Run svZeroD simulation with new parameters - interface.update_state(init_state_y, init_state_ydot); + interface.update_state(init_state_y, init_state_ydot); interface.run_simulation(0.0, times, solutions, error_code); // Parse and print output @@ -240,7 +273,7 @@ int main(int argc, char** argv) mean_bc_lca1_outlet_pressure = 0.0; for (int tstep = 0; tstep < interface.num_output_steps_; tstep++) { for (int state = 0; state < interface.system_size_; state++) { - sol_idx = interface.system_size_*tstep + state; + sol_idx = interface.system_size_ * tstep + state; if (state == aortic_inlet_flow_id) { aortic_flow[tstep] = solutions[sol_idx]; mean_aortic_flow += solutions[sol_idx]; @@ -260,45 +293,56 @@ int main(int argc, char** argv) mean_aortic_pressure /= (double)interface.num_output_steps_; mean_bc_lca1_outlet_flow /= (double)interface.num_output_steps_; mean_bc_lca1_outlet_pressure /= (double)interface.num_output_steps_; - + std::cout << "After parameter update: " << std::endl; std::cout << "Mean aortic flow = " << mean_aortic_flow << std::endl; std::cout << "Mean aortic pressure = " << mean_aortic_pressure << std::endl; - std::cout << "Mean BC_lca1 outlet flow = " << mean_bc_lca1_outlet_flow << std::endl; - std::cout << "Mean BC_lca1 outlet pressure = " << mean_bc_lca1_outlet_pressure << std::endl; - + std::cout << "Mean BC_lca1 outlet flow = " << mean_bc_lca1_outlet_flow + << std::endl; + std::cout << "Mean BC_lca1 outlet pressure = " << mean_bc_lca1_outlet_pressure + << std::endl; + // Check if outputs are correct is_wrong = false; - error_msg = check_simulation_results(mean_aortic_flow, mean_aortic_pressure, mean_bc_lca1_outlet_flow, mean_bc_lca1_outlet_pressure, 63.2715, 101.232, 0.0691725, 3.17279, is_wrong); + error_msg = check_simulation_results(mean_aortic_flow, mean_aortic_pressure, + mean_bc_lca1_outlet_flow, + mean_bc_lca1_outlet_pressure, 63.2715, + 101.232, 0.0691725, 3.17279, is_wrong); if (is_wrong) { - throw std::runtime_error("After parameter update, error in the following quantities: "+error_msg); + throw std::runtime_error( + "After parameter update, error in the following quantities: " + + error_msg); } - - // Test restarting simulation by overwriting state with prescribed initial condition and restoring initial parameters - // Check if correct parameters are read and then update block parameters + + // Test restarting simulation by overwriting state with prescribed initial + // condition and restoring initial parameters Check if correct parameters are + // read and then update block parameters for (int i = 0; i < interface.block_names_.size(); i++) { std::string block_name = interface.block_names_[i]; - if ((block_name.substr(0,6) == "BC_lca") || (block_name.substr(0,6) == "BC_rca")) { + if ((block_name.substr(0, 6) == "BC_lca") || + (block_name.substr(0, 6) == "BC_rca")) { // Read current block parameters and compare with data above interface.read_block_params(block_name, coronary_params); auto initial_params = initial_coronary_params[block_name]; // Update parameters by dividing Ra, Ram and Rv by param_update_factor - coronary_params[0] /= param_update_factor; - coronary_params[1] /= param_update_factor; - coronary_params[2] /= param_update_factor; + coronary_params[0] /= param_update_factor; + coronary_params[1] /= param_update_factor; + coronary_params[2] /= param_update_factor; interface.update_block_params(block_name, coronary_params); // Verify that parameters were correctly updated interface.read_block_params(block_name, coronary_params); for (int j = 0; j < 3; j++) { - if(abs(coronary_params[j]-initial_params[j])/initial_params[j] > 0.01) { - throw std::runtime_error("Wrong parameters after update for block " + block_name); + if (abs(coronary_params[j] - initial_params[j]) / initial_params[j] > + 0.01) { + throw std::runtime_error("Wrong parameters after update for block " + + block_name); } } } } // Restore initial state and run simulation - interface.update_state(init_state_y, init_state_ydot); + interface.update_state(init_state_y, init_state_ydot); interface.run_simulation(0.0, times, solutions, error_code); // Parse and print output @@ -309,7 +353,7 @@ int main(int argc, char** argv) mean_bc_lca1_outlet_pressure = 0.0; for (int tstep = 0; tstep < interface.num_output_steps_; tstep++) { for (int state = 0; state < interface.system_size_; state++) { - sol_idx = interface.system_size_*tstep + state; + sol_idx = interface.system_size_ * tstep + state; if (state == aortic_inlet_flow_id) { aortic_flow[tstep] = solutions[sol_idx]; mean_aortic_flow += solutions[sol_idx]; @@ -329,19 +373,26 @@ int main(int argc, char** argv) mean_aortic_pressure /= (double)interface.num_output_steps_; mean_bc_lca1_outlet_flow /= (double)interface.num_output_steps_; mean_bc_lca1_outlet_pressure /= (double)interface.num_output_steps_; - - std::cout <<"After restart with same parameters: " << std::endl; - std::cout <<"Mean aortic flow = " << mean_aortic_flow << std::endl; - std::cout <<"Mean aortic pressure = " << mean_aortic_pressure << std::endl; - std::cout <<"Mean BC_lca1 outlet flow = " << mean_bc_lca1_outlet_flow << std::endl; - std::cout <<"Mean BC_lca1 outlet pressure = " << mean_bc_lca1_outlet_pressure << std::endl; + + std::cout << "After restart with same parameters: " << std::endl; + std::cout << "Mean aortic flow = " << mean_aortic_flow << std::endl; + std::cout << "Mean aortic pressure = " << mean_aortic_pressure << std::endl; + std::cout << "Mean BC_lca1 outlet flow = " << mean_bc_lca1_outlet_flow + << std::endl; + std::cout << "Mean BC_lca1 outlet pressure = " << mean_bc_lca1_outlet_pressure + << std::endl; // Check if outputs are correct is_wrong = false; - error_msg = check_simulation_results(mean_aortic_flow, mean_aortic_pressure, mean_bc_lca1_outlet_flow, mean_bc_lca1_outlet_pressure, 63.329, 101.158, 0.135971, 3.17321, is_wrong); + error_msg = check_simulation_results(mean_aortic_flow, mean_aortic_pressure, + mean_bc_lca1_outlet_flow, + mean_bc_lca1_outlet_pressure, 63.329, + 101.158, 0.135971, 3.17321, is_wrong); if (is_wrong) { - throw std::runtime_error("After restart simulation, error in the following quantities: "+error_msg); + throw std::runtime_error( + "After restart simulation, error in the following quantities: " + + error_msg); } - + return 0; } diff --git a/tests/test_interface/test_03/main.cpp b/tests/test_interface/test_03/main.cpp index be47ab45d..e0a3c62ff 100644 --- a/tests/test_interface/test_03/main.cpp +++ b/tests/test_interface/test_03/main.cpp @@ -1,31 +1,33 @@ // Test interfacing to svZeroSolver. -// This test mimics an external 3D solver (svSolver/svFSI) interfacing with svZeroDSolver -// The model consists of an RCR BC which acts as a Neumann BC for an external solver -// It mimics two consecutive time steps of an external solver +// This test mimics an external 3D solver (svSolver/svFSI) interfacing with +// svZeroDSolver The model consists of an RCR BC which acts as a Neumann BC for +// an external solver It mimics two consecutive time steps of an external solver -#include "../LPNSolverInterface/LPNSolverInterface.h" +#include +#include #include #include -#include -#include + +#include "../LPNSolverInterface/LPNSolverInterface.h" namespace fs = std::filesystem; //------ // main //------ // -int main(int argc, char** argv) -{ +int main(int argc, char** argv) { LPNSolverInterface interface; if (argc != 3) { - std::runtime_error("Usage: svZeroD_interface_test03 "); + std::runtime_error( + "Usage: svZeroD_interface_test03 " + ""); } - + // Load shared library and get interface functions. // File extension of the shared library depends on the system - fs::path build_dir = argv[1]; - fs::path iface_dir = build_dir / "src" / "interface"; + 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"; @@ -36,13 +38,15 @@ int main(int argc, char** argv) } 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() + " !"); + throw std::runtime_error("Could not find shared libraries " + + lib_so.string() + " or " + lib_dylib.string() + + " or " + lib_dll.string() + " !"); } // Set up the svZeroD model std::string file_name = std::string(argv[2]); interface.initialize(file_name); - + // Check number of variables and blocks if (interface.system_size_ != 3) { throw std::runtime_error("interface.system_size_ != 3"); @@ -54,61 +58,68 @@ int main(int argc, char** argv) // Set external time step size (flow solver step size) double external_step_size = 0.005; interface.set_external_step_size(external_step_size); - + // Save the initial condition - std::vector init_state_y = {-6.2506662304695681e+01, -3.8067539421845140e+04, -3.0504233282976966e+04}; - std::vector init_state_ydot = {-3.0873806830951793e+01, -2.5267653962355386e+05, -2.4894080899699836e+05}; + std::vector init_state_y = {-6.2506662304695681e+01, + -3.8067539421845140e+04, + -3.0504233282976966e+04}; + std::vector init_state_ydot = {-3.0873806830951793e+01, + -2.5267653962355386e+05, + -2.4894080899699836e+05}; // Get variable IDs for inlet to RCR block std::vector IDs; - std::string block_name = "RCR"; + std::string block_name = "RCR"; interface.get_block_node_IDs(block_name, IDs); int num_inlet_nodes = IDs[0]; - int num_outlet_nodes = IDs[1+num_inlet_nodes*2]; + int num_outlet_nodes = IDs[1 + num_inlet_nodes * 2]; if ((num_inlet_nodes != 1) || (num_outlet_nodes != 0)) { throw std::runtime_error("Wrong number of inlets/outlets for RCR"); } int rcr_inlet_flow_id = IDs[1]; int rcr_inlet_pressure_id = IDs[2]; - + // Update block parameters with current flow from 3D solver std::vector new_params(5); - std::vector params = {-6.2506662041472836e+01, -6.2599344518688739e+01}; - std::vector interface_times = {1.9899999999999796e+00, 1.9949999999999795e+00}; - // Format of new_params for flow/pressure blocks: + std::vector params = {-6.2506662041472836e+01, + -6.2599344518688739e+01}; + std::vector interface_times = {1.9899999999999796e+00, + 1.9949999999999795e+00}; + // Format of new_params for flow/pressure blocks: // [N, time_1, time_2, ..., time_N, value1, value2, ..., value_N] // where N is number of time points and value* is flow/pressure new_params[0] = 2.0; for (int i = 0; i < 2; i++) { - new_params[1+i] = interface_times[i]; - new_params[3+i] = params[i]; + new_params[1 + i] = interface_times[i]; + new_params[3 + i] = params[i]; } - interface.update_block_params("RCR_coupling", new_params); - + interface.update_block_params("RCR_coupling", new_params); + // Set up vectors to run svZeroD simulation - std::vector solutions(interface.system_size_*interface.num_output_steps_); + std::vector solutions(interface.system_size_ * + interface.num_output_steps_); std::vector times(interface.num_output_steps_); int error_code = 0; - + // Run svZeroD simulation - interface.update_state(init_state_y, init_state_ydot); + interface.update_state(init_state_y, init_state_ydot); interface.run_simulation(interface_times[0], times, solutions, error_code); - + // Parse output and calculate mean flow/pressure in aorta and coronary int sol_idx = 0; double mean_pressure = 0.0; for (int tstep = 0; tstep < interface.num_output_steps_; tstep++) { for (int state = 0; state < interface.system_size_; state++) { - sol_idx = interface.system_size_*tstep + state; + sol_idx = interface.system_size_ * tstep + state; if (state == rcr_inlet_pressure_id) { mean_pressure += solutions[sol_idx]; } } } mean_pressure /= (double)interface.num_output_steps_; - std::cout <<"Simulation output: " << std::endl; - std::cout <<"Mean pressure = " << mean_pressure << std::endl; - + std::cout << "Simulation output: " << std::endl; + std::cout << "Mean pressure = " << mean_pressure << std::endl; + // Compare mean pressure with pre-computed ("correct") values double error_limit = 0.05; if (abs(-mean_pressure / 38690.2 - 1.0) > error_limit) { @@ -123,29 +134,29 @@ int main(int argc, char** argv) params = {-6.2599344283486374e+01, -6.2630248751732964e+01}; interface_times = {1.9949999999999795e+00, 1.9999999999999793e+00}; for (int i = 0; i < 2; i++) { - new_params[1+i] = interface_times[i]; - new_params[3+i] = params[i]; + new_params[1 + i] = interface_times[i]; + new_params[3 + i] = params[i]; } - interface.update_block_params("RCR_coupling", new_params); - + interface.update_block_params("RCR_coupling", new_params); + // Run svZeroD simulation - interface.update_state(init_state_y, init_state_ydot); + interface.update_state(init_state_y, init_state_ydot); interface.run_simulation(interface_times[0], times, solutions, error_code); - + // Parse output and calculate mean flow/pressure in aorta and coronary sol_idx = 0; mean_pressure = 0.0; for (int tstep = 0; tstep < interface.num_output_steps_; tstep++) { for (int state = 0; state < interface.system_size_; state++) { - sol_idx = interface.system_size_*tstep + state; + sol_idx = interface.system_size_ * tstep + state; if (state == rcr_inlet_pressure_id) { mean_pressure += solutions[sol_idx]; } } } mean_pressure /= (double)interface.num_output_steps_; - std::cout <<"Simulation output: " << std::endl; - std::cout <<"Mean pressure = " << mean_pressure << std::endl; + std::cout << "Simulation output: " << std::endl; + std::cout << "Mean pressure = " << mean_pressure << std::endl; // Compare mean pressure with pre-computed ("correct") values if (abs(-mean_pressure / 39911.3 - 1.0) > error_limit) { diff --git a/tests/test_solver.py b/tests/test_solver.py index 74bbf27be..d38315c9b 100644 --- a/tests/test_solver.py +++ b/tests/test_solver.py @@ -61,8 +61,8 @@ def test_solver(testfile): rtol_pres = RTOL_PRES rtol_flow = RTOL_FLOW if 'coupledBlock_closedLoopHeart_withCoronaries.json' in testfile: - rtol_pres = 1.0e-1 - rtol_flow = 1.0e-1 + rtol_pres = 2.0e-1 + rtol_flow = 2.0e-1 this_file_dir = os.path.abspath(os.path.dirname(__file__))