From 7d8afbf5a90226403dde6a228aa4df5a2281b57d Mon Sep 17 00:00:00 2001 From: Jakob Richter <103637279+richterjakob@users.noreply.github.com> Date: Wed, 20 Sep 2023 23:45:16 +0200 Subject: [PATCH] Remove python (#34) * Remove python folder * Delete python test * Make simulate function return a pandas dataframe * Cleanup tests and dependencies * Remove coverage settings for pytest * Remove python doc * Move cpp docs one level up * Update doxyfile * Update Readme * Update github actions * Remove unnecessary entries from gitignore * Fix bug * Add stub file * Fix tests * Test cases for shared library interface (#33) * added test case and included in github/workflows/test.yml * modified test * included dylib extension for macos * clean up * added test2 * updated library path * add cmake build to github workflows * updated docs * more docs updates * added svzerdcalibrator to features in doc * cmake build for interface tests * tests use common LPNSolverInterface * added to github workflows --------- Co-authored-by: Karthik Menon Co-authored-by: Karthik Menon Co-authored-by: Karthik Menon * Remove python folder * Delete python test * Make simulate function return a pandas dataframe * Cleanup tests and dependencies * Remove coverage settings for pytest * Remove python doc * Move cpp docs one level up * Update doxyfile * Update Readme * Update github actions * Remove unnecessary entries from gitignore * Fix bug * Add stub file * Fix tests * Fix bug * Update documentation * Fix bug --------- Co-authored-by: Jakob Richter Co-authored-by: Karthik Menon <40070586+menon-karthik@users.noreply.github.com> Co-authored-by: Karthik Menon Co-authored-by: Karthik Menon Co-authored-by: Karthik Menon Co-authored-by: Martin R. Pfaller --- .github/workflows/documentation.yml | 21 +- .github/workflows/test.yml | 37 +- .gitignore | 6 +- README.md | 11 +- applications/svzerodplus.cpp | 22 +- docs/{cpp => }/Doxyfile | 16 +- docs/{cpp => }/DoxygenLayout.xml | 0 docs/{cpp => }/doxygen-awesome.css | 0 docs/{cpp => }/pages/developer_guide.md | 2 +- docs/{cpp => }/pages/main.md | 64 ++- docs/python/source/conf.py | 41 -- docs/python/source/index.rst | 98 ---- docs/python/source/modules.rst | 7 - docs/python/source/refs.bib | 23 - docs/python/source/rst/bibliography.rst | 7 - docs/{cpp => }/references.bib | 0 pyproject.toml | 3 - setup.cfg | 16 +- src/interface/interface.cpp | 2 +- src/io/csvwriter.hpp | 2 +- svzerodplus.pyi | 95 ++++ svzerodsolver/__init__.py | 34 -- svzerodsolver/algebra.py | 210 -------- svzerodsolver/model/__init__.py | 29 -- svzerodsolver/model/block.py | 164 ------ svzerodsolver/model/bloodvessel.py | 110 ----- svzerodsolver/model/bloodvesseljunction.py | 45 -- svzerodsolver/model/dofhandler.py | 77 --- svzerodsolver/model/flowreferencebc.py | 94 ---- svzerodsolver/model/internaljunction.py | 94 ---- svzerodsolver/model/node.py | 73 --- svzerodsolver/model/openloopcoronarybc.py | 175 ------- svzerodsolver/model/pressurereferencebc.py | 94 ---- svzerodsolver/model/resistancebc.py | 104 ---- svzerodsolver/model/windkesselbc.py | 128 ----- svzerodsolver/runner.py | 158 ------ svzerodsolver/runnercpp.py | 48 -- svzerodsolver/utils.py | 270 ---------- tests/conftest.py | 11 - tests/test_integration.py | 488 ++++++++++++------ tests/test_integration_cpp.py | 550 --------------------- 41 files changed, 529 insertions(+), 2900 deletions(-) rename docs/{cpp => }/Doxyfile (96%) rename docs/{cpp => }/DoxygenLayout.xml (100%) rename docs/{cpp => }/doxygen-awesome.css (100%) rename docs/{cpp => }/pages/developer_guide.md (99%) rename docs/{cpp => }/pages/main.md (84%) delete mode 100644 docs/python/source/conf.py delete mode 100644 docs/python/source/index.rst delete mode 100644 docs/python/source/modules.rst delete mode 100644 docs/python/source/refs.bib delete mode 100644 docs/python/source/rst/bibliography.rst rename docs/{cpp => }/references.bib (100%) create mode 100644 svzerodplus.pyi delete mode 100644 svzerodsolver/__init__.py delete mode 100644 svzerodsolver/algebra.py delete mode 100644 svzerodsolver/model/__init__.py delete mode 100644 svzerodsolver/model/block.py delete mode 100644 svzerodsolver/model/bloodvessel.py delete mode 100644 svzerodsolver/model/bloodvesseljunction.py delete mode 100644 svzerodsolver/model/dofhandler.py delete mode 100644 svzerodsolver/model/flowreferencebc.py delete mode 100644 svzerodsolver/model/internaljunction.py delete mode 100644 svzerodsolver/model/node.py delete mode 100644 svzerodsolver/model/openloopcoronarybc.py delete mode 100644 svzerodsolver/model/pressurereferencebc.py delete mode 100644 svzerodsolver/model/resistancebc.py delete mode 100644 svzerodsolver/model/windkesselbc.py delete mode 100644 svzerodsolver/runner.py delete mode 100644 svzerodsolver/runnercpp.py delete mode 100644 svzerodsolver/utils.py delete mode 100644 tests/conftest.py delete mode 100644 tests/test_integration_cpp.py diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 130544f4b..f8bb537ca 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -1,4 +1,4 @@ -# This workflow builds and deploys the html documentation for svZeroDSolver. +# This workflow builds and deploys the html documentation for svZeroDPlus. name: Documentation on: [push, pull_request] permissions: @@ -8,32 +8,23 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v2 - - name: Make html directory - run: mkdir docs/html + - name: Make build directory + run: mkdir docs/build - name: Build doxygen documentation continue-on-error: true uses: mattnotmitt/doxygen-action@v1.9.5 with: working-directory: '.' - doxyfile-path: 'docs/cpp/Doxyfile' + doxyfile-path: 'docs/Doxyfile' enable-latex: true - - name: Install sphinx dependencies - run: | - conda create -n zerod python=3.9 - conda run -n zerod pip install -e .[dev] - - name: Build sphinx documentation - continue-on-error: true - run: | - conda run -n zerod sphinx-apidoc -o docs/python/source svzerodsolver - conda run -n zerod sphinx-build -b html -d docs/python/build/doctrees docs/python/source docs/html/python - name: Save documentation uses: actions/upload-artifact@v3 with: name: documentation - path: ./docs/html + 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/html + publish_dir: ./docs/build/html diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 3e7656422..2598e078b 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -1,10 +1,9 @@ -# This workflow builds and tests svZeroDSolver. The Python version is tested on -# ubuntu and the C++ version is built and tested on different versions of -# ubuntu and macOS. +# This workflow builds and tests svZeroDPlus. It is built and tested on +# different versions of ubuntu and macOS. name: Build and test on: [push, pull_request] jobs: - cpp: + test: runs-on: ${{ matrix.os }} strategy: matrix: @@ -18,11 +17,11 @@ jobs: - name: Install svZeroDPlus run: | conda create -n zerod python=3.9 - conda run -n zerod pip install pytest pytest-cov pytest-mock - conda run -n zerod pip install -e . + conda run -n zerod pip install -e ".[dev]" - name: Test the build run: | - conda run -n zerod pytest tests/test_integration_cpp.py + cd tests + conda run -n zerod pytest - name: Build using CMake run: | mkdir Release @@ -39,26 +38,4 @@ jobs: cd test_01 ./svZeroD_interface_test01 ../../../../Release ../../test_01/svzerod_3Dcoupling.json cd ../test_02 - ./svZeroD_interface_test02 ../../../../Release ../../test_02/svzerod_tuned.json - python: - runs-on: ubuntu-latest - steps: - - uses: actions/checkout@v2 - - name: Set up environment - uses: actions/setup-python@v2 - with: - python-version: 3.9 - - name: Install dependencies - run: | - conda create -n zerod python=3.9 - conda run -n zerod pip install pytest pytest-cov pytest-mock - conda run -n zerod pip install -e . - - name: Run pytest - run: | - conda run -n zerod pytest tests/test_integration.py - - name: Save coverage report - uses: actions/upload-artifact@v3 - with: - name: coverage_report - path: htmlcov - + ./svZeroD_interface_test02 ../../../../Release ../../test_02/svzerod_tuned.json \ No newline at end of file diff --git a/.gitignore b/.gitignore index 1d1207475..1f2ecab7f 100644 --- a/.gitignore +++ b/.gitignore @@ -18,11 +18,7 @@ __pycache__ public # Documentation -docs/html -docs/python/source/*.rst -!docs/python/source/index.rst -!docs/python/source/modules.rst -docs/python/build +docs/build # C++ build and externals directories Release*/ diff --git a/README.md b/README.md index 4df48b519..3ebbaac52 100755 --- a/README.md +++ b/README.md @@ -9,14 +9,11 @@ svZeroDPlus is a fast simulation tool for modeling the hemodynamics of -vascular networks using zero-dimensional (0D) lumped parameter models. There -are two different implementations of the solver, a high performance C++ -version and an easy to use Python version. You can find more information about -both under the following links: +vascular networks using zero-dimensional (0D) lumped parameter models. +You can find more information about svZeroDPlus under the following links: -* [**C++ Documentation**](https://StanfordCBCL.github.io/svZeroDPlus/cpp) -* [**C++ Developer Guide**](https://StanfordCBCL.github.io/svZeroDPlus/cpp/md_developer_guide.html) -* [**Python Documentation**](https://StanfordCBCL.github.io/svZeroDPlus/python) +* [**Documentation**](https://StanfordCBCL.github.io/svZeroDPlus) +* [**Developer Guide**](https://StanfordCBCL.github.io/svZeroDPlus/md_developer_guide.html) * [**Wiki**](https://github.com/StanfordCBCL/svZeroDPlus/wiki) * [**Bug Reports**](https://github.com/StanfordCBCL/svZeroDPlus/issues) * [**Forum**](https://github.com/StanfordCBCL/svZeroDPlus/discussions) diff --git a/applications/svzerodplus.cpp b/applications/svzerodplus.cpp index 4450ce2ab..682e1bcbf 100644 --- a/applications/svzerodplus.cpp +++ b/applications/svzerodplus.cpp @@ -44,7 +44,6 @@ namespace py = pybind11; PYBIND11_MODULE(svzerodplus, m) { using Solver = SOLVE::Solver; - m.doc() = "svZeroDSolver"; py::class_(m, "Solver") .def(py::init([](py::dict& config) { const nlohmann::json& config_json = config; @@ -55,17 +54,32 @@ PYBIND11_MODULE(svzerodplus, m) { const auto& config_json = nlohmann::json::parse(ifs); return Solver(config_json); })) - .def("copy", [](Solver& solver) { return Solver(solver); }) .def("run", &Solver::run) .def("get_single_result", &Solver::get_single_result) .def("get_single_result_avg", &Solver::get_single_result_avg) - .def("update_block_params", &Solver::update_block_params); + .def("get_full_result", [](Solver& solver) { + py::module_ pd = py::module_::import("pandas"); + py::module_ io = py::module_::import("io"); + auto result = solver.get_full_result(); + return pd.attr("read_csv")(io.attr("StringIO")(result)); + }); m.def("simulate", [](py::dict& config) { + py::module_ pd = py::module_::import("pandas"); + py::module_ io = py::module_::import("io"); const nlohmann::json& config_json = config; auto solver = Solver(config_json); solver.run(); - return solver.get_full_result(); + return pd.attr("read_csv")(io.attr("StringIO")(solver.get_full_result())); + }); + m.def("simulate", [](std::string config_file) { + py::module_ pd = py::module_::import("pandas"); + py::module_ io = py::module_::import("io"); + std::ifstream ifs(config_file); + const auto& config_json = nlohmann::json::parse(ifs); + auto solver = Solver(config_json); + solver.run(); + return pd.attr("read_csv")(io.attr("StringIO")(solver.get_full_result())); }); m.def("calibrate", [](py::dict& config) { const nlohmann::json& config_json = config; diff --git a/docs/cpp/Doxyfile b/docs/Doxyfile similarity index 96% rename from docs/cpp/Doxyfile rename to docs/Doxyfile index 7ab866dfc..cf1fd6669 100644 --- a/docs/cpp/Doxyfile +++ b/docs/Doxyfile @@ -8,7 +8,7 @@ PROJECT_NAME = svZeroDPlus PROJECT_NUMBER = PROJECT_BRIEF = PROJECT_LOGO = -OUTPUT_DIRECTORY = docs/html +OUTPUT_DIRECTORY = docs/build CREATE_SUBDIRS = NO OUTPUT_LANGUAGE = English BRIEF_MEMBER_DESC = YES @@ -27,7 +27,7 @@ ABBREVIATE_BRIEF = "The $name class" \ ALWAYS_DETAILED_SEC = YES INLINE_INHERITED_MEMB = NO FULL_PATH_NAMES = YES -STRIP_FROM_PATH = "src" "docs/cpp/pages" +STRIP_FROM_PATH = "src" "docs/pages" STRIP_FROM_INC_PATH = SHORT_NAMES = NO JAVADOC_AUTOBRIEF = NO @@ -90,8 +90,8 @@ SHOW_USED_FILES = YES SHOW_FILES = YES SHOW_NAMESPACES = YES FILE_VERSION_FILTER = -LAYOUT_FILE = ./docs/cpp/DoxygenLayout.xml -CITE_BIB_FILES = ./docs/cpp/references.bib +LAYOUT_FILE = ./docs/DoxygenLayout.xml +CITE_BIB_FILES = ./docs/references.bib #--------------------------------------------------------------------------- # Configuration options related to warning and progress messages #--------------------------------------------------------------------------- @@ -106,7 +106,7 @@ WARN_LOGFILE = # Configuration options related to the input files #--------------------------------------------------------------------------- INPUT = "src" \ - "docs/cpp/pages" + "docs/pages" INPUT_ENCODING = UTF-8 FILE_PATTERNS = *.cpp \ *.hpp \ @@ -122,7 +122,7 @@ EXCLUDE_SYMBOLS = EXAMPLE_PATH = EXAMPLE_PATTERNS = * EXAMPLE_RECURSIVE = NO -IMAGE_PATH = docs/cpp +IMAGE_PATH = images INPUT_FILTER = FILTER_PATTERNS = FILTER_SOURCE_FILES = NO @@ -149,12 +149,12 @@ IGNORE_PREFIX = # Configuration options related to the HTML output #--------------------------------------------------------------------------- GENERATE_HTML = YES -HTML_OUTPUT = cpp +HTML_OUTPUT = html HTML_FILE_EXTENSION = .html HTML_HEADER = HTML_FOOTER = HTML_STYLESHEET = -HTML_EXTRA_STYLESHEET = ./docs/cpp/doxygen-awesome.css +HTML_EXTRA_STYLESHEET = ./docs/doxygen-awesome.css HTML_EXTRA_FILES = HTML_COLORSTYLE_HUE = 220 HTML_COLORSTYLE_SAT = 100 diff --git a/docs/cpp/DoxygenLayout.xml b/docs/DoxygenLayout.xml similarity index 100% rename from docs/cpp/DoxygenLayout.xml rename to docs/DoxygenLayout.xml diff --git a/docs/cpp/doxygen-awesome.css b/docs/doxygen-awesome.css similarity index 100% rename from docs/cpp/doxygen-awesome.css rename to docs/doxygen-awesome.css diff --git a/docs/cpp/pages/developer_guide.md b/docs/pages/developer_guide.md similarity index 99% rename from docs/cpp/pages/developer_guide.md rename to docs/pages/developer_guide.md index 8f1f6e6d0..876c75b82 100644 --- a/docs/cpp/pages/developer_guide.md +++ b/docs/pages/developer_guide.md @@ -95,7 +95,7 @@ The documentation is automatically build in the GitHub CI/CD and published on GitHub pages. If you want to build the documentation locally, you can use: ``` -doxygen docs/cpp/Doxyfile +doxygen docs/Doxyfile ``` If you do not have Doxygen install you can do that with `brew install doxygen` diff --git a/docs/cpp/pages/main.md b/docs/pages/main.md similarity index 84% rename from docs/cpp/pages/main.md rename to docs/pages/main.md index 1d7138948..21e916bfe 100644 --- a/docs/cpp/pages/main.md +++ b/docs/pages/main.md @@ -125,35 +125,35 @@ you installed svZerodPlus via pip to enable this feature. We start by importing svzerodplus: ```python -import svzerodplus +>>> import svzerodplus ``` Next, we create a solver from our configuration. The configuration can be specified by either a path to a JSON file: ```python -solver = svzerodplus.Solver("tests/cases/steadyFlow_RLC_R.json") +>>> solver = svzerodplus.Solver("tests/cases/steadyFlow_RLC_R.json") ``` or as a Python dictionary: ```python -my_config = {...} -solver = svzerodplus.Solver(my_config) +>>> my_config = {...} +>>> solver = svzerodplus.Solver(my_config) ``` To run the simulation we add: ```python -solver.run() +>>> solver.run() ``` The simulation result is now saved in the solver instance. We can obtain results for individual degrees-of-freedom (DOFs) as ```python -solver.get_single_result("flow:INFLOW:branch0_seg0") +>>> solver.get_single_result("flow:INFLOW:branch0_seg0") ->>> array([5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., +array([5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., @@ -164,11 +164,57 @@ The naming of the DOFs is similar to how results are written if the simulation option `output_variable_based` is activated (see below). We can also obtain the mean result for a DOF over time with: ```python -solver.get_single_result_avg("flow:INFLOW:branch0_seg0") +>>> solver.get_single_result_avg("flow:INFLOW:branch0_seg0") + +5.0 +``` + +Or the result of the full simulation as a pandas data frame: + +```python +>>> solver.get_full_result() + + name time flow_in flow_out pressure_in pressure_out +0 branch0_seg0 0.00 5.0 5.0 1100.0 600.0 +1 branch0_seg0 0.01 5.0 5.0 1100.0 600.0 +2 branch0_seg0 0.02 5.0 5.0 1100.0 600.0 +3 branch0_seg0 0.03 5.0 5.0 1100.0 600.0 +4 branch0_seg0 0.04 5.0 5.0 1100.0 600.0 +.. ... ... ... ... ... ... +96 branch0_seg0 0.96 5.0 5.0 1100.0 600.0 +97 branch0_seg0 0.97 5.0 5.0 1100.0 600.0 +98 branch0_seg0 0.98 5.0 5.0 1100.0 600.0 +99 branch0_seg0 0.99 5.0 5.0 1100.0 600.0 +100 branch0_seg0 1.00 5.0 5.0 1100.0 600.0 + +[101 rows x 6 columns] +``` + +There is also a function to retrieve the full result directly based on a given configuration: + +```python + +>>> my_config = {...} +>>> svzerodplus.simulate(my_config) + + name time flow_in flow_out pressure_in pressure_out +0 branch0_seg0 0.00 5.0 5.0 1100.0 600.0 +1 branch0_seg0 0.01 5.0 5.0 1100.0 600.0 +2 branch0_seg0 0.02 5.0 5.0 1100.0 600.0 +3 branch0_seg0 0.03 5.0 5.0 1100.0 600.0 +4 branch0_seg0 0.04 5.0 5.0 1100.0 600.0 +.. ... ... ... ... ... ... +96 branch0_seg0 0.96 5.0 5.0 1100.0 600.0 +97 branch0_seg0 0.97 5.0 5.0 1100.0 600.0 +98 branch0_seg0 0.98 5.0 5.0 1100.0 600.0 +99 branch0_seg0 0.99 5.0 5.0 1100.0 600.0 +100 branch0_seg0 1.00 5.0 5.0 1100.0 600.0 + +[101 rows x 6 columns] ->>> 5.0 ``` + ## Configuration svZeroDSolver is configured using either a JSON file or a Python diff --git a/docs/python/source/conf.py b/docs/python/source/conf.py deleted file mode 100644 index 718e6a951..000000000 --- a/docs/python/source/conf.py +++ /dev/null @@ -1,41 +0,0 @@ -"""Configuration file for the Sphinx documentation builder.""" -import svzerodsolver - -project = svzerodsolver.NAME -release = svzerodsolver.VERSION -copyright = svzerodsolver.COPYRIGHT - -extensions = [ - "sphinx.ext.autodoc", - "sphinx.ext.doctest", - "sphinx.ext.todo", - "sphinx.ext.coverage", - "sphinx.ext.mathjax", - "sphinx.ext.viewcode", - "sphinx.ext.inheritance_diagram", - "sphinx.ext.napoleon", - "m2r2", - "sphinxcontrib.bibtex", - "sphinx_autodoc_typehints", -] - -source_suffix = [".rst", ".md"] -master_doc = "index" -exclude_patterns = [] -html_theme = "pydata_sphinx_theme" -# html_logo = "img/logo.png" -# html_favicon = "img/favicon.png" -bibtex_bibfiles = ["refs.bib"] -bibtex_default_style = "unsrt" -set_type_checking_flag = False - -html_theme_options = { - "logo": { - "text": "svZeroDPlus Python Documentation", - }, - "collapse_navigation": True, - "show_prev_next": False, - "github_url": "https://github.com/StanfordCBCL/svZeroDSolver", - "navbar_end": ["navbar-icon-links"], -} -html_title = "%s %s Manual" % (project, release) diff --git a/docs/python/source/index.rst b/docs/python/source/index.rst deleted file mode 100644 index a7e212e5e..000000000 --- a/docs/python/source/index.rst +++ /dev/null @@ -1,98 +0,0 @@ -.. _mainpage: - -********************************** -svZeroDSolver Python Documentation -********************************** - -**Version**: |release| - -**Useful links**: -`Source Repository `_ | -`Issue Tracker `_ | -`SimVascular `_ - -svZeroDSolver is a Python package that simulates the hemodynamics in zero-dimensional -(0D) lumped parameter models of vascular networks. These 0D models are governed -by differential algebraic equations (DAEs). - -The solver uses a highly modular framework to model the vascular anatomy, using -individual 0D elements to represent different parts of the vascular anatomy -(and boundary conditions). The individual 0D elements and their associated -governing equations defined in `models` subpackage. -In the `algebra` module, the blocks are assembled and simulated using the -generalized-alpha time-stepping method. - -Installation -************ - -svZeroDSolver and all its dependencies can be installed easily via pip. - -.. code-block:: bash - - pip install git+https://github.com/StanfordCBCL/svZeroDSolver.git - - -For Contributers -================ - -The following guide provides all necessary steps to install your local -svZeroDSolver repository via pip in editable mode to allow for local code changes -to reflect in the package. - -If you are contributing to svZeroDSolver, it is highly recommended to use a virtual -environment like [Miniconda](https://docs.conda.io/en/latest/miniconda.html). -After installing Miniconda you can create a new environment and enter it using: - -.. code-block:: bash - - conda create -n zerodsolver python=3.9 - conda activate zerodsolver - -After that, enter the repository folder and install the svZeroDSolver -**with development related dependencies** using: - -.. code-block:: bash - - pip install -e .[dev] - - -*If you are using the `zsh` shell, enter: `pip install -e ".[dev]"`* - -Usage -***** - -Command line -============ - -To run svZeroDSolver form the command line, run: - -.. code-block:: bash - - zerod SOLVER_INPUT_FILE SOLVER_OUTPUT_FILE - - -For more information about command line options, enter: - -.. code-block:: bash - - zerod --help - -As a python module -================== - -.. code-block:: python - - import svzerodsolver - svzerodsolver.runner.run_from_file('input.json', 'output.json') - -This variant enables running svZeroDSolver within a user-defined Python code -(e.g. parameter optimization, uncertainty quantification) - -More information -**************** - -.. toctree:: - :maxdepth: 1 - - Documentation - References diff --git a/docs/python/source/modules.rst b/docs/python/source/modules.rst deleted file mode 100644 index 81c812cae..000000000 --- a/docs/python/source/modules.rst +++ /dev/null @@ -1,7 +0,0 @@ -Documentation -============= - -.. toctree:: - :maxdepth: 3 - - svzerodsolver diff --git a/docs/python/source/refs.bib b/docs/python/source/refs.bib deleted file mode 100644 index e12f64014..000000000 --- a/docs/python/source/refs.bib +++ /dev/null @@ -1,23 +0,0 @@ -@article{JANSEN2000305, -title = {A generalized-α method for integrating the filtered Navier–Stokes equations with a stabilized finite element method}, -journal = {Computer Methods in Applied Mechanics and Engineering}, -volume = {190}, -number = {3}, -pages = {305-319}, -year = {2000}, -issn = {0045-7825}, -doi = {https://doi.org/10.1016/S0045-7825(00)00203-6}, -url = {https://www.sciencedirect.com/science/article/pii/S0045782500002036}, -author = {Kenneth E. Jansen and Christian H. Whiting and Gregory M. Hulbert}, -abstract = {A generalized-α method is developed and analyzed for linear, first-order systems. The method is then extended to the filtered Navier–Stokes equations within the context of a stabilized finite element method. The formulation is studied through the application to laminar flow past a circular cylinder and turbulent flow past a long, transverse groove. The method is formulated to obtain a second-order accurate family of time integrators whose high frequency amplification factor is the sole free parameter. Such an approach allows the replication of midpoint rule (zero damping), Gear's method (maximal damping), or anything in between.} -} -@article{kim_coronary, -author = {Kim, H and Vignon-Clementel, Irene and Coogan, Jessica and Figueroa, Carlos and Jansen, Kenneth and Taylor, CA}, -year = {2010}, -month = {10}, -pages = {3195-209}, -title = {Patient-Specific Modeling of Blood Flow and Pressure in Human Coronary Arteries}, -volume = {38}, -journal = {Annals of biomedical engineering}, -doi = {10.1007/s10439-010-0083-6} -} \ No newline at end of file diff --git a/docs/python/source/rst/bibliography.rst b/docs/python/source/rst/bibliography.rst deleted file mode 100644 index 07fc32015..000000000 --- a/docs/python/source/rst/bibliography.rst +++ /dev/null @@ -1,7 +0,0 @@ -.. _references: - -********** -References -********** - -.. bibliography:: \ No newline at end of file diff --git a/docs/cpp/references.bib b/docs/references.bib similarity index 100% rename from docs/cpp/references.bib rename to docs/references.bib diff --git a/pyproject.toml b/pyproject.toml index 40b5fef85..7fa0a2bfd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,9 +4,6 @@ build-backend = "setuptools.build_meta" [tool.pytest.ini_options] addopts = """ - --cov-report=html - --cov-report=term - --cov=svzerodsolver --ignore=Release --ignore=Debug """ diff --git a/setup.cfg b/setup.cfg index 46f92998e..b458f3f89 100644 --- a/setup.cfg +++ b/setup.cfg @@ -9,28 +9,14 @@ long_description_content_type = text/markdown packages = find: python_requires = >=3.0 install_requires = - numpy - scipy - orjson - click - orjson # Faster than json - pyarrow # For faster csv reading pandas + numpy [options.extras_require] dev = pytest - pytest-cov - pytest-mock - sphinx - sphinx-autodoc-typehints - m2r2 - pydata-sphinx-theme - sphinxcontrib-bibtex - sphinxcontrib-mermaid [options.entry_points] console_scripts = svzerodsolver = svzerodplus:run_simulation_cli svzerodcalibrator = svzerodplus:run_calibration_cli - diff --git a/src/interface/interface.cpp b/src/interface/interface.cpp index 29b77b311..55928632a 100644 --- a/src/interface/interface.cpp +++ b/src/interface/interface.cpp @@ -35,7 +35,7 @@ #include -//#include "io/configreader.hpp" +// #include "io/configreader.hpp" typedef double T; diff --git a/src/io/csvwriter.hpp b/src/io/csvwriter.hpp index 4c9e017cb..244be40ff 100644 --- a/src/io/csvwriter.hpp +++ b/src/io/csvwriter.hpp @@ -90,7 +90,7 @@ std::string to_vessel_csv(std::vector ×, continue; } - std::string name = block->get_name(); + std::string name = model.get_block_name(i); inflow_dof = block->inlet_nodes[0]->flow_dof; outflow_dof = block->outlet_nodes[0]->flow_dof; inpres_dof = block->inlet_nodes[0]->pres_dof; diff --git a/svzerodplus.pyi b/svzerodplus.pyi new file mode 100644 index 000000000..5ed5dce6c --- /dev/null +++ b/svzerodplus.pyi @@ -0,0 +1,95 @@ +# This file contains the stubs for the public Python interface of svZeroDPlus. +# This is necessary for supporting auto-completion in Visual Studio Code or +# PyCharm. +"""svZeroDPlus Python interface.""" +from __future__ import annotations +import numpy +import typing +import pandas + +__all__ = ["Solver", "calibrate", "simulate"] + +class Solver: + """Lumped-parameter solver.""" + + @typing.overload + def __init__(self, arg0: dict) -> None: + """Create a new lumped-parameter solver. + + Args: + arg0: Solver configuration dictionary. + """ + ... + @typing.overload + def __init__(self, arg0: str) -> None: + """Create a new lumped-parameter solver. + + Args: + arg0: Path to solver configuration file. + """ + ... + def get_full_result(self) -> pandas.DataFrame: + """Get the full result of the simulation. + + Returns: + Simulation result as a dataframe. + """ + ... + def get_single_result(self, arg0: str) -> numpy.ndarray: + """Get the simulation result for a single degree-of-freedom (DOF). + + Args: + arg0: Name of the DOF. + + Returns: + Time-dependent simulation result for DOF. + """ + ... + def get_single_result_avg(self, arg0: str) -> float: + """Get the mean simulation result for a single degree-of-freedom (DOF). + + Args: + arg0: Name of the DOF. + + Returns: + Mean simulation result for the DOF. + """ + ... + def run(self) -> None: + """Run the simulation.""" + ... + +def calibrate(arg0: dict) -> dict: + """Run a Levenberg-Marquardt calibration. + + Args: + arg0: Calibration configuration dictionary. + + Returns: + Calibrated 0D solver input file. + """ + ... + +@typing.overload +def simulate(arg0: dict) -> pandas.DataFrame: + """Run a lumped-parameter simulation. + + Args: + arg0: Simulation configuration file. + + Returns: + Simulation result as a dataframe. + """ + ... + +@typing.overload +def simulate(arg0: str) -> pandas.DataFrame: + """Run a lumped-parameter simulation. + + Args: + arg0: Path to solver configuration file. + + Returns: + Simulation result as a dataframe. + """ + ... diff --git a/svzerodsolver/__init__.py b/svzerodsolver/__init__.py deleted file mode 100644 index f79b7a425..000000000 --- a/svzerodsolver/__init__.py +++ /dev/null @@ -1,34 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - - -NAME = "svZeroDPlus" -VERSION = "v1.0" -COPYRIGHT = "Stanford University, The Regents of the University of California, and others." diff --git a/svzerodsolver/algebra.py b/svzerodsolver/algebra.py deleted file mode 100644 index 904fe3333..000000000 --- a/svzerodsolver/algebra.py +++ /dev/null @@ -1,210 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""This module holds all algebra related classes and methods.""" -import numpy as np -from scipy.sparse import linalg - - -class GeneralizedAlpha: - """Generalized alpha time integrator. :cite:`JANSEN2000305` - - Solves system E*ydot + F*y + C = 0 with generalized alpha and Newton-Raphson - for non-linear residual. - - Args: - rho: Generalized-alpha parameter rho. - n: System size. - time_step_size: Time step size. - atol: Absolute tolerance. - max_iter: Maximum nonlinear iterations. - """ - - def __init__(self, rho, n, time_step_size, atol=1e-8, max_iter=30): - - # Setup constants for generalized alpha time integration - self.alpha_m = 0.5 * (3.0 - rho) / (1.0 + rho) - self.alpha_f = 1.0 / (1.0 + rho) - self.gamma = 0.5 + self.alpha_m - self.alpha_f - - self.fac = self.alpha_m / (self.alpha_f * self.gamma) - - # problem dimension - self.n = n - self.time_step_size = time_step_size - self.time_step_size_inv = 1.0 / time_step_size - - self.atol = atol - self.max_iter = max_iter - - # jacobian matrix - if self.n > 800: - self.solve = linalg.spsolve - else: - self.solve = np.linalg.solve - - # stores matrices E, F, vector C, and tangent matrices dE, dF, dC - self.mat = {} - self.mats = ["E", "F", "dE", "dF", "dC"] - self.vecs = ["C"] - for m in self.mats: - self.mat[m] = np.zeros((self.n, self.n)) - for v in self.vecs: - self.mat[v] = np.zeros(self.n) - - def assemble(self, block_list): - """ - Assemble block matrices into global matrices - """ - for bl in block_list: - bl.assemble(self.mat) - - def step(self, y, ydot, time, block_list): - """ - Perform one time step - """ - - # Make copies to prevent overwriting - y = y.copy() - ydot = ydot.copy() - - # initial guess for time step - curr_y = y + 0.5 * self.time_step_size * ydot - curr_ydot = ydot * ((self.gamma - 0.5) / self.gamma) - - # Substep level quantities - yaf = y + self.alpha_f * (curr_y - y) - ydotam = ydot + self.alpha_m * (curr_ydot - ydot) - - # initialize solution - time = time + self.alpha_f * self.time_step_size - - # initialize blocks - for b in block_list: - b.update_time(time) - - fac_ydotam = self.fac * self.time_step_size_inv - for iter in range(self.max_iter): - # update solution-dependent blocks - for b in block_list: - b.update_solution(yaf) - - # Assemble - self.assemble(block_list) - res = ( - -self.mat["E"].dot(ydotam) - - self.mat["F"].dot(yaf) - - self.mat["C"] - ) - - # Check termination criteria - if np.abs(res).max() <= self.atol: - break - - lhs = self.mat["F"] + ( - self.mat["dE"] - + self.mat["dF"] - + self.mat["dC"] - + self.mat["E"] * self.fac * self.time_step_size_inv - ) - - # solve for Newton increment - dy = self.solve(lhs, res) - - # update solution - yaf += dy - ydotam += dy * fac_ydotam - - if iter == self.max_iter - 1: - print( - f"Max NR iterations reached at time: {time:.3f}s with, max error: {max(abs(res))}" - ) - - # update time step - curr_y = y + (yaf - y) / self.alpha_f - curr_ydot = ydot + (ydotam - ydot) / self.alpha_m - - return curr_y, curr_ydot - - -def run_integrator( - block_list, - dofhandler, - num_time_steps, - time_step_size, - y_initial=None, - ydot_initial=None, - rho=0.1, - method="genalpha", - atol=10 - 8, - max_iter=30, -): - """Run time integration. - - Args: - block_list: List of model blocks. - dofhandler: Degree-of-freedom handler. - num_time_steps: Number of time steps. - time_step_size: Time step size. - y_initial: Initial y. - ydot_initial: Initial ydot. - rho: Generalized-alpha parameter rho. - method: Time integration method to use. - atol: Absolute tolerance. - max_iter: Maximum nonlinear iterations. - """ - - y = np.zeros(dofhandler.n) if y_initial is None else y_initial.copy() - ydot = ( - np.zeros(dofhandler.n) if ydot_initial is None else ydot_initial.copy() - ) - - y_over_time = [y.copy()] - ydot_over_time = [ydot.copy()] - time_steps = np.arange( - 0.0, num_time_steps * time_step_size, time_step_size, dtype=float - ) - if method == "genalpha": - integrator = GeneralizedAlpha( - rho, y.shape[0], time_step_size, atol=atol, max_iter=max_iter - ) - else: - raise ValueError(f"Unknown integration method {method}.") - - for time in time_steps[:-1]: - y, ydot = integrator.step( - y, - ydot, - time, - block_list, - ) - y_over_time.append(y.copy()) - ydot_over_time.append(ydot.copy()) - - return time_steps, y_over_time, ydot_over_time diff --git a/svzerodsolver/model/__init__.py b/svzerodsolver/model/__init__.py deleted file mode 100644 index c128e39dd..000000000 --- a/svzerodsolver/model/__init__.py +++ /dev/null @@ -1,29 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/svzerodsolver/model/block.py b/svzerodsolver/model/block.py deleted file mode 100644 index e4d174f3b..000000000 --- a/svzerodsolver/model/block.py +++ /dev/null @@ -1,164 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""This module holds the Block class.""" -from __future__ import annotations -from abc import ABC, abstractclassmethod - -import numpy as np -from scipy.interpolate import CubicSpline - -from .dofhandler import DOFHandler - - -class Block(ABC): - """Base class for 0D model components. - - A block is the base class of 0D model elements. It is the place where the - contribution of an element to the global system is controlled. It defines - all relevant attributes and methods of an element and a few common helpers - for setting it up. - - Attributes: - name: Name of the block. - inflow_nodes: Inflow nodes of the element. - outflow_nodes: Outflow nodes of the element. - """ - - # Number of equations the element contirbutes at assembly - _NUM_EQUATIONS = None - - # Number of internal variables of the block - _NUM_INTERNAL_VARS = 0 - - def __init__(self, params: dict = None, name: str = None): - """Create a new Block. - - Args: - params: The configuration paramaters of the block. Mostly comprised - of constants for element contribution calculation. - name: Optional name of the block. - """ - self.name = name - self._params = params - - # Inflow and outflow wires of the block. Will be set by wires - self.inflow_nodes = [] - self.outflow_nodes = [] - - # Element contribution matrices - self._mat = {} - self._vec = {} - - # row and column indices of block in global matrix - self._global_row_id = None - self._flat_row_ids = None - self._flat_col_ids = None - - @abstractclassmethod - def from_config(cls, config: dict) -> "Block": - """Create block from config dictionary. - - Args: - config: The configuration dict for the block. - - Returns: - The block instance. - """ - pass - - def setup_dofs(self, dofhandler: DOFHandler) -> None: - """Setup degree of freedoms of the block. - - Registers equations and internal variables at a DOF handler. - - Args: - dofhandler: The DOF handler to register the variables and equations - at. - """ - - # Register internal variables - internal_vars = [ - dofhandler.register_variable(f"var_{i}_{self.name}") - for i in range(self._NUM_INTERNAL_VARS) - ] - - # Collect assembly column ids based from inflow/outflow wires and - # internal variables - global_col_id = [] - for wire in self.inflow_nodes: - global_col_id += [wire.pres_dof, wire.flow_dof] - for wire in self.outflow_nodes: - global_col_id += [wire.pres_dof, wire.flow_dof] - global_col_id += internal_vars - - # Register equations of the block - self._global_row_id = [ - dofhandler.register_equation() for _ in range(self._NUM_EQUATIONS) - ] - - # Create flat indices to assemble matrices as flattend array (faster) - meshgrid = np.array( - np.meshgrid(self._global_row_id, global_col_id) - ).T.reshape(-1, 2) - self._flat_row_ids, self._flat_col_ids = meshgrid[:, 0], meshgrid[:, 1] - - def assemble(self, mat: dict[str, np.ndarray]) -> None: - """Assemble block to global system. - - Args: - mat: Global system. - """ - for key, value in self._vec.items(): - mat[key][self._global_row_id] = value - for key, value in self._mat.items(): - mat[key][self._flat_row_ids, self._flat_col_ids] = value.ravel() - - def update_time(self, time: float) -> None: - """Update time dependent element contributions. - - Args: - time: Current time. - """ - pass - - def update_solution(self, y: np.ndarray) -> None: - """Update solution dependent element contributions. - - Args: - y: Current solution. - """ - pass - - def _interpolate(self, times, values): - if times is None: - raise ValueError("No time sequence provided for interpolation.") - return CubicSpline( - np.array(times), np.array(values), bc_type="periodic" - ) diff --git a/svzerodsolver/model/bloodvessel.py b/svzerodsolver/model/bloodvessel.py deleted file mode 100644 index 44890f05c..000000000 --- a/svzerodsolver/model/bloodvessel.py +++ /dev/null @@ -1,110 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""This module holds the BloodVessel class.""" -import numpy as np - -from .block import Block - - -class BloodVessel(Block): - r"""Resistor-capacitor-inductor blood vessel with optional stenosis. - - Models the mechanical behavior of a bloodvessel with optional stenosis. - - Attributes: - name: Name of the block. - inflow_nodes: Inflow nodes of the element. - outflow_nodes: Outflow nodes of the element. - """ - - _NUM_EQUATIONS = 3 - _NUM_INTERNAL_VARS = 1 - - def __init__(self, params: dict = None, name: str = None): - """Create a new BloodVessel instance. - - Args: - params: The configuration paramaters of the block. Mostly comprised - of constants for element contribution calculation. - name: Optional name of the block. - """ - super().__init__(params=params, name=name) - - # the ordering of the solution variables is : (P_in, Q_in, P_out, Q_out) - self._mat["E"] = np.zeros((3, 5), dtype=float) - self._mat["E"][0, 3] = -self._params["L"] - self._mat["E"][1, 4] = -self._params["C"] - self._mat["F"] = np.array( - [ - [1.0, -self._params["R"], -1.0, 0.0, 0.0], - [0.0, 1.0, 0.0, -1.0, 0.0], - [1.0, -self._params["R"], 0.0, 0.0, -1.0], - ], - dtype=float, - ) - self._mat["dF"] = np.zeros((3, 5), dtype=float) - - # This class's update_solution is only needed if class is configured with - # non-zero stenosis coefficient. If stenosis coefficient is 0.0 the - # emtpy parent update_solution can be used and is faster. - if self._params["stenosis_coefficient"] == 0.0: - self.update_solution = super().update_solution - - @classmethod - def from_config(cls, config: dict) -> "BloodVessel": - """Create block from config dictionary. - - Args: - config: The configuration dict for the block. - - Returns: - The block instance. - """ - params = dict( - R=config["zero_d_element_values"].get("R_poiseuille"), - C=config["zero_d_element_values"].get("C", 0.0), - L=config["zero_d_element_values"].get("L", 0.0), - stenosis_coefficient=config["zero_d_element_values"].get( - "stenosis_coefficient", 0.0 - ), - ) - return cls(params=params, name="V" + str(config["vessel_id"])) - - def update_solution(self, y: np.ndarray) -> None: - """Update solution dependent element contributions. - - Args: - y: Current solution. - """ - Q_in = np.abs(y[self.inflow_nodes[0].flow_dof]) - fac1 = -self._params["stenosis_coefficient"] * Q_in - fac2 = fac1 - self._params["R"] - self._mat["F"][[0, 2], 1] = fac2 - self._mat["dF"][[0, 2], 1] = fac1 diff --git a/svzerodsolver/model/bloodvesseljunction.py b/svzerodsolver/model/bloodvesseljunction.py deleted file mode 100644 index ebfdba7a9..000000000 --- a/svzerodsolver/model/bloodvesseljunction.py +++ /dev/null @@ -1,45 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""This module holds the BloodVesselJunction class.""" -from .internaljunction import InternalJunction - - -class BloodVesselJunction(InternalJunction): - """Blood vessel junction. - - (dummy for future implementation of blood pressure losses at junctions) - - Attributes: - name: Name of the block. - inflow_nodes: Inflow nodes of the element. - outflow_nodes: Outflow nodes of the element. - """ - - pass diff --git a/svzerodsolver/model/dofhandler.py b/svzerodsolver/model/dofhandler.py deleted file mode 100644 index 9d40a0e31..000000000 --- a/svzerodsolver/model/dofhandler.py +++ /dev/null @@ -1,77 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""This module holds the DOFHandler class.""" - - -class DOFHandler: - """Degree-of-freedom handler. - - This class handles degrees-of-freedom for model variables and equations. - It assigns each element with row and column indices which it can use to - assemble it's local contributions into the local contributions into the - global system. - - Attributes: - variables: List of variable names corresponding to the global IDs. - Variables without a name have an entry None. - """ - - def __init__(self): - """Create a new DOFHandler instance.""" - self._var_counter = -1 - self._eqn_counter = -1 - self.variables: list[str] = [] - - @property - def n(self) -> int: - """Size of the system.""" - return self._eqn_counter + 1 - - def register_variable(self, name: str = None) -> int: - """Register a new variable and get global index. - - Args: - name: Optional name of the variables. - - Returns: - global_id: Global id of the variable. - """ - self._var_counter += 1 - self.variables.append(name) - return self._var_counter - - def register_equation(self) -> int: - """Register a new equation and get global index. - - Returns: - global_id: Global id of the equation. - """ - self._eqn_counter += 1 - return self._eqn_counter diff --git a/svzerodsolver/model/flowreferencebc.py b/svzerodsolver/model/flowreferencebc.py deleted file mode 100644 index 729057021..000000000 --- a/svzerodsolver/model/flowreferencebc.py +++ /dev/null @@ -1,94 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""This module holds the FlowReferenceBC class.""" -from typing import Sequence - -import numpy as np - -from .block import Block - - -class FlowReferenceBC(Block): - r"""Flow reference boundary condition. - - Applies a prescribed flow to a boundary. - - Attributes: - name: Name of the block. - inflow_nodes: Inflow nodes of the element. - outflow_nodes: Outflow nodes of the element. - """ - - _NUM_EQUATIONS = 1 - - def __init__(self, params: dict = None, name: str = None): - """Create a new BloodVessel instance. - - Args: - params: The configuration paramaters of the block. Mostly comprised - of constants for element contribution calculation. - name: Optional name of the block. - """ - super().__init__(params=params, name=name) - - if isinstance(self._params["Q"], Sequence): - self._q_func = self._interpolate( - self._params["time"], self._params["Q"] - ) - self._vec["C"] = np.zeros(1, dtype=float) - else: - self._vec["C"] = np.array([-self._params["Q"]], dtype=float) - self.update_time = super().update_time - - self._mat["F"] = np.array([[0.0, 1.0]], dtype=float) - - @classmethod - def from_config(cls, config): - """Create block from config dictionary. - - Args: - config: The configuration dict for the block. - - Returns: - The block instance. - """ - params = dict( - time=config["bc_values"].get("t", None), - Q=config["bc_values"].get("Q"), - ) - return cls(params=params, name=config["name"]) - - def update_time(self, time): - """Update time dependent element contributions. - - Args: - time: Current time. - """ - self._vec["C"][0] = -self._q_func(time) diff --git a/svzerodsolver/model/internaljunction.py b/svzerodsolver/model/internaljunction.py deleted file mode 100644 index 17d31c872..000000000 --- a/svzerodsolver/model/internaljunction.py +++ /dev/null @@ -1,94 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""This module holds the InternalJunction class.""" -import numpy as np - -from .block import Block - - -class InternalJunction(Block): - r"""Junction block. - - Represents a basic model junction without special mechanical behavior. - Across all inlets and outlets of the junction, mass is conserved and - pressure is continuous. - - Attributes: - name: Name of the block. - inflow_nodes: Inflow nodes of the element. - outflow_nodes: Outflow nodes of the element. - """ - - def setup_dofs(self, dofhandler): - """Setup degree of freedoms of the block. - - Registers equations and internal variables at a DOF handler. - - Args: - dofhandler: The DOF handler to register the variables and equations - at. - """ - # Derive number of inlets and outlets - num_inlets = len(self.inflow_nodes) - num_outlets = len(self.outflow_nodes) - - # Set number of equations of a junction block based on number of - # inlets/outlets. Must be set before calling parent constructor - self._NUM_EQUATIONS = num_inlets + num_outlets - super().setup_dofs(dofhandler) - - # Set some constant element element contributions that needed - # _NUM_EQUATIONS - self._mat["F"] = np.zeros( - (self._NUM_EQUATIONS, self._NUM_EQUATIONS * 2) - ) - for i in range(self._NUM_EQUATIONS - 1): - self._mat["F"][i, [0, 2 * i + 2]] = [1.0, -1.0] - self._mat["F"][-1, np.arange(1, 2 * self._NUM_EQUATIONS, 2)] = [ - 1.0 - ] * num_inlets + [-1.0] * num_outlets - - @classmethod - def from_config(cls, config: dict) -> "InternalJunction": - """Create block from config dictionary. - - Args: - config: The configuration dict for the block. - - Returns: - The block instance. - """ - name = config["junction_name"] - if not name.startswith("J") and not name[1].isnumeric(): - raise ValueError( - f"Invalid junction name {name}. Junction names must " - "start with J following by a number." - ) - return cls(name=name) diff --git a/svzerodsolver/model/node.py b/svzerodsolver/model/node.py deleted file mode 100644 index 5c350b4a4..000000000 --- a/svzerodsolver/model/node.py +++ /dev/null @@ -1,73 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""This module holds the Node class.""" -from .block import Block -from .dofhandler import DOFHandler - - -class Node: - """Node. - - Nodes connect two blocks with each other. Each node corresponds to a - flow and pressure value of the system. - - Attributes: - name: Name of the node. - flow_dof: Global ID of the flow value associated with the node. - pres_dof: Global ID of the pressure value associated with the node. - """ - - def __init__(self, ele1: Block, ele2: Block, name: str = None): - """Create a new Node instance. - - Args: - ele1: First element for the node to connect. - ele2: Second element for the node to connect. - name: Optional name of the node. - """ - self.name = name - self.flow_dof: int = None - self.pres_dof: int = None - - # Make the node the ouflow node of the first element and the inflow - # node of the second element - ele1.outflow_nodes.append(self) - ele2.inflow_nodes.append(self) - - def setup_dofs(self, dofhandler: DOFHandler): - """Setup degree of freedoms of the node. - - Registers the pressure and the flow variable at a DOF handler. - - Args: - dofhandler: The DOF handler to register the variables at. - """ - self.flow_dof = dofhandler.register_variable("Q_" + self.name) - self.pres_dof = dofhandler.register_variable("P_" + self.name) diff --git a/svzerodsolver/model/openloopcoronarybc.py b/svzerodsolver/model/openloopcoronarybc.py deleted file mode 100644 index 7f0b76733..000000000 --- a/svzerodsolver/model/openloopcoronarybc.py +++ /dev/null @@ -1,175 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""This module holds the OpenLoopCoronaryBC class.""" -from typing import Sequence - -import numpy as np - -from .block import Block - - -class OpenLoopCoronaryBC(Block): - """Open-loop coronary boundary condition. :cite:`kim_coronary` - - Attributes: - name: Name of the block. - inflow_nodes: Inflow nodes of the element. - outflow_nodes: Outflow nodes of the element. - """ - - _NUM_EQUATIONS = 2 - _NUM_INTERNAL_VARS = 1 - - def __init__(self, params: dict = None, name: str = None): - """Create a new OpenLoopCoronaryBC instance. - - Args: - params: The configuration paramaters of the block. Mostly comprised - of constants for element contribution calculation. - name: Optional name of the block. - """ - super().__init__(params=params, name=name) - - self._pv_func = None - if isinstance(self._params["Pv"], Sequence): - self._pv_func = self._interpolate( - self._params["time"], self._params["Pv"] - ) - - self._pim_func = None - if isinstance(self._params["Pim"], Sequence): - self._pim_func = self._interpolate( - self._params["time"], self._params["Pim"] - ) - - self._vec["C"] = np.zeros(2) - self._mat["F"] = np.zeros((2, 3)) - - self._mat["F"][0, 2] = -1.0 - - Cim_Rv = self._params["Cim"] * self._params["Rv"] - self._mat["E"] = np.zeros((2, 3)) - self._mat["E"][0, 0] = -self._params["Ca"] * Cim_Rv - self._mat["E"][0, 1] = self._params["Ra"] * self._params["Ca"] * Cim_Rv - self._mat["E"][0, 2] = -Cim_Rv - self._mat["E"][1, 2] = -Cim_Rv * self._params["Ram"] - self._mat["F"][0, 1] = Cim_Rv - self._mat["F"][1, 0] = Cim_Rv - self._mat["F"][1, 1] = -Cim_Rv * self._params["Ra"] - self._mat["F"][1, 2] = -(self._params["Rv"] + self._params["Ram"]) - - if self._pv_func is None and self._pim_func is None: - self._vec["C"][0] = ( - -self._params["Cim"] * self._params["Pim"] - + self._params["Cim"] * self._params["Pv"] - ) - # Pa is assumed to be 0.0 - self._vec["C"][1] = ( - -self._params["Cim"] - * (self._params["Rv"] + self._params["Ram"]) - * self._params["Pim"] - + self._params["Ram"] - * self._params["Cim"] - * self._params["Pv"] - ) - self.update_time = super().update_time - elif self._pv_func is None: - self._pv_func = lambda _: self._params["Pv"] - elif self._pim_func is None: - self._pim_func = lambda _: self._params["Pim"] - - if self._params["steady"]: - del self._mat["E"] - self._mat["F"] = np.array( - [ - [ - -self._params["Cim"], - self._params["Cim"] - * (self._params["Ra"] + self._params["Ram"]), - 1.0, - ], - [ - -1.0, - self._params["Ra"] - + self._params["Ram"] - + self._params["Rv"], - 0.0, - ], - ] - ) - self._vec["C"][0] = -self._params["Cim"] * self._params["Pim"] - self._vec["C"][1] = self._params["Pv"] - - @classmethod - def from_config(cls, config): - """Create block from config dictionary. - - Args: - config: The configuration dict for the block. - - Returns: - The block instance. - """ - params = dict( - time=config["bc_values"].get("t", None), - Ra=config["bc_values"]["Ra1"], - Ca=config["bc_values"]["Ca"], - Ram=config["bc_values"]["Ra2"], - Cim=config["bc_values"]["Cc"], - Rv=config["bc_values"]["Rv1"], - Pim=config["bc_values"]["Pim"], - Pv=config["bc_values"]["P_v"], - steady=config["steady"], - ) - - return cls(params=params, name=config["name"]) - - def update_time(self, time): - """Update time dependent element contributions. - - Args: - time: Current time. - """ - # For this open-loop coronary BC, the ordering of solution unknowns is : - # (P_in, Q_in, V_im) where V_im is the volume of the second capacitor, Cim - # Q_in is the flow through the first resistor - # and P_in is the pressure at the inlet of the first resistor - Pim_value = self._pim_func(time) - Pv_value = self._pv_func(time) - self._vec["C"][0] = ( - -self._params["Cim"] * Pim_value + self._params["Cim"] * Pv_value - ) - # Pa is assumed to be 0.0 - self._vec["C"][1] = ( - -self._params["Cim"] - * (self._params["Rv"] + self._params["Ram"]) - * Pim_value - + self._params["Ram"] * self._params["Cim"] * Pv_value - ) diff --git a/svzerodsolver/model/pressurereferencebc.py b/svzerodsolver/model/pressurereferencebc.py deleted file mode 100644 index eb561b32b..000000000 --- a/svzerodsolver/model/pressurereferencebc.py +++ /dev/null @@ -1,94 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""This module holds the PressureReferenceBC class.""" -from typing import Sequence - -import numpy as np - -from .block import Block - - -class PressureReferenceBC(Block): - """Pressure reference boundary condition. - - Applies a predefined pressure at a boundary. - - Attributes: - name: Name of the block. - inflow_nodes: Inflow nodes of the element. - outflow_nodes: Outflow nodes of the element. - """ - - _NUM_EQUATIONS = 1 - - def __init__(self, params: dict = None, name: str = None): - """Create a new PressureReferenceBC instance. - - Args: - params: The configuration paramaters of the block. Mostly comprised - of constants for element contribution calculation. - name: Optional name of the block. - """ - super().__init__(params=params, name=name) - - if isinstance(self._params["P"], Sequence): - self._p_func = self._interpolate( - self._params["time"], self._params["P"] - ) - self._vec["C"] = np.zeros(1, dtype=float) - else: - self._vec["C"] = np.array([-self._params["P"]], dtype=float) - self.update_time = super().update_time - - self._mat["F"] = np.array([[1.0, 0.0]], dtype=float) - - @classmethod - def from_config(cls, config): - """Create block from config dictionary. - - Args: - config: The configuration dict for the block. - - Returns: - The block instance. - """ - params = dict( - time=config["bc_values"].get("t", None), - P=config["bc_values"].get("P"), - ) - return cls(params=params, name=config["name"]) - - def update_time(self, time): - """Update time dependent element contributions. - - Args: - time: Current time. - """ - self._vec["C"][0] = -self._p_func(time) diff --git a/svzerodsolver/model/resistancebc.py b/svzerodsolver/model/resistancebc.py deleted file mode 100644 index 079f48040..000000000 --- a/svzerodsolver/model/resistancebc.py +++ /dev/null @@ -1,104 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""This module holds the ResistanceBC class.""" -from typing import Sequence - -import numpy as np - -from .block import Block - - -class ResistanceBC(Block): - """Resistance boundary condition. - - Attributes: - name: Name of the block. - inflow_nodes: Inflow nodes of the element. - outflow_nodes: Outflow nodes of the element. - """ - - _NUM_EQUATIONS = 1 - - def __init__(self, params: dict = None, name: str = None): - """Create a new ResistanceBC instance. - - Args: - params: The configuration paramaters of the block. Mostly comprised - of constants for element contribution calculation. - name: Optional name of the block. - """ - super().__init__(params=params, name=name) - - self._r_func = None - if isinstance(self._params["R"], Sequence): - self._r_func = self._interpolate( - self._params["time"], self._params["R"] - ) - self._mat["F"] = np.array([[1.0, 0.0]], dtype=float) - else: - self._mat["F"] = np.array([[1.0, -self._params["R"]]], dtype=float) - - self._pd_func = None - if isinstance(self._params["Pd"], Sequence): - self._pd_func = self._interpolate( - self._params["time"], self._params["Pd"] - ) - self._vec["C"] = np.array([0.0], dtype=float) - else: - self._vec["C"] = np.array([-self._params["Pd"]], dtype=float) - - if self._r_func is None and self._pd_func is None: - self.update_time = super().update_time - - @classmethod - def from_config(cls, config): - """Create block from config dictionary. - - Args: - config: The configuration dict for the block. - - Returns: - The block instance. - """ - params = dict( - time=config["bc_values"].get("t", None), - Pd=config["bc_values"].get("Pd"), - R=config["bc_values"].get("R"), - ) - return cls(params=params, name=config["name"]) - - def update_time(self, time): - """Update time dependent element contributions. - - Args: - time: Current time. - """ - self._mat["F"][0, 1] = -self._r_func(time) - self._vec["C"][0] = -self._pd_func(time) diff --git a/svzerodsolver/model/windkesselbc.py b/svzerodsolver/model/windkesselbc.py deleted file mode 100644 index b0c75023c..000000000 --- a/svzerodsolver/model/windkesselbc.py +++ /dev/null @@ -1,128 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""This module holds the WindkesselBC class.""" -from typing import Sequence - -import numpy as np - -from .block import Block - - -class WindkesselBC(Block): - """Windkessel RCR boundary condition. - - Attributes: - name: Name of the block. - inflow_nodes: Inflow nodes of the element. - outflow_nodes: Outflow nodes of the element. - """ - - _NUM_EQUATIONS = 2 - _NUM_INTERNAL_VARS = 1 - - def __init__(self, params: dict = None, name: str = None): - """Create a new WindkesselBC instance. - - Args: - params: The configuration paramaters of the block. Mostly comprised - of constants for element contribution calculation. - name: Optional name of the block. - """ - super().__init__(params=params, name=name) - - self._mat["E"] = np.zeros((2, 3), dtype=float) - self._mat["F"] = np.array( - [[1.0, 0.0, -1.0], [0.0, 0.0, -1.0]], dtype=float - ) - self._vec["C"] = np.array([0.0, 0.0], dtype=float) - - self._rp_func = None - if isinstance(self._params["Rp"], Sequence): - self._rp_func = self._interpolate( - self._params["time"], self._params["Rp"] - ) - - self._c_func = None - if isinstance(self._params["C"], Sequence): - self._c_func = self._interpolate( - self._params["time"], self._params["C"] - ) - - self._rd_func = None - if isinstance(self._params["Rd"], Sequence): - self._rd_func = self._interpolate( - self._params["time"], self._params["Rd"] - ) - - self._pd_func = None - if isinstance(self._params["Pd"], Sequence): - self._pd_func = self._interpolate( - self._params["time"], self._params["Pd"] - ) - - if [self._rp_func, self._c_func, self._rd_func, self._pd_func] == [ - None - ] * 4: - self._mat["E"][1, 2] = -self._params["Rd"] * self._params["C"] - self._mat["F"][0, 1] = -self._params["Rp"] - self._mat["F"][1, 1] = self._params["Rd"] - self._vec["C"][1] = self._params["Pd"] - self.update_time = super().update_time - - @classmethod - def from_config(cls, config): - """Create block from config dictionary. - - Args: - config: The configuration dict for the block. - - Returns: - The block instance. - """ - params = dict( - time=config["bc_values"].get("t", None), - Rp=config["bc_values"].get("Rp"), - C=config["bc_values"].get("C"), - Rd=config["bc_values"].get("Rd"), - Pd=config["bc_values"].get("Pd"), - ) - return cls(params=params, name=config["name"]) - - def update_time(self, time): - """Update time dependent element contributions. - - Args: - time: Current time. - """ - rd_t = self._rd_func(time) - self._mat["E"][1, 2] = -rd_t * self._c_func(time) - self._mat["F"][0, 1] = -self._rp_func(time) - self._mat["F"][1, 1] = rd_t - self._vec["C"][1] = self._pd_func(time) diff --git a/svzerodsolver/runner.py b/svzerodsolver/runner.py deleted file mode 100644 index 7e7287f94..000000000 --- a/svzerodsolver/runner.py +++ /dev/null @@ -1,158 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""This module holds the main execution routines of svZeroDSolver.""" -import json -import pandas as pd -import numpy as np - -import click - -from .algebra import run_integrator -from .utils import ( - convert_unsteady_bcs_to_steady, - create_blocks, - format_results_to_dict, - get_solver_params, -) - - -def run_from_config(parameters): - """Run the svZeroDSolver. - - Args: - config: Python dict of the configuration. - - Returns: - Python dict with results. - """ - - y_initial = None - ydot_initial = None - sim_params = parameters["simulation_parameters"] - if sim_params.get("steady_initial", True): - steady_parameters = convert_unsteady_bcs_to_steady(parameters) - - # to run the 0d model with steady BCs to steady-state, simulate this model - # with large time step size for an arbitrarily small number of cardiac cycles - block_list, dofhandler = create_blocks(steady_parameters, steady=True) - time_step_size, num_time_steps = get_solver_params(steady_parameters) - (_, y_list, ydot_list) = run_integrator( - block_list, - dofhandler, - num_time_steps, - time_step_size, - atol=sim_params.get("absolute_tolerance", 1e-8), - max_iter=sim_params.get("maximum_nonlinear_iterations", 30), - ) - y_initial = y_list[-1] - ydot_initial = ydot_list[-1] - - block_list, dofhandler = create_blocks(parameters) - time_step_size, num_time_steps = get_solver_params(parameters) - (time_steps, y_list, ydot_list) = run_integrator( - block_list, - dofhandler, - num_time_steps, - time_step_size, - y_initial, - ydot_initial, - atol=sim_params.get("absolute_tolerance", 1e-8), - max_iter=sim_params.get("maximum_nonlinear_iterations", 30), - ) - - branch_result = format_results_to_dict(time_steps, y_list, block_list) - - result = pd.DataFrame( - columns=[ - "name", - "time", - "flow_in", - "flow_out", - "pressure_in", - "pressure_out", - ] - ) - for branch_id, name in enumerate(branch_result["names"]): - result = pd.concat( - [ - result, - pd.DataFrame.from_dict( - { - "name": [name] * len(branch_result["time"]), - "time": np.array(branch_result["time"]), - "flow_in": np.array( - branch_result["flow_in"][branch_id] - ), - "flow_out": np.array( - branch_result["flow_out"][branch_id] - ), - "pressure_in": np.array( - branch_result["pressure_in"][branch_id] - ), - "pressure_out": np.array( - branch_result["pressure_out"][branch_id] - ), - } - ), - ] - ) - return result - - -def run_from_file(input_file, output_file): - """Run the svZeroDSolver from file. - - Args: - input_file: Input file with configuration. - output_file: Output file with configuration. - """ - with open(input_file) as ff: - config = json.load(ff) - result = run_from_config(config) - with open(output_file, "w") as ff: - json.dump(result, ff) - - -@click.command() -@click.option( - "--input_file", - help="Path to the svZeroDSolver input file.", - required=True, - type=str, -) -@click.option( - "--output_file", - help="Path to the svZeroDSolver output file.", - required=True, - type=str, -) -def _run_from_from_sys_args(input_file, output_file): - """Run the svZeroDSolver.""" - run_from_file(input_file, output_file) diff --git a/svzerodsolver/runnercpp.py b/svzerodsolver/runnercpp.py deleted file mode 100644 index f7d081cdd..000000000 --- a/svzerodsolver/runnercpp.py +++ /dev/null @@ -1,48 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""Module to call the C++ version of svZeroDSolver. - -Expects to find a build in the Release folder. -""" -from io import StringIO -import svzerodplus - -from pandas import read_csv - - -def run_from_config(config): - """Run the C++ svZeroDSolver. - - Args: - config: Python dict of the configuration. - - Returns: - Pandas dataframe with the results.""" - return read_csv(StringIO(svzerodplus.simulate(config))) diff --git a/svzerodsolver/utils.py b/svzerodsolver/utils.py deleted file mode 100644 index 4271c561e..000000000 --- a/svzerodsolver/utils.py +++ /dev/null @@ -1,270 +0,0 @@ -# Copyright (c) Stanford University, The Regents of the University of -# California, and others. -# -# All Rights Reserved. -# -# See Copyright-SimVascular.txt for additional details. -# -# Permission is hereby granted, free of charge, to any person obtaining -# a copy of this software and associated documentation files (the -# "Software"), to deal in the Software without restriction, including -# without limitation the rights to use, copy, modify, merge, publish, -# distribute, sublicense, and/or sell copies of the Software, and to -# permit persons to whom the Software is furnished to do so, subject -# to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""This module holds helper functions for svZeroDSolver.""" -from copy import deepcopy - -import numpy as np - -from .model.bloodvessel import BloodVessel -from .model.dofhandler import DOFHandler -from .model.flowreferencebc import FlowReferenceBC -from .model.internaljunction import InternalJunction -from .model.node import Node -from .model.openloopcoronarybc import OpenLoopCoronaryBC -from .model.pressurereferencebc import PressureReferenceBC -from .model.resistancebc import ResistanceBC -from .model.windkesselbc import WindkesselBC - - -def get_solver_params(config): - """Determine time step size and number of timesteps. - - Args: - config: svZeroDSolver configuration. - - Returns: - time_step_size: Time step size: - num_time_steps: Number of time steps. - """ - sim_params = config["simulation_parameters"] - cardiac_cycle_period = sim_params.get("cardiac_cycle_period", 1.0) - num_cycles = sim_params.get("number_of_cardiac_cycles") - num_pts_per_cycle = sim_params.get("number_of_time_pts_per_cardiac_cycle") - time_step_size = cardiac_cycle_period / (num_pts_per_cycle - 1) - num_time_steps = int((num_pts_per_cycle - 1) * num_cycles + 1) - return time_step_size, num_time_steps - - -def convert_unsteady_bcs_to_steady(config): - """Convert unsteady boundary conditions to their steady equivalent. - - Args: - config: svZeroDSolver configuration. - - Returns: - steady_config: Configuration of the steady equivalent. - """ - steady_config = deepcopy(config) - steady_config["simulation_parameters"][ - "number_of_time_pts_per_cardiac_cycle" - ] = 11 - steady_config["simulation_parameters"]["number_of_cardiac_cycles"] = 3 - bc_identifiers = {"FLOW": "Q", "PRESSURE": "P", "CORONARY": "Pim"} - for i, bc in enumerate(config["boundary_conditions"]): - if bc["bc_type"] in bc_identifiers: - bc_values = bc["bc_values"][bc_identifiers[bc["bc_type"]]] - # Time averaged value for a single cariadic_cycle - del steady_config["boundary_conditions"][i]["bc_values"]["t"] - steady_config["boundary_conditions"][i]["bc_values"][ - bc_identifiers[bc["bc_type"]] - ] = np.mean(bc_values) - if bc["bc_type"] == "RCR": - steady_config["boundary_conditions"][i]["bc_values"]["C"] = 0.0 - - return steady_config - - -def create_blocks(config, steady=False): - """Create blocks. - - Args: - config: svZeroDSolver configuration. - steady: Toggle if blocks should be created with steady behavior. - - Returns: - blocks: List of created blocks. - dofhandler: Degree-of-freedom handler. - """ - block_dict = {} - connections = [] - - # Create junctions - for junction_config in config["junctions"]: - if junction_config["junction_type"] in [ - "NORMAL_JUNCTION", - "internal_junction", - ]: - junction = InternalJunction.from_config(junction_config) - else: - raise ValueError( - f"Unknown junction type: {junction_config['junction_type']}" - ) - connections += [ - (f"V{vid}", junction.name) - for vid in junction_config["inlet_vessels"] - ] - connections += [ - (junction.name, f"V{vid}") - for vid in junction_config["outlet_vessels"] - ] - if junction.name in block_dict: - raise RuntimeError(f"Junction {junction.name} already exists.") - block_dict[junction.name] = junction - - # Create vessels and boundary conditions - for vessel_config in config["vessels"]: - if vessel_config["zero_d_element_type"] == "BloodVessel": - vessel = BloodVessel.from_config(vessel_config) - else: - raise NotImplementedError - if "boundary_conditions" in vessel_config: - if "inlet" in vessel_config["boundary_conditions"]: - connections.append( - ( - "BC" + str(vessel_config["vessel_id"]) + "_inlet", - vessel.name, - ) - ) - if "outlet" in vessel_config["boundary_conditions"]: - connections.append( - ( - vessel.name, - "BC" + str(vessel_config["vessel_id"]) + "_outlet", - ) - ) - if vessel.name in block_dict: - raise RuntimeError(f"Vessel {vessel.name} already exists.") - block_dict[vessel.name] = vessel - if "boundary_conditions" in vessel_config: - locations = [ - loc - for loc in ("inlet", "outlet") - if loc in vessel_config["boundary_conditions"] - ] - for location in locations: - vessel_id = vessel_config["vessel_id"] - for bc_config in config["boundary_conditions"]: - if ( - bc_config["bc_name"] - == vessel_config["boundary_conditions"][location] - ): - bc_config = dict( - name="BC" + str(vessel_id) + "_" + location, - steady=steady, - **bc_config, - ) - break - - if ( - "t" in bc_config["bc_values"] - and len(bc_config["bc_values"]["t"]) >= 2 - ): - cardiac_cycle_period = ( - bc_config["bc_values"]["t"][-1] - - bc_config["bc_values"]["t"][0] - ) - if ( - "cardiac_cycle_period" - in config["simulation_parameters"] - and cardiac_cycle_period - != config["simulation_parameters"][ - "cardiac_cycle_period" - ] - ): - raise RuntimeError( - f"The time series of the boundary condition for segment {vessel_id} does " - "not have the same cardiac cycle period as the other boundary conditions." - ) - elif ( - "cardiac_cycle_period" - not in config["simulation_parameters"] - ): - config["simulation_parameters"][ - "cardiac_cycle_period" - ] = cardiac_cycle_period - - if bc_config["bc_type"] == "RESISTANCE": - bc = ResistanceBC.from_config(bc_config) - elif bc_config["bc_type"] == "RCR": - bc = WindkesselBC.from_config(bc_config) - elif bc_config["bc_type"] == "FLOW": - bc = FlowReferenceBC.from_config(bc_config) - elif bc_config["bc_type"] == "PRESSURE": - bc = PressureReferenceBC.from_config(bc_config) - elif bc_config["bc_type"] == "CORONARY": - bc = OpenLoopCoronaryBC.from_config(bc_config) - else: - raise NotImplementedError - block_dict[bc.name] = bc - - # Create nodes - dofhandler = DOFHandler() - for ele1_name, ele2_name in connections: - node = Node( - block_dict[ele1_name], - block_dict[ele2_name], - name=ele1_name + "_" + ele2_name, - ) - node.setup_dofs(dofhandler) - - # Setup degrees of freedom of the model - for key in sorted(block_dict): - block_dict[key].setup_dofs(dofhandler) - return list(block_dict.values()), dofhandler - - -def format_results_to_dict(time_steps, result_array, block_list): - """Format result array to dict format. - - Args: - time_steps: List of time steps. - - Returns: - results: Resuluts in dict format. - """ - - result_array = np.array(result_array) - vessels = [block for block in block_list if isinstance(block, BloodVessel)] - results = { - "flow_in": [], - "flow_out": [], - "names": [], - "pressure_in": [], - "pressure_out": [], - "time": list(time_steps), - } - - for vessel in vessels: - - results["names"].append(vessel.name) - results["flow_in"].append( - list(result_array[:, vessel.inflow_nodes[0].flow_dof]) - ) - results["flow_out"].append( - list(result_array[:, vessel.outflow_nodes[0].flow_dof]) - ) - results["pressure_in"].append( - list(result_array[:, vessel.inflow_nodes[0].pres_dof]) - ) - results["pressure_out"].append( - list(result_array[:, vessel.outflow_nodes[0].pres_dof]) - ) - - return results diff --git a/tests/conftest.py b/tests/conftest.py deleted file mode 100644 index e41d4114a..000000000 --- a/tests/conftest.py +++ /dev/null @@ -1,11 +0,0 @@ -from tempfile import TemporaryDirectory - -import numpy as np -import pytest - - -@pytest.fixture -def tempdir(): - """Temporary directory for test purposes.""" - with TemporaryDirectory() as tempdir: - yield tempdir diff --git a/tests/test_integration.py b/tests/test_integration.py index 645c4381b..34a2d28f5 100644 --- a/tests/test_integration.py +++ b/tests/test_integration.py @@ -2,371 +2,545 @@ import os import numpy as np +import pytest -from svzerodsolver.runner import run_from_config +import svzerodplus + +this_file_dir = os.path.abspath(os.path.dirname(__file__)) RTOL_PRES = 1.0e-7 RTOL_FLOW = 1.0e-8 -def run_test_case_by_name(name, testdir): +def run_test_case_by_name(name, output_variable_based=False, folder="."): """Run a test case by its case name. Args: name: Name of the test case. testdir: Directory for performing the simulation. """ - testfile = os.path.join(os.path.dirname(__file__), "cases", name + ".json") + testfile = os.path.join(this_file_dir, "cases", folder, name + ".json") with open(testfile) as ff: config = json.load(ff) - result = run_from_config(config) - - output = { - "pressure_in": {}, - "pressure_out": {}, - "flow_in": {}, - "flow_out": {}, - } - - last_seg_id = 0 - - for vessel in config["vessels"]: - name = vessel["vessel_name"] - branch_id, seg_id = name.split("_") - branch_id, seg_id = int(branch_id[6:]), int(seg_id[3:]) - vessel_id = vessel["vessel_id"] - - if seg_id == 0: - output["pressure_in"][branch_id] = np.array( - result[result.name == f"V{vessel_id}"]["pressure_in"] - ) - output["flow_in"][branch_id] = np.array( - result[result.name == f"V{vessel_id}"]["flow_in"] - ) - output["pressure_out"][branch_id] = np.array( - result[result.name == f"V{vessel_id}"]["pressure_out"] - ) - output["flow_out"][branch_id] = np.array( - result[result.name == f"V{vessel_id}"]["flow_out"] - ) - elif seg_id > last_seg_id: - output["pressure_out"][branch_id] = np.array( - result[result.name == f"V{vessel_id}"]["pressure_out"] - ) - output["flow_out"][branch_id] = np.array( - result[result.name == f"V{vessel_id}"]["flow_out"] - ) - - last_seg_id = seg_id + result = svzerodplus.simulate(config) + + if output_variable_based == False: + output = { + "pressure_in": {}, + "pressure_out": {}, + "flow_in": {}, + "flow_out": {}, + } + + last_seg_id = 0 + + for vessel in config["vessels"]: + name = vessel["vessel_name"] + branch_id, seg_id = name.split("_") + branch_id, seg_id = int(branch_id[6:]), int(seg_id[3:]) + vessel_id = vessel["vessel_id"] + + if seg_id == 0: + output["pressure_in"][branch_id] = np.array( + result[result.name == name]["pressure_in"] + ) + output["flow_in"][branch_id] = np.array( + result[result.name == name]["flow_in"] + ) + output["pressure_out"][branch_id] = np.array( + result[result.name == name]["pressure_out"] + ) + output["flow_out"][branch_id] = np.array( + result[result.name == name]["flow_out"] + ) + elif seg_id > last_seg_id: + output["pressure_out"][branch_id] = np.array( + result[result.name == name]["pressure_out"] + ) + output["flow_out"][branch_id] = np.array( + result[result.name == name]["flow_out"] + ) + + last_seg_id = seg_id + + elif output_variable_based == True: + output = result return output -def get_result(result_array, field, branch, branch_node, time_step): +def get_result(result_array, field, branch, time_step): """ "Get results at specific field, branch, branch_node and time step.""" # extract result - node_conversion = ["in", "out"] - field_new = f"{field}_" + node_conversion[branch_node] - return result_array[field_new][branch][time_step] + return result_array[field][branch][time_step] -def test_steady_flow_R_R(tmpdir): - results = run_test_case_by_name("steadyFlow_R_R", tmpdir) +def test_steady_flow_R_R(): + results = run_test_case_by_name("steadyFlow_R_R") assert np.isclose( - get_result(results, "pressure", 0, 0, -1), 1100.0, rtol=RTOL_PRES + get_result(results, "pressure_in", 0, -1), 1100.0, rtol=RTOL_PRES ) # inlet pressure assert np.isclose( - get_result(results, "pressure", 0, -1, -1), 600.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 0, -1), 600.0, rtol=RTOL_PRES ) # outlet pressure assert np.isclose( - get_result(results, "flow", 0, 0, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW ) # inlet flow assert np.isclose( - get_result(results, "flow", 0, -1, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW ) # outlet flow -def test_steady_flow_r_coronary(tmpdir): - results = run_test_case_by_name("steadyFlow_R_coronary", tmpdir) +def test_steady_flow_r_coronary(): + results = run_test_case_by_name("steadyFlow_R_coronary") assert np.isclose( - get_result(results, "pressure", 0, 0, -1), 2000.0, rtol=RTOL_PRES + get_result(results, "pressure_in", 0, -1), 2000.0, rtol=RTOL_PRES ) # inlet pressure assert np.isclose( - get_result(results, "pressure", 0, -1, -1), 1500.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 0, -1), 1500.0, rtol=RTOL_PRES ) # outlet pressure assert np.isclose( - get_result(results, "flow", 0, 0, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW ) # inlet flow assert np.isclose( - get_result(results, "flow", 0, -1, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW ) # outlet flow -def test_steady_flow_rlc_r(tmpdir): - results = run_test_case_by_name("steadyFlow_RLC_R", tmpdir) +def test_steady_flow_rlc_r(): + results = run_test_case_by_name("steadyFlow_RLC_R") assert np.isclose( - get_result(results, "pressure", 0, 0, -1), 1100.0, rtol=RTOL_PRES + get_result(results, "pressure_in", 0, -1), 1100.0, rtol=RTOL_PRES ) # inlet pressure assert np.isclose( - get_result(results, "pressure", 0, -1, -1), 600.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 0, -1), 600.0, rtol=RTOL_PRES ) # outlet pressure assert np.isclose( - get_result(results, "flow", 0, 0, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW ) # inlet flow assert np.isclose( - get_result(results, "flow", 0, -1, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW ) # outlet flow -def test_steady_flow_rc_r(tmpdir): - results = run_test_case_by_name("steadyFlow_RC_R", tmpdir) +def test_steady_flow_rc_r(): + results = run_test_case_by_name("steadyFlow_RC_R") assert np.isclose( - get_result(results, "pressure", 0, 0, -1), 1100.0, rtol=RTOL_PRES + get_result(results, "pressure_in", 0, -1), 1100.0, rtol=RTOL_PRES ) # inlet pressure assert np.isclose( - get_result(results, "pressure", 0, -1, -1), 600.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 0, -1), 600.0, rtol=RTOL_PRES ) # outlet pressure assert np.isclose( - get_result(results, "flow", 0, 0, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW ) # inlet flow assert np.isclose( - get_result(results, "flow", 0, -1, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW ) # outlet flow -def test_steady_flow_rl_r(tmpdir): - results = run_test_case_by_name("steadyFlow_RL_R", tmpdir) +def test_steady_flow_rl_r(): + results = run_test_case_by_name("steadyFlow_RL_R") assert np.isclose( - get_result(results, "pressure", 0, 0, -1), 1100.0, rtol=RTOL_PRES + get_result(results, "pressure_in", 0, -1), 1100.0, rtol=RTOL_PRES ) # inlet pressure assert np.isclose( - get_result(results, "pressure", 0, -1, -1), 600.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 0, -1), 600.0, rtol=RTOL_PRES ) # outlet pressure assert np.isclose( - get_result(results, "flow", 0, 0, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW ) # inlet flow assert np.isclose( - get_result(results, "flow", 0, -1, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW ) # outlet flow -def test_steady_flow_r_rcr(tmpdir): - results = run_test_case_by_name("steadyFlow_R_RCR", tmpdir) +def test_steady_flow_r_rcr(): + results = run_test_case_by_name("steadyFlow_R_RCR") assert np.isclose( - get_result(results, "pressure", 0, 0, -1), 10500.0, rtol=RTOL_PRES + get_result(results, "pressure_in", 0, -1), 10500.0, rtol=RTOL_PRES ) # inlet pressure assert np.isclose( - get_result(results, "pressure", 0, -1, -1), 10000.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 0, -1), 10000.0, rtol=RTOL_PRES ) # outlet pressure assert np.isclose( - get_result(results, "flow", 0, 0, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW ) # inlet flow assert np.isclose( - get_result(results, "flow", 0, -1, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW ) # outlet flow -def test_steady_flow_r_steady_pressure(tmpdir): - results = run_test_case_by_name("steadyFlow_R_steadyPressure", tmpdir) +def test_steady_flow_r_steady_pressure(): + results = run_test_case_by_name("steadyFlow_R_steadyPressure") assert np.isclose( - get_result(results, "pressure", 0, 0, -1), 1500.0, rtol=RTOL_PRES + get_result(results, "pressure_in", 0, -1), 1500.0, rtol=RTOL_PRES ) # inlet pressure assert np.isclose( - get_result(results, "pressure", 0, -1, -1), 1000.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 0, -1), 1000.0, rtol=RTOL_PRES ) # outlet pressure assert np.isclose( - get_result(results, "flow", 0, 0, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW ) # inlet flow assert np.isclose( - get_result(results, "flow", 0, -1, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW ) # outlet flow -def test_steady_flow_stenosis_r(tmpdir): - results = run_test_case_by_name("steadyFlow_stenosis_R", tmpdir) +def test_steady_flow_stenosis_r(): + results = run_test_case_by_name("steadyFlow_stenosis_R") assert np.isclose( - get_result(results, "pressure", 0, 0, -1), 3600.0, rtol=RTOL_PRES + get_result(results, "pressure_in", 0, -1), 3600.0, rtol=RTOL_PRES ) # inlet pressure assert np.isclose( - get_result(results, "pressure", 0, -1, -1), 600.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 0, -1), 600.0, rtol=RTOL_PRES ) # outlet pressure assert np.isclose( - get_result(results, "flow", 0, 0, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW ) # inlet flow assert np.isclose( - get_result(results, "flow", 0, -1, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW ) # outlet flow -def test_steady_flow_bifurcationr_r1(tmpdir): - results = run_test_case_by_name("steadyFlow_bifurcationR_R1", tmpdir) +def test_steady_flow_bifurcationr_r1(): + results = run_test_case_by_name("steadyFlow_bifurcationR_R1") assert np.isclose( - get_result(results, "pressure", 0, 0, -1), 1100.0, rtol=RTOL_PRES + get_result(results, "pressure_in", 0, -1), 1100.0, rtol=RTOL_PRES ) # parent inlet pressure assert np.isclose( - get_result(results, "pressure", 0, -1, -1), 600.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 0, -1), 600.0, rtol=RTOL_PRES ) # parent outlet pressure assert np.isclose( - get_result(results, "pressure", 1, 0, -1), 600.0, rtol=RTOL_PRES + get_result(results, "pressure_in", 1, -1), 600.0, rtol=RTOL_PRES ) # daughter1 inlet pressure assert np.isclose( - get_result(results, "pressure", 1, -1, -1), 350.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 1, -1), 350.0, rtol=RTOL_PRES ) # daughter1 outlet pressure assert np.isclose( - get_result(results, "pressure", 2, 0, -1), 600.0, rtol=RTOL_PRES + get_result(results, "pressure_in", 2, -1), 600.0, rtol=RTOL_PRES ) # daughter2 inlet pressure assert np.isclose( - get_result(results, "pressure", 2, -1, -1), 350.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 2, -1), 350.0, rtol=RTOL_PRES ) # daughter2 outlet pressure assert np.isclose( - get_result(results, "flow", 0, 0, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW ) # parent inlet flow assert np.isclose( - get_result(results, "flow", 0, -1, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW ) # parent outlet flow assert np.isclose( - get_result(results, "flow", 1, 0, -1), 2.5, rtol=RTOL_FLOW + get_result(results, "flow_in", 1, -1), 2.5, rtol=RTOL_FLOW ) # daughter1 inlet flow assert np.isclose( - get_result(results, "flow", 1, -1, -1), 2.5, rtol=RTOL_FLOW + get_result(results, "flow_out", 1, -1), 2.5, rtol=RTOL_FLOW ) # daughter1 outlet flow assert np.isclose( - get_result(results, "flow", 2, 0, -1), 2.5, rtol=RTOL_FLOW + get_result(results, "flow_in", 2, -1), 2.5, rtol=RTOL_FLOW ) # daughter2 inlet flow assert np.isclose( - get_result(results, "flow", 2, -1, -1), 2.5, rtol=RTOL_FLOW + get_result(results, "flow_out", 2, -1), 2.5, rtol=RTOL_FLOW ) # daughter2 outlet flow -def test_steady_flow_bifurcationr_r2(tmpdir): - results = run_test_case_by_name("steadyFlow_bifurcationR_R2", tmpdir) +def test_steady_flow_bifurcationr_r2(): + results = run_test_case_by_name("steadyFlow_bifurcationR_R2") assert np.isclose( - get_result(results, "pressure", 0, 0, -1), 3462.5, rtol=RTOL_PRES + get_result(results, "pressure_in", 0, -1), 3462.5, rtol=RTOL_PRES ) # parent inlet pressure assert np.isclose( - get_result(results, "pressure", 0, -1, -1), 1962.5, rtol=RTOL_PRES + get_result(results, "pressure_out", 0, -1), 1962.5, rtol=RTOL_PRES ) # parent outlet pressure assert np.isclose( - get_result(results, "pressure", 1, 0, -1), 1962.5, rtol=RTOL_PRES + get_result(results, "pressure_in", 1, -1), 1962.5, rtol=RTOL_PRES ) # daughter1 inlet pressure assert np.isclose( - get_result(results, "pressure", 1, -1, -1), 432.5, rtol=RTOL_PRES + get_result(results, "pressure_out", 1, -1), 432.5, rtol=RTOL_PRES ) # daughter1 outlet pressure assert np.isclose( - get_result(results, "pressure", 2, 0, -1), 1962.5, rtol=RTOL_PRES + get_result(results, "pressure_in", 2, -1), 1962.5, rtol=RTOL_PRES ) # daughter2 inlet pressure assert np.isclose( - get_result(results, "pressure", 2, -1, -1), 1375.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 2, -1), 1375.0, rtol=RTOL_PRES ) # daughter2 outlet pressure assert np.isclose( - get_result(results, "flow", 0, 0, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW ) # parent inlet flow assert np.isclose( - get_result(results, "flow", 0, -1, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW ) # parent outlet flow assert np.isclose( - get_result(results, "flow", 1, 0, -1), 3.825, rtol=RTOL_FLOW + get_result(results, "flow_in", 1, -1), 3.825, rtol=RTOL_FLOW ) # daughter1 inlet flow assert np.isclose( - get_result(results, "flow", 1, -1, -1), 3.825, rtol=RTOL_FLOW + get_result(results, "flow_out", 1, -1), 3.825, rtol=RTOL_FLOW ) # daughter1 outlet flow assert np.isclose( - get_result(results, "flow", 2, 0, -1), 1.175, rtol=RTOL_FLOW + get_result(results, "flow_in", 2, -1), 1.175, rtol=RTOL_FLOW ) # daughter2 inlet flow assert np.isclose( - get_result(results, "flow", 2, -1, -1), 1.175, rtol=RTOL_FLOW + get_result(results, "flow_out", 2, -1), 1.175, rtol=RTOL_FLOW ) # daughter2 outlet flow -def test_pulsatile_flow_r_rcr(tmpdir): - results = run_test_case_by_name("pulsatileFlow_R_RCR", tmpdir) +def test_steadyFlow_blood_vessel_junction(): + results = run_test_case_by_name("steadyFlow_blood_vessel_junction") + assert np.isclose( + get_result(results, "pressure_in", 0, -1), 3462.5, rtol=RTOL_PRES + ) # parent inlet pressure + assert np.isclose( + get_result(results, "pressure_out", 0, -1), 1962.5, rtol=RTOL_PRES + ) # parent outlet pressure + assert np.isclose( + get_result(results, "pressure_in", 1, -1), 1580.0, rtol=RTOL_PRES + ) # daughter1 inlet pressure + assert np.isclose( + get_result(results, "pressure_out", 1, -1), 432.5, rtol=RTOL_PRES + ) # daughter1 outlet pressure + assert np.isclose( + get_result(results, "pressure_in", 2, -1), 1727.5, rtol=RTOL_PRES + ) # daughter2 inlet pressure assert np.isclose( - get_result(results, "pressure", 0, 0, 0), 4620.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 2, -1), 1375.0, rtol=RTOL_PRES + ) # daughter2 outlet pressure + assert np.isclose( + get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW + ) # parent inlet flow + assert np.isclose( + get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW + ) # parent outlet flow + assert np.isclose( + get_result(results, "flow_in", 1, -1), 3.825, rtol=RTOL_FLOW + ) # daughter1 inlet flow + assert np.isclose( + get_result(results, "flow_out", 1, -1), 3.825, rtol=RTOL_FLOW + ) # daughter1 outlet flow + assert np.isclose( + get_result(results, "flow_in", 2, -1), 1.175, rtol=RTOL_FLOW + ) # daughter2 inlet flow + assert np.isclose( + get_result(results, "flow_out", 2, -1), 1.175, rtol=RTOL_FLOW + ) # daughter2 outlet flow + + +def test_pulsatile_flow_r_rcr(): + results = run_test_case_by_name("pulsatileFlow_R_RCR") + assert np.isclose( + get_result(results, "pressure_in", 0, 0), 4620.0, rtol=RTOL_PRES ) # inlet pressure assert np.isclose( - get_result(results, "pressure", 0, -1, 0), 4400.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 0, 0), 4400.0, rtol=RTOL_PRES ) # outlet pressure assert np.isclose( - get_result(results, "flow", 0, 0, 0), 2.2, rtol=RTOL_FLOW + get_result(results, "flow_in", 0, 0), 2.2, rtol=RTOL_FLOW ) # inlet flow assert np.isclose( - get_result(results, "flow", 0, -1, 0), 2.2, rtol=RTOL_FLOW + get_result(results, "flow_out", 0, 0), 2.2, rtol=RTOL_FLOW ) # outlet flow -def test_pulsatile_flow_r_coronary(tmpdir): - results = run_test_case_by_name("pulsatileFlow_R_coronary", tmpdir) +def test_pulsatile_flow_r_coronary(): + results = run_test_case_by_name("pulsatileFlow_R_coronary") assert np.isclose( - get_result(results, "pressure", 0, 0, 0), 880.0, rtol=RTOL_PRES + get_result(results, "pressure_in", 0, 0), 880.0, rtol=RTOL_PRES ) # inlet pressure assert np.isclose( - get_result(results, "pressure", 0, -1, 0), 660.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 0, 0), 660.0, rtol=RTOL_PRES ) # outlet pressure assert np.isclose( - get_result(results, "flow", 0, 0, 0), 2.2, rtol=RTOL_FLOW + get_result(results, "flow_in", 0, 0), 2.2, rtol=RTOL_FLOW ) # inlet flow assert np.isclose( - get_result(results, "flow", 0, -1, 0), 2.2, rtol=RTOL_FLOW + get_result(results, "flow_out", 0, 0), 2.2, rtol=RTOL_FLOW ) # outlet flow -def test_pulsatile_flow_cstenosis_steady_pressure(tmpdir): - results = run_test_case_by_name("pulsatileFlow_CStenosis_steadyPressure", tmpdir) +def test_pulsatile_flow_cstenosis_steady_pressure(): + results = run_test_case_by_name("pulsatileFlow_CStenosis_steadyPressure") assert np.isclose( - get_result(results, "pressure", 0, 0, -439), - 0.5937166987044968, + get_result(results, "pressure_in", 0, -439), + 0.5933049197138334, rtol=1.0e-5, ) # inlet pressure assert np.isclose( - get_result(results, "pressure", 0, -1, -439), 0.1, rtol=1.0e-5 + get_result(results, "pressure_out", 0, -439), 0.1, rtol=1.0e-5 ) # outlet pressure assert np.isclose( - get_result(results, "flow", 0, 0, -439), - 0.7026499697830042, + get_result(results, "flow_in", 0, -439), + 0.7023611813029965, rtol=1.0e-5, ) # inlet flow assert np.isclose( - get_result(results, "flow", 0, -1, -439), - 0.7026499697830042, + get_result(results, "flow_out", 0, -439), + 0.7018707542098627, rtol=1.0e-5, ) # outlet flow -def test_steady_flow_confluencer_r(tmpdir): - results = run_test_case_by_name("steadyFlow_confluenceR_R", tmpdir) +def test_steady_flow_confluencer_r(): + results = run_test_case_by_name("steadyFlow_confluenceR_R") assert np.isclose( - get_result(results, "pressure", 0, 0, -1), 6600.0, rtol=RTOL_PRES + get_result(results, "pressure_in", 0, -1), 6600.0, rtol=RTOL_PRES ) # parent inlet pressure assert np.isclose( - get_result(results, "pressure", 0, -1, -1), 6100.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 0, -1), 6100.0, rtol=RTOL_PRES ) # parent outlet pressure assert np.isclose( - get_result(results, "pressure", 1, 0, -1), 8100.0, rtol=RTOL_PRES + get_result(results, "pressure_in", 1, -1), 8100.0, rtol=RTOL_PRES ) # daughter1 inlet pressure assert np.isclose( - get_result(results, "pressure", 1, -1, -1), 6100.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 1, -1), 6100.0, rtol=RTOL_PRES ) # daughter1 outlet pressure assert np.isclose( - get_result(results, "pressure", 2, 0, -1), 6100.0, rtol=RTOL_PRES + get_result(results, "pressure_in", 2, -1), 6100.0, rtol=RTOL_PRES ) # daughter2 inlet pressure assert np.isclose( - get_result(results, "pressure", 2, -1, -1), 1600.0, rtol=RTOL_PRES + get_result(results, "pressure_out", 2, -1), 1600.0, rtol=RTOL_PRES ) # daughter2 outlet pressure assert np.isclose( - get_result(results, "flow", 0, 0, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW ) # parent inlet flow assert np.isclose( - get_result(results, "flow", 0, -1, -1), 5.0, rtol=RTOL_FLOW + get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW ) # parent outlet flow assert np.isclose( - get_result(results, "flow", 1, 0, -1), 10.0, rtol=RTOL_FLOW + get_result(results, "flow_in", 1, -1), 10.0, rtol=RTOL_FLOW ) # daughter1 inlet flow assert np.isclose( - get_result(results, "flow", 1, -1, -1), 10.0, rtol=RTOL_FLOW + get_result(results, "flow_out", 1, -1), 10.0, rtol=RTOL_FLOW ) # daughter1 outlet flow assert np.isclose( - get_result(results, "flow", 2, 0, -1), 15.0, rtol=RTOL_FLOW + get_result(results, "flow_in", 2, -1), 15.0, rtol=RTOL_FLOW ) # daughter2 inlet flow assert np.isclose( - get_result(results, "flow", 2, -1, -1), 15.0, rtol=RTOL_FLOW + get_result(results, "flow_out", 2, -1), 15.0, rtol=RTOL_FLOW ) # daughter2 outlet flow + + +def test_closed_loop_heart_single_vessel(): + results = run_test_case_by_name("closedLoopHeart_singleVessel") + assert np.isclose( + np.mean(np.array(results["pressure_in"][0])), 55.703345704742844, rtol=RTOL_PRES + ) # mean aortic pressure + assert np.isclose( + np.amax(np.array(results["pressure_in"][0])), 73.97450170686889, rtol=RTOL_PRES + ) # max aortic pressure + assert np.isclose( + np.amin(np.array(results["pressure_in"][0])), 0.0, rtol=RTOL_PRES + ) # min aortic pressure + assert np.isclose( + np.mean(np.array(results["flow_in"][0])), 43.21028819256006, rtol=RTOL_FLOW + ) # aortic inflow + + +def test_closed_loop_heart_with_coronaries(): + results = run_test_case_by_name("closedLoopHeart_withCoronaries") + assert np.isclose( + np.mean(np.array(results["pressure_in"][0])), 50.162313086833805, rtol=RTOL_PRES + ) # mean aortic pressure + assert np.isclose( + np.amax(np.array(results["pressure_in"][0])), 69.01524513715958, rtol=RTOL_PRES + ) # max aortic pressure + assert np.isclose( + np.mean(np.array(results["flow_in"][0])), 38.05442038841015, rtol=RTOL_FLOW + ) # mean aortic flow + assert np.isclose( + np.amax(np.array(results["flow_in"][0])), 171.35198346122127, rtol=RTOL_FLOW + ) # max aortic flow + + +def test_coupled_block_heart_single_vessel(): + result = run_test_case_by_name( + "coupledBlock_closedLoopHeart_singleVessel", output_variable_based=True + ) + aortic_pressure = result[result.name == "pressure:J_heart_outlet:external_outlet"][ + "y" + ].to_numpy() + assert np.isclose( + np.mean(aortic_pressure[-50:]), 69.92379300168665, rtol=RTOL_PRES + ) # mean aortic pressure + assert np.isclose( + np.amax(aortic_pressure[-50:]), 91.44472791507646, rtol=RTOL_PRES + ) # max aortic pressure + assert np.isclose( + np.amin(aortic_pressure[-50:]), 49.246695924657494, rtol=RTOL_PRES + ) # min aortic pressure + + +def test_coupled_block_heart_with_coronaries(): + result = run_test_case_by_name( + "coupledBlock_closedLoopHeart_withCoronaries", output_variable_based=True + ) + aortic_pressure = result[result.name == "pressure:J_heart_outlet:external_outlet"][ + "y" + ].to_numpy() + assert np.isclose( + np.mean(aortic_pressure[-50:]), 59.52487958523876, rtol=RTOL_PRES + ) # mean aortic pressure + assert np.isclose( + np.amax(aortic_pressure[-50:]), 81.0040824877808, rtol=RTOL_PRES + ) # max aortic pressure + assert np.isclose( + np.amin(aortic_pressure[-50:]), 38.80066561075395, rtol=RTOL_PRES + ) # min aortic pressure + + +def test_steady_flow_calibration(): + with open( + os.path.join(this_file_dir, "cases", "steadyFlow_calibration.json") + ) as ff: + config = json.load(ff) + + result = svzerodplus.calibrate(config) + + calibrated_parameters = result["vessels"][0]["zero_d_element_values"] + + assert np.isclose( + np.mean(calibrated_parameters["R_poiseuille"]), 100, rtol=RTOL_PRES + ) + assert np.isclose(np.mean(calibrated_parameters["C"]), 0.0001, rtol=RTOL_PRES) + assert np.isclose(np.mean(calibrated_parameters["L"]), 1.0, rtol=RTOL_PRES) + assert np.isclose( + np.mean(calibrated_parameters["stenosis_coefficient"]), 0.0, rtol=RTOL_PRES + ) + + +@pytest.mark.parametrize("model_id", ["0080_0001", "0104_0001", "0140_2001"]) +def test_calibration_vmr(model_id): + """Test actual models from the vascular model repository.""" + with open( + os.path.join( + this_file_dir, "cases", "vmr", "input", f"{model_id}_calibrate_from_0d.json" + ) + ) as ff: + config = json.load(ff) + + with open( + os.path.join( + this_file_dir, + "cases", + "vmr", + "reference", + f"{model_id}_optimal_from_0d.json", + ) + ) as ff: + reference = json.load(ff) + + result = svzerodplus.calibrate(config) + + for i, vessel in enumerate(reference["vessels"]): + for key, value in vessel["zero_d_element_values"].items(): + np.isclose( + result["vessels"][i]["zero_d_element_values"][key], + value, + rtol=RTOL_PRES, + ) + + for i, junction in enumerate(reference["junctions"]): + if "junction_values" in junction: + for key, value in junction["junction_values"].items(): + np.allclose( + result["junctions"][i]["junction_values"][key], + value, + rtol=RTOL_PRES, + ) diff --git a/tests/test_integration_cpp.py b/tests/test_integration_cpp.py deleted file mode 100644 index e10bf8637..000000000 --- a/tests/test_integration_cpp.py +++ /dev/null @@ -1,550 +0,0 @@ -import json -import os -from io import StringIO -from tempfile import TemporaryDirectory - -import numpy as np -import pandas as pd -import pytest - -import svzerodplus - -this_file_dir = os.path.abspath(os.path.dirname(__file__)) - -RTOL_PRES = 1.0e-7 -RTOL_FLOW = 1.0e-8 - - -def run_test_case_by_name(name, output_variable_based=False, folder="."): - """Run a test case by its case name. - - Args: - name: Name of the test case. - testdir: Directory for performing the simulation. - """ - testfile = os.path.join(this_file_dir, "cases", folder, name + ".json") - with open(testfile) as ff: - config = json.load(ff) - output = svzerodplus.simulate(config) - result = pd.read_csv(StringIO(output)) - - if output_variable_based == False: - output = { - "pressure_in": {}, - "pressure_out": {}, - "flow_in": {}, - "flow_out": {}, - } - - last_seg_id = 0 - - for vessel in config["vessels"]: - name = vessel["vessel_name"] - branch_id, seg_id = name.split("_") - branch_id, seg_id = int(branch_id[6:]), int(seg_id[3:]) - vessel_id = vessel["vessel_id"] - - if seg_id == 0: - output["pressure_in"][branch_id] = np.array( - result[result.name == name]["pressure_in"] - ) - output["flow_in"][branch_id] = np.array( - result[result.name == name]["flow_in"] - ) - output["pressure_out"][branch_id] = np.array( - result[result.name == name]["pressure_out"] - ) - output["flow_out"][branch_id] = np.array( - result[result.name == name]["flow_out"] - ) - elif seg_id > last_seg_id: - output["pressure_out"][branch_id] = np.array( - result[result.name == name]["pressure_out"] - ) - output["flow_out"][branch_id] = np.array( - result[result.name == name]["flow_out"] - ) - - last_seg_id = seg_id - - elif output_variable_based == True: - output = result - - return output - - -def get_result(result_array, field, branch, time_step): - """ "Get results at specific field, branch, branch_node and time step.""" - # extract result - return result_array[field][branch][time_step] - - -def test_steady_flow_R_R(): - results = run_test_case_by_name("steadyFlow_R_R") - assert np.isclose( - get_result(results, "pressure_in", 0, -1), 1100.0, rtol=RTOL_PRES - ) # inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 0, -1), 600.0, rtol=RTOL_PRES - ) # outlet pressure - assert np.isclose( - get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW - ) # inlet flow - assert np.isclose( - get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW - ) # outlet flow - - -def test_steady_flow_r_coronary(): - results = run_test_case_by_name("steadyFlow_R_coronary") - assert np.isclose( - get_result(results, "pressure_in", 0, -1), 2000.0, rtol=RTOL_PRES - ) # inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 0, -1), 1500.0, rtol=RTOL_PRES - ) # outlet pressure - assert np.isclose( - get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW - ) # inlet flow - assert np.isclose( - get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW - ) # outlet flow - - -def test_steady_flow_rlc_r(): - results = run_test_case_by_name("steadyFlow_RLC_R") - assert np.isclose( - get_result(results, "pressure_in", 0, -1), 1100.0, rtol=RTOL_PRES - ) # inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 0, -1), 600.0, rtol=RTOL_PRES - ) # outlet pressure - assert np.isclose( - get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW - ) # inlet flow - assert np.isclose( - get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW - ) # outlet flow - - -def test_steady_flow_rc_r(): - results = run_test_case_by_name("steadyFlow_RC_R") - assert np.isclose( - get_result(results, "pressure_in", 0, -1), 1100.0, rtol=RTOL_PRES - ) # inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 0, -1), 600.0, rtol=RTOL_PRES - ) # outlet pressure - assert np.isclose( - get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW - ) # inlet flow - assert np.isclose( - get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW - ) # outlet flow - - -def test_steady_flow_rl_r(): - results = run_test_case_by_name("steadyFlow_RL_R") - assert np.isclose( - get_result(results, "pressure_in", 0, -1), 1100.0, rtol=RTOL_PRES - ) # inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 0, -1), 600.0, rtol=RTOL_PRES - ) # outlet pressure - assert np.isclose( - get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW - ) # inlet flow - assert np.isclose( - get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW - ) # outlet flow - - -def test_steady_flow_r_rcr(): - results = run_test_case_by_name("steadyFlow_R_RCR") - assert np.isclose( - get_result(results, "pressure_in", 0, -1), 10500.0, rtol=RTOL_PRES - ) # inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 0, -1), 10000.0, rtol=RTOL_PRES - ) # outlet pressure - assert np.isclose( - get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW - ) # inlet flow - assert np.isclose( - get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW - ) # outlet flow - - -def test_steady_flow_r_steady_pressure(): - results = run_test_case_by_name("steadyFlow_R_steadyPressure") - assert np.isclose( - get_result(results, "pressure_in", 0, -1), 1500.0, rtol=RTOL_PRES - ) # inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 0, -1), 1000.0, rtol=RTOL_PRES - ) # outlet pressure - assert np.isclose( - get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW - ) # inlet flow - assert np.isclose( - get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW - ) # outlet flow - - -def test_steady_flow_stenosis_r(): - results = run_test_case_by_name("steadyFlow_stenosis_R") - assert np.isclose( - get_result(results, "pressure_in", 0, -1), 3600.0, rtol=RTOL_PRES - ) # inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 0, -1), 600.0, rtol=RTOL_PRES - ) # outlet pressure - assert np.isclose( - get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW - ) # inlet flow - assert np.isclose( - get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW - ) # outlet flow - - -def test_steady_flow_bifurcationr_r1(): - results = run_test_case_by_name("steadyFlow_bifurcationR_R1") - assert np.isclose( - get_result(results, "pressure_in", 0, -1), 1100.0, rtol=RTOL_PRES - ) # parent inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 0, -1), 600.0, rtol=RTOL_PRES - ) # parent outlet pressure - assert np.isclose( - get_result(results, "pressure_in", 1, -1), 600.0, rtol=RTOL_PRES - ) # daughter1 inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 1, -1), 350.0, rtol=RTOL_PRES - ) # daughter1 outlet pressure - assert np.isclose( - get_result(results, "pressure_in", 2, -1), 600.0, rtol=RTOL_PRES - ) # daughter2 inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 2, -1), 350.0, rtol=RTOL_PRES - ) # daughter2 outlet pressure - assert np.isclose( - get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW - ) # parent inlet flow - assert np.isclose( - get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW - ) # parent outlet flow - assert np.isclose( - get_result(results, "flow_in", 1, -1), 2.5, rtol=RTOL_FLOW - ) # daughter1 inlet flow - assert np.isclose( - get_result(results, "flow_out", 1, -1), 2.5, rtol=RTOL_FLOW - ) # daughter1 outlet flow - assert np.isclose( - get_result(results, "flow_in", 2, -1), 2.5, rtol=RTOL_FLOW - ) # daughter2 inlet flow - assert np.isclose( - get_result(results, "flow_out", 2, -1), 2.5, rtol=RTOL_FLOW - ) # daughter2 outlet flow - - -def test_steady_flow_bifurcationr_r2(): - results = run_test_case_by_name("steadyFlow_bifurcationR_R2") - assert np.isclose( - get_result(results, "pressure_in", 0, -1), 3462.5, rtol=RTOL_PRES - ) # parent inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 0, -1), 1962.5, rtol=RTOL_PRES - ) # parent outlet pressure - assert np.isclose( - get_result(results, "pressure_in", 1, -1), 1962.5, rtol=RTOL_PRES - ) # daughter1 inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 1, -1), 432.5, rtol=RTOL_PRES - ) # daughter1 outlet pressure - assert np.isclose( - get_result(results, "pressure_in", 2, -1), 1962.5, rtol=RTOL_PRES - ) # daughter2 inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 2, -1), 1375.0, rtol=RTOL_PRES - ) # daughter2 outlet pressure - assert np.isclose( - get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW - ) # parent inlet flow - assert np.isclose( - get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW - ) # parent outlet flow - assert np.isclose( - get_result(results, "flow_in", 1, -1), 3.825, rtol=RTOL_FLOW - ) # daughter1 inlet flow - assert np.isclose( - get_result(results, "flow_out", 1, -1), 3.825, rtol=RTOL_FLOW - ) # daughter1 outlet flow - assert np.isclose( - get_result(results, "flow_in", 2, -1), 1.175, rtol=RTOL_FLOW - ) # daughter2 inlet flow - assert np.isclose( - get_result(results, "flow_out", 2, -1), 1.175, rtol=RTOL_FLOW - ) # daughter2 outlet flow - - -def test_steadyFlow_blood_vessel_junction(): - results = run_test_case_by_name("steadyFlow_blood_vessel_junction") - assert np.isclose( - get_result(results, "pressure_in", 0, -1), 3462.5, rtol=RTOL_PRES - ) # parent inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 0, -1), 1962.5, rtol=RTOL_PRES - ) # parent outlet pressure - assert np.isclose( - get_result(results, "pressure_in", 1, -1), 1580.0, rtol=RTOL_PRES - ) # daughter1 inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 1, -1), 432.5, rtol=RTOL_PRES - ) # daughter1 outlet pressure - assert np.isclose( - get_result(results, "pressure_in", 2, -1), 1727.5, rtol=RTOL_PRES - ) # daughter2 inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 2, -1), 1375.0, rtol=RTOL_PRES - ) # daughter2 outlet pressure - assert np.isclose( - get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW - ) # parent inlet flow - assert np.isclose( - get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW - ) # parent outlet flow - assert np.isclose( - get_result(results, "flow_in", 1, -1), 3.825, rtol=RTOL_FLOW - ) # daughter1 inlet flow - assert np.isclose( - get_result(results, "flow_out", 1, -1), 3.825, rtol=RTOL_FLOW - ) # daughter1 outlet flow - assert np.isclose( - get_result(results, "flow_in", 2, -1), 1.175, rtol=RTOL_FLOW - ) # daughter2 inlet flow - assert np.isclose( - get_result(results, "flow_out", 2, -1), 1.175, rtol=RTOL_FLOW - ) # daughter2 outlet flow - - -def test_pulsatile_flow_r_rcr(): - results = run_test_case_by_name("pulsatileFlow_R_RCR") - assert np.isclose( - get_result(results, "pressure_in", 0, 0), 4620.0, rtol=RTOL_PRES - ) # inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 0, 0), 4400.0, rtol=RTOL_PRES - ) # outlet pressure - assert np.isclose( - get_result(results, "flow_in", 0, 0), 2.2, rtol=RTOL_FLOW - ) # inlet flow - assert np.isclose( - get_result(results, "flow_out", 0, 0), 2.2, rtol=RTOL_FLOW - ) # outlet flow - - -def test_pulsatile_flow_r_coronary(): - results = run_test_case_by_name("pulsatileFlow_R_coronary") - assert np.isclose( - get_result(results, "pressure_in", 0, 0), 880.0, rtol=RTOL_PRES - ) # inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 0, 0), 660.0, rtol=RTOL_PRES - ) # outlet pressure - assert np.isclose( - get_result(results, "flow_in", 0, 0), 2.2, rtol=RTOL_FLOW - ) # inlet flow - assert np.isclose( - get_result(results, "flow_out", 0, 0), 2.2, rtol=RTOL_FLOW - ) # outlet flow - - -def test_pulsatile_flow_cstenosis_steady_pressure(): - results = run_test_case_by_name("pulsatileFlow_CStenosis_steadyPressure") - assert np.isclose( - get_result(results, "pressure_in", 0, -439), - 0.5933049197138334, - rtol=1.0e-5, - ) # inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 0, -439), 0.1, rtol=1.0e-5 - ) # outlet pressure - assert np.isclose( - get_result(results, "flow_in", 0, -439), - 0.7023611813029965, - rtol=1.0e-5, - ) # inlet flow - assert np.isclose( - get_result(results, "flow_out", 0, -439), - 0.7018707542098627, - rtol=1.0e-5, - ) # outlet flow - - -def test_steady_flow_confluencer_r(): - results = run_test_case_by_name("steadyFlow_confluenceR_R") - assert np.isclose( - get_result(results, "pressure_in", 0, -1), 6600.0, rtol=RTOL_PRES - ) # parent inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 0, -1), 6100.0, rtol=RTOL_PRES - ) # parent outlet pressure - assert np.isclose( - get_result(results, "pressure_in", 1, -1), 8100.0, rtol=RTOL_PRES - ) # daughter1 inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 1, -1), 6100.0, rtol=RTOL_PRES - ) # daughter1 outlet pressure - assert np.isclose( - get_result(results, "pressure_in", 2, -1), 6100.0, rtol=RTOL_PRES - ) # daughter2 inlet pressure - assert np.isclose( - get_result(results, "pressure_out", 2, -1), 1600.0, rtol=RTOL_PRES - ) # daughter2 outlet pressure - assert np.isclose( - get_result(results, "flow_in", 0, -1), 5.0, rtol=RTOL_FLOW - ) # parent inlet flow - assert np.isclose( - get_result(results, "flow_out", 0, -1), 5.0, rtol=RTOL_FLOW - ) # parent outlet flow - assert np.isclose( - get_result(results, "flow_in", 1, -1), 10.0, rtol=RTOL_FLOW - ) # daughter1 inlet flow - assert np.isclose( - get_result(results, "flow_out", 1, -1), 10.0, rtol=RTOL_FLOW - ) # daughter1 outlet flow - assert np.isclose( - get_result(results, "flow_in", 2, -1), 15.0, rtol=RTOL_FLOW - ) # daughter2 inlet flow - assert np.isclose( - get_result(results, "flow_out", 2, -1), 15.0, rtol=RTOL_FLOW - ) # daughter2 outlet flow - - -def test_closed_loop_heart_single_vessel(tmpdir): - results = run_test_case_by_name("closedLoopHeart_singleVessel") - assert np.isclose( - np.mean(np.array(results["pressure_in"][0])), 55.703345704742844, rtol=RTOL_PRES - ) # mean aortic pressure - assert np.isclose( - np.amax(np.array(results["pressure_in"][0])), 73.97450170686889, rtol=RTOL_PRES - ) # max aortic pressure - assert np.isclose( - np.amin(np.array(results["pressure_in"][0])), 0.0, rtol=RTOL_PRES - ) # min aortic pressure - assert np.isclose( - np.mean(np.array(results["flow_in"][0])), 43.21028819256006, rtol=RTOL_FLOW - ) # aortic inflow - - -def test_closed_loop_heart_with_coronaries(tmpdir): - results = run_test_case_by_name("closedLoopHeart_withCoronaries") - assert np.isclose( - np.mean(np.array(results["pressure_in"][0])), 50.162313086833805, rtol=RTOL_PRES - ) # mean aortic pressure - assert np.isclose( - np.amax(np.array(results["pressure_in"][0])), 69.01524513715958, rtol=RTOL_PRES - ) # max aortic pressure - assert np.isclose( - np.mean(np.array(results["flow_in"][0])), 38.05442038841015, rtol=RTOL_FLOW - ) # mean aortic flow - assert np.isclose( - np.amax(np.array(results["flow_in"][0])), 171.35198346122127, rtol=RTOL_FLOW - ) # max aortic flow - - -def test_coupled_block_heart_single_vessel(tmpdir): - result = run_test_case_by_name( - "coupledBlock_closedLoopHeart_singleVessel", output_variable_based=True - ) - aortic_pressure = result[result.name == "pressure:J_heart_outlet:external_outlet"][ - "y" - ].to_numpy() - assert np.isclose( - np.mean(aortic_pressure[-50:]), 69.92379300168665, rtol=RTOL_PRES - ) # mean aortic pressure - assert np.isclose( - np.amax(aortic_pressure[-50:]), 91.44472791507646, rtol=RTOL_PRES - ) # max aortic pressure - assert np.isclose( - np.amin(aortic_pressure[-50:]), 49.246695924657494, rtol=RTOL_PRES - ) # min aortic pressure - - -def test_coupled_block_heart_with_coronaries(tmpdir): - result = run_test_case_by_name( - "coupledBlock_closedLoopHeart_withCoronaries", output_variable_based=True - ) - aortic_pressure = result[result.name == "pressure:J_heart_outlet:external_outlet"][ - "y" - ].to_numpy() - assert np.isclose( - np.mean(aortic_pressure[-50:]), 59.52487958523876, rtol=RTOL_PRES - ) # mean aortic pressure - assert np.isclose( - np.amax(aortic_pressure[-50:]), 81.0040824877808, rtol=RTOL_PRES - ) # max aortic pressure - assert np.isclose( - np.amin(aortic_pressure[-50:]), 38.80066561075395, rtol=RTOL_PRES - ) # min aortic pressure - - -def test_steady_flow_calibration(tmpdir): - with open( - os.path.join(this_file_dir, "cases", "steadyFlow_calibration.json") - ) as ff: - config = json.load(ff) - - result = svzerodplus.calibrate(config) - - calibrated_parameters = result["vessels"][0]["zero_d_element_values"] - - assert np.isclose( - np.mean(calibrated_parameters["R_poiseuille"]), 100, rtol=RTOL_PRES - ) - assert np.isclose(np.mean(calibrated_parameters["C"]), 0.0001, rtol=RTOL_PRES) - assert np.isclose(np.mean(calibrated_parameters["L"]), 1.0, rtol=RTOL_PRES) - assert np.isclose( - np.mean(calibrated_parameters["stenosis_coefficient"]), 0.0, rtol=RTOL_PRES - ) - - -@pytest.mark.parametrize("model_id", ["0080_0001", "0104_0001", "0140_2001"]) -def test_calibration_vmr(model_id): - """Test actual models from the vascular model repository.""" - with open( - os.path.join( - this_file_dir, "cases", "vmr", "input", f"{model_id}_calibrate_from_0d.json" - ) - ) as ff: - config = json.load(ff) - - with open( - os.path.join( - this_file_dir, - "cases", - "vmr", - "reference", - f"{model_id}_optimal_from_0d.json", - ) - ) as ff: - reference = json.load(ff) - - result = svzerodplus.calibrate(config) - - for i, vessel in enumerate(reference["vessels"]): - for key, value in vessel["zero_d_element_values"].items(): - np.isclose( - result["vessels"][i]["zero_d_element_values"][key], - value, - rtol=RTOL_PRES, - ) - - for i, junction in enumerate(reference["junctions"]): - if "junction_values" in junction: - for key, value in junction["junction_values"].items(): - np.allclose( - result["junctions"][i]["junction_values"][key], - value, - rtol=RTOL_PRES, - )