diff --git a/.github/workflows/cd.yml b/.github/workflows/cd.yml index 45846e9cb..bea75c0c2 100644 --- a/.github/workflows/cd.yml +++ b/.github/workflows/cd.yml @@ -1,6 +1,14 @@ name: CD -on: [push, workflow_dispatch] +on: + pull_request: + push: + branches: + - master + release: + types: + - published + env: REGISTRY: ghcr.io diff --git a/.github/workflows/pip.yml b/.github/workflows/pip.yml index ebca15d74..cf8de5c9d 100644 --- a/.github/workflows/pip.yml +++ b/.github/workflows/pip.yml @@ -6,6 +6,9 @@ on: push: branches: - master + release: + types: + - published jobs: build: @@ -14,9 +17,14 @@ jobs: strategy: fail-fast: false matrix: - platform: [windows-latest, macos-13, macos-14, ubuntu-latest] -# python-version: ["3.7", "3.12", "pypy-3.9"] - python-version: ["3.12"] + #platform: [ windows-latest, macos-13, macos-14, ubuntu-latest ] + platform: [ macos-13, macos-14, ubuntu-latest ] + python-version: ["3.9","3.10","3.11","3.12"] +# exclude: +# - platform: windows-latest +# python-version: "3.9" +# - platform: windows-latest +# python-version: "3.12" steps: - uses: actions/checkout@v4 @@ -79,49 +87,94 @@ jobs: mingw-w64-clang-x86_64-libaec mingw-w64-clang-x86_64-python-pip-tools - - - name: Install Linux Dependencies - if: matrix.platform == 'ubuntu-latest' - run: | - sudo apt-get update - sudo apt-get install -y libboost-all-dev - sudo apt-get install -y libhdf5-dev - sudo apt-get install -y ninja-build - - gcc --version - gfortran --version - cmake --version - dpkg -s libboost-all-dev - dpkg -s libhdf5-dev - - - uses: actions/setup-python@v5 with: python-version: "${{ matrix.python-version }}" - name: Build and install macos-13 run: | + export MACOSX_DEPLOYMENT_TARGET=13.0 PATH="/usr/local/opt/llvm/bin:$PATH" - pip install --verbose . + pip wheel . -w dist + ls -al dist + pip install delocate + delocate-listdeps dist/*.whl + delocate-wheel -v dist/*.whl + pip install dist/*.whl if: matrix.platform == 'macos-13' - name: Build and install macos-14 run: | PATH="/opt/homebrew/opt/llvm/bin:$PATH" - pip install --verbose . + pip wheel . -w dist + ls -al dist + pip install delocate + delocate-listdeps dist/*.whl + delocate-wheel -v dist/*.whl + pip install dist/*.whl if: matrix.platform == 'macos-14' - name: Build and install windows-latest shell: msys2 {0} run: | - pip install --verbose -Ccmake.define.CMAKE_C_COMPILER="clang.exe" -Ccmake.define.CMAKE_CXX_COMPILER="clang++.exe" . + pip wheel --verbose -Ccmake.define.CMAKE_C_COMPILER="clang.exe" -Ccmake.define.CMAKE_CXX_COMPILER="clang++.exe" . -w dist + pip install delvewheel + for wheel_file in dist/*.whl; do + echo "not inspecting $wheel_file" + echo "delvewheel needed -v $wheel_file" + echo "not repairing $wheel_file" + echo "delvewheel repair -v $wheel_file" + done + pip install dist/*.whl if: matrix.platform == 'windows-latest' - - name: Build and install ubuntu-latest + - name: Build and install manylinux for python 3.9 + run: | + docker run --rm -v $(pwd):/io quay.io/pypa/manylinux2014_x86_64 \ + /bin/bash -c \ + "yum install -y boost-devel hdf5-devel ninja-build && \ + /opt/python/cp39-cp39/bin/pip wheel /io/ -w /io/dist && \ + auditwheel repair /io/dist/*.whl -w /io/dist/" + echo "keeping only the manylinux wheels, remove those with -linux_x86_64.whl in the name" + sudo rm -rf dist/*-linux_x86_64.whl + pip install dist/*.whl + if: matrix.platform == 'ubuntu-latest' && matrix.python-version == '3.9' + + - name: Build and install manylinux for python 3.10 + run: | + docker run --rm -v $(pwd):/io quay.io/pypa/manylinux2014_x86_64 \ + /bin/bash -c \ + "yum install -y boost-devel hdf5-devel ninja-build && \ + /opt/python/cp310-cp310/bin/pip wheel /io/ -w /io/dist && \ + auditwheel repair /io/dist/*.whl -w /io/dist/" + echo "keeping only the manylinux wheels, remove those with -linux_x86_64.whl in the name" + sudo rm -rf dist/*-linux_x86_64.whl + pip install dist/*.whl + if: matrix.platform == 'ubuntu-latest' && matrix.python-version == '3.10' + + - name: Build and install manylinux for python 3.11 + run: | + docker run --rm -v $(pwd):/io quay.io/pypa/manylinux2014_x86_64 \ + /bin/bash -c \ + "yum install -y boost-devel hdf5-devel ninja-build && \ + /opt/python/cp311-cp311/bin/pip wheel /io/ -w /io/dist && \ + auditwheel repair /io/dist/*.whl -w /io/dist/" + echo "keeping only the manylinux wheels, remove those with -linux_x86_64.whl in the name" + sudo rm -rf dist/*-linux_x86_64.whl + pip install dist/*.whl + if: matrix.platform == 'ubuntu-latest' && matrix.python-version == '3.11' + + - name: Build and install manylinux for python 3.12 run: | - platform=linux - pip install --verbose . - if: matrix.platform == 'ubuntu-latest' + docker run --rm -v $(pwd):/io quay.io/pypa/manylinux2014_x86_64 \ + /bin/bash -c \ + "yum install -y boost-devel hdf5-devel ninja-build && \ + /opt/python/cp312-cp312/bin/pip wheel /io/ -w /io/dist && \ + auditwheel repair /io/dist/*.whl -w /io/dist/" + echo "keeping only the manylinux wheels, remove those with -linux_x86_64.whl in the name" + sudo rm -rf dist/*-linux_x86_64.whl + pip install dist/*.whl + if: matrix.platform == 'ubuntu-latest' && matrix.python-version == '3.12' - name: Test not windows run: | @@ -136,8 +189,75 @@ jobs: pytest if: matrix.platform == 'windows-latest' + - name: set wheel_file and artifact_name environment vars (non-windows) + if: matrix.platform != 'windows-latest' + run: | + for wheel_file in dist/*.whl; do + echo "Uploading $wheel_file" + echo "wheel_file=$wheel_file" >> $GITHUB_ENV + echo "artifact_name=$(basename $wheel_file)" >> $GITHUB_ENV + done + + - name: set wheel_file and artifact_name environment vars (windows) + if: matrix.platform == 'windows-latest' + shell: pwsh + run: | + foreach ($wheel_file in Get-ChildItem -Path dist -Filter *.whl) { + Write-Output "Uploading $wheel_file" + echo "wheel_file=$wheel_file" | Out-File -Append -FilePath $env:GITHUB_ENV -Encoding utf8 + echo "artifact_name=$(Split-Path $wheel_file -Leaf)" | Out-File -Append -FilePath $env:GITHUB_ENV -Encoding utf8 + } + + - name: Upload wheel file + uses: actions/upload-artifact@v4 + with: + name: ${{ env.artifact_name }} + path: ${{ env.wheel_file }} + - name: Setup tmate session uses: mxschmitt/action-tmate@v3 if: ${{ failure() }} with: github-token: ${{ secrets.GITHUB_TOKEN }} + + publish-to-pypi: + name: Publish to PyPI + needs: build + runs-on: ubuntu-latest + if: github.event_name == 'release' && github.event.action == 'published' + environment: pypi + permissions: + id-token: write + + steps: + - name: download all artifacts (wheels and executables) + uses: actions/download-artifact@v4 + with: + path: dist_temp/ + + - name: Install Python + uses: actions/setup-python@v2 + with: + python-version: '3.12' + + - name: Publish wheels to PyPI + env: + TWINE_USERNAME: __token__ + TWINE_PASSWORD: ${{ secrets.PYPI_TOKEN }} + run: | + pip install twine + mkdir dist + echo "current directory is $(pwd)" + ls -al dist* + echo "absolute paths of wheels:" + find dist_temp -type f -name '*.whl' -exec sh -c 'echo "$(realpath {})"/' \; + echo "moving wheels to dist directory" + find dist_temp -type f -name '*.whl' -exec sh -c 'mv "$(realpath {})" dist/' \; + twine upload --repository pypi dist/*.whl + + - name: Setup tmate session for PyPI upload failure + uses: mxschmitt/action-tmate@v3 + if: ${{ failure() }} + with: + github-token: ${{ secrets.GITHUB_TOKEN }} + diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml index 327dcb707..cc2b606e2 100644 --- a/.github/workflows/wheels.yml +++ b/.github/workflows/wheels.yml @@ -1,4 +1,4 @@ -name: Wheels +name: "Wheels" on: workflow_dispatch: @@ -107,7 +107,7 @@ jobs: mingw-w64-clang-x86_64-libzip mingw-w64-clang-x86_64-zlib mingw-w64-clang-x86_64-libaec - + mingw-w64-clang-x86_64-python-pip-tools - name: Install Linux Dependencies if: matrix.os == 'ubuntu-latest' @@ -123,9 +123,40 @@ jobs: dpkg -s libboost-all-dev dpkg -s libhdf5-dev + # Used to host cibuildwheel + - uses: actions/setup-python@v5 + + - name: Install cibuildwheel + run: python -m pip install cibuildwheel==2.18.1 + + - name: build wheel on macos-13 + if: matrix.os == 'macos-13' + run: | + export PATH="/usr/local/opt/llvm/bin:$PATH" + python -m cibuildwheel --output-dir wheelhouse + env: + CIBW_ARCHS: native + + - name: build wheel on macos-14 + if: matrix.os == 'macos-14' + run: | + export PATH="/opt/homebrew/opt/llvm/bin:$PATH" + python -m cibuildwheel --output-dir wheelhouse + env: + CIBW_ARCHS: native + + - name: build wheel on windows-latest + if: matrix.os == 'windows-latest' + shell: msys2 {0} + run: | + python -m cibuildwheel --output-dir wheelhouse + env: + CIBW_ARCHS: native - name: build wheel - uses: pypa/cibuildwheel@v2.17 + if: matrix.os == 'ubuntu-latest' + run: | + python -m cibuildwheel --output-dir wheelhouse env: CIBW_ARCHS: native diff --git a/.gitignore b/.gitignore index 2c0fede2f..b9c043b76 100644 --- a/.gitignore +++ b/.gitignore @@ -17,5 +17,12 @@ /.idea/ /nfsim/ /all_solvers/* +/.venv/ +/venv/ .DS_Store +/dist/ + +tests/__pycache__/ + +*.whl diff --git a/CMakeLists.txt b/CMakeLists.txt index 2e8b5b683..bb050d33e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,6 +3,7 @@ option(OPTION_TARGET_PYTHON_BINDING "build Python bindings" ON) option(OPTION_TARGET_MESSAGING "Messaging (requires libcurl)" off) option(OPTION_TARGET_DOCS "Generate Doxygen documentation" off) option(OPTION_TARGET_FV_SOLVER on) +option(OPTION_TARGET_TESTS "Build tests" off) option(OPTION_VCELL "Compile Smoldyn for VCell" ON) option(OPTION_NSV "Compile Smoldyn with NextSubvolume functionality" OFF) option(OPTION_USE_OPENGL "Build with OpenGL support" OFF) @@ -217,7 +218,7 @@ if (OPTION_TARGET_PYTHON_BINDING) find_package(pybind11 CONFIG REQUIRED) # Add a library using FindPython's tooling (pybind11 also provides a helper like this) - python_add_library(_core MODULE src/main.cpp src/SolverMain.cpp src/SolverMain.h WITH_SOABI) + python_add_library(_core MODULE src/main.cpp WITH_SOABI) target_link_libraries(_core PRIVATE vcell pybind11::headers) # This is passing in the version as a define just as an example @@ -226,35 +227,4 @@ if (OPTION_TARGET_PYTHON_BINDING) # The install directory is the output (wheel) directory install(TARGETS _core DESTINATION pyvcell_fvsolver) -else(OPTION_TARGET_PYTHON_BINDING) - - include(FetchContent) - FetchContent_Declare( - googletest - URL https://github.com/google/googletest/archive/03597a01ee50ed33e9dfd640b249b4be3799d395.zip - ) - if (WINDOWS) - # For Windows: Prevent overriding the parent project's compiler/linker settings - set(gtest_force_shared_crt ON CACHE BOOL "" FORCE) - endif() - - FetchContent_MakeAvailable(googletest) - enable_testing() - - get_cmake_property(_variableNames VARIABLES) - list (SORT _variableNames) - foreach (_variableName ${_variableNames}) - message(STATUS "${_variableName}=${${_variableName}}") - endforeach() - - include(CMakePrintHelpers) - cmake_print_variables(OPTION_TARGET_MESSAGING - OPTION_TARGET_DOCS - OPTION_TARGET_FV_SOLVER - OPTION_TARGET_LIBSMOLDYN OPTION_TARGET_VCELL OPTION_TARGET_NSV ) - cmake_print_variables(CMAKE_CXX_FLAGS CMAKE_C_FLAGS CMAKE_Fortran_FLAGS) - cmake_print_variables(CMAKE_SYSTEM_NAME WINDOWS WIN32 MINGW APPLE ARCH_64bit ARCH_32bit) - cmake_print_variables(CMAKE_CPP_COMPILER CMAKE_C_COMPILER CMAKE_CXX_COMPILER CMAKE_Fortran_COMPILER) - cmake_print_variables(HAVE_ZLIB) - endif (OPTION_TARGET_PYTHON_BINDING) \ No newline at end of file diff --git a/README.md b/README.md index 337cf1b1b..c5dc82680 100644 --- a/README.md +++ b/README.md @@ -8,9 +8,11 @@ ![CI](https://github.com/virtualcell/vcell-solvers/actions/workflows/cd.yml/badge.svg) # vcell-fvsolver -Virtual Cell solvers [virtualcell/vcell-solvers](https://github.com/virtualcell/vcell-solvers) is a collection of numerical simulation codes used in the Virtual Cell framework [virtualcell/vcell](https://github.com/virtualcell/vcell)). +Virtual Cell Finite Volume solver [virtualcell/vcell-fvsolver](https://github.com/virtualcell/vcell-fvsolver) +is a reaction-diffusion-advection PDE solver for computational cell biology. +This solver is used within the Virtual Cell modeling and simulation application [virtualcell/vcell](https://github.com/virtualcell/vcell) +and as a component in the Virtual Cell Python API [virtualcell/pyvcell](https://github.com/virtualcell/pyvcell) (coming soon). - ## The Virtual Cell Project The Virtual Cell is a modeling and simulation framework for computational biology. For details see http://vcell.org and http://github.com/virtualcell. @@ -27,4 +29,4 @@ and generates the output files (.log, .zip, .mesh, .meshmetrics, .hdf5). The .functions file is not used by the solver, but is helpful for interpreting the results in the context of the original model. -This package is intented to be used by the Virtual Cell Python API, pyvcell (coming soon). +This package is intended to be used by the Virtual Cell Python API [virtualcell/pyvcell](https://github.com/virtualcell/pyvcell) (coming soon). diff --git a/VCell/CMakeLists.txt b/VCell/CMakeLists.txt index c038c29bb..fbd2bab89 100644 --- a/VCell/CMakeLists.txt +++ b/VCell/CMakeLists.txt @@ -62,6 +62,7 @@ set(HEADER_FILES include/VCELL/SimulationExpression.h # include/VCELL/SimulationMessaging.h include/VCELL/Solver.h + include/VCELL/SolverMain.h include/VCELL/SparseLinearSolver.h # include/VCELL/SparseMatrix.h include/VCELL/SparseMatrixEqnBuilder.h @@ -141,6 +142,7 @@ set(SRC_FILES src/SimulationExpression.cpp # src/SimulationMessaging.cpp src/Solver.cpp + src/SolverMain.cpp src/SparseLinearSolver.cpp # src/SparseMatrix.cpp src/SparseMatrixEqnBuilder.cpp @@ -238,3 +240,8 @@ if (OPTION_TARGET_FV_SOLVER) RUNTIME DESTINATION ${OPTION_EXE_DIRECTORY}) endif() + +# add google tests +if (OPTION_TARGET_TESTS) + add_subdirectory(tests) +endif() \ No newline at end of file diff --git a/VCell/include/VCELL/AlgebraicSystem.h b/VCell/include/VCELL/AlgebraicSystem.h index d41546808..508fe17cb 100644 --- a/VCell/include/VCELL/AlgebraicSystem.h +++ b/VCell/include/VCELL/AlgebraicSystem.h @@ -5,14 +5,17 @@ #ifndef ALGEBRAICSYSTEM_H #define ALGEBRAICSYSTEM_H +class SimTool; + class AlgebraicSystem { public: + virtual ~AlgebraicSystem() = default; void solveSystem(); - virtual void initVars()=0; - inline double getX(int index) {return x[index];} - inline int getDimension() {return dimension;} - inline void setTolerance(double tol){tolerance = tol;} + virtual void initVars(SimTool* sim_tool)=0; + double getX(int index) const {return x[index];} + int getDimension() const {return dimension;} + void setTolerance(double tol) {tolerance = tol;} protected: AlgebraicSystem(int dimension); diff --git a/VCell/include/VCELL/CartesianMesh.h b/VCell/include/VCELL/CartesianMesh.h index 843e1a9f3..39535d64f 100644 --- a/VCell/include/VCELL/CartesianMesh.h +++ b/VCell/include/VCELL/CartesianMesh.h @@ -8,9 +8,8 @@ #include #include #include -#include #include -#include + using std::vector; using std::istream; @@ -18,6 +17,7 @@ class Geometry; class VolumeRegion; class MembraneRegion; class VolumeVariable; +class VCellModel; class Particle; class SparseMatrixPCG; template class IncidenceMatrix; @@ -31,68 +31,65 @@ enum BoundaryLocation {BL_Xm = 0, BL_Xp, BL_Ym, BL_Yp, BL_Zm, BL_Zp}; class CartesianMesh : public Mesh { public: - CartesianMesh(double captureNeighborhood=0); - void initialize(istream& ifs); + explicit CartesianMesh(double captureNeighborhood=0); + void initialize(VCellModel* model, istream& ifs); - virtual WorldCoord getVolumeWorldCoord(long volumeIndex); - virtual WorldCoord getMembraneWorldCoord(long membraneIndex); - virtual WorldCoord getMembraneWorldCoord(MembraneElement *element); + WorldCoord getVolumeWorldCoord(long volumeIndex) override; + WorldCoord getMembraneWorldCoord(long membraneIndex) override; + WorldCoord getMembraneWorldCoord(MembraneElement *element) override; virtual long getVolumeIndex(WorldCoord coord); - virtual double getVolumeOfElement_cu(long volumeIndex); + double getVolumeOfElement_cu(long volumeIndex) override; - virtual void showSummary(FILE *fp); - virtual void write(FILE *fp); - virtual void writeMeshMetrics(FILE* fp); + void showSummary(FILE *fp) override; + void write(FILE *fp) override; + void writeMeshMetrics(FILE* fp) override; - inline double getDomainSizeX() { return domainSizeX; } - inline double getDomainSizeY() { return domainSizeY; } - inline double getDomainSizeZ() { return domainSizeZ; } - inline double getDomainOriginX() { return domainOriginX; } - inline double getDomainOriginY() { return domainOriginY; } - inline double getDomainOriginZ() { return domainOriginZ; } + double getDomainSizeX() const { return domainSizeX; } + double getDomainSizeY() const { return domainSizeY; } + double getDomainSizeZ() const { return domainSizeZ; } + double getDomainOriginX() const { return domainOriginX; } + double getDomainOriginY() const { return domainOriginY; } + double getDomainOriginZ() const { return domainOriginZ; } - inline int getNumVolumeX() { return numX; } - inline int getNumVolumeY() { return numY; } - inline int getNumVolumeZ() { return numZ; } + int getNumVolumeX() const { return numX; } + int getNumVolumeY() const { return numY; } + int getNumVolumeZ() const { return numZ; } - double getXScale_um() { return scaleX_um; } - double getYScale_um() { return scaleY_um; } - double getZScale_um() { return scaleZ_um; } + double getXScale_um() const { return scaleX_um; } + double getYScale_um() const { return scaleY_um; } + double getZScale_um() const { return scaleZ_um; } - double getXArea_squm() { return areaX_squm; } - double getYArea_squm() { return areaY_squm; } - double getZArea_squm() { return areaZ_squm; } + double getXArea_squm() const { return areaX_squm; } + double getYArea_squm() const { return areaY_squm; } + double getZArea_squm() const { return areaZ_squm; } - double getVolume_cu() { return volume_cu; } + double getVolume_cu() const { return volume_cu; } VolumeRegion *getVolumeRegion(int i); MembraneRegion *getMembraneRegion(int i); - int getNumVolumeRegions() { return static_cast(pVolumeRegions.size()); } - int getNumMembraneRegions() { return static_cast(pMembraneRegions.size()); } - + int getNumVolumeRegions() const { return static_cast(pVolumeRegions.size()); } + int getNumMembraneRegions() const { return static_cast(pMembraneRegions.size()); } - int getMembraneNeighborMask(long meindex); - int getMembraneNeighborMask(MembraneElement* element); - double* getMembraneFluxArea(long index); + int getMembraneNeighborMask(long meindex) override; + int getMembraneNeighborMask(MembraneElement* element) override; + double* getMembraneFluxArea(long index) override; - MeshCoord getMeshCoord(long index) { - void asSigned(); + MeshCoord getMeshCoord(long index) const + { MeshCoord mc; - mc.x = index % numX; + mc.x = index % numX; mc.y = (index / numX) % numY; mc.z = index/ numXY; return mc; } -private: - //void setVolumeLists(); - void readGeometryFile(istream& ifs); +private: + void readGeometryFile(VCellModel* model, istream& ifs); void setBoundaryConditions(); void initScale(); - //long getVolumeIndex(MeshCoord); void findMembraneNeighbors(); typedef StatusIndex NeighborIndex; NeighborIndex orthoIndex(long memIndex, long insideIndex, long outsideIndex, long indexer, int boundMask); @@ -100,33 +97,33 @@ class CartesianMesh : public Mesh void findMembranePointInCurve(int n, long index, int neighborDir, int& leftOverN, int& returnNeighbor); - double domainSizeX; - double domainSizeY; - double domainSizeZ; - double domainOriginX; - double domainOriginY; - double domainOriginZ; + double domainSizeX{}; + double domainSizeY{}; + double domainSizeZ{}; + double domainOriginX{}; + double domainOriginY{}; + double domainOriginZ{}; - int numX; - int numY; - int numZ; - int numXY; + int numX{}; + int numY{}; + int numZ{}; + int numXY{}; - double scaleX_um; - double scaleY_um; - double scaleZ_um; + double scaleX_um{}; + double scaleY_um{}; + double scaleZ_um{}; - double areaX_squm; - double areaY_squm; - double areaZ_squm; + double areaX_squm{}; + double areaY_squm{}; + double areaZ_squm{}; - double volume_cu; + double volume_cu{}; struct { int oppositeDirection; int nDirections; - } membraneInfo; + } membraneInfo{}; /** * return index of direction opposite "d"; value depends on number of dimensions */ @@ -139,11 +136,10 @@ class CartesianMesh : public Mesh void writeVolumeElementsMapVolumeRegion(FILE *fp); void writeMembraneRegionMapVolumeRegion(FILE *fp); void writeMembraneElements_Connectivity_Region(FILE *fp); - //void writeContourElements(FILE *fp); vector pVolumeRegions; vector pMembraneRegions; - void computeMembraneCoupling(void); + void computeMembraneCoupling() override; void computeExactNormals(); WorldCoord computeExactNormal(long meIndex); @@ -169,7 +165,7 @@ class CartesianMesh : public Mesh * it is a function of the curvature of the membrane */ //review - ArrayHolder getNormalApproximationHops(const long index); + ArrayHolder getNormalApproximationHops(long index); int computeNormalApproximationHops(int startingIndex, CurvePlane curvePlane, vector curvex, vector curvey, int currentMeIndexInVector, bool bClose); bool findCurve(int startingIndex, CurvePlane curvePlane, vector& curvex, vector& curvey, int& currentMEInVector); diff --git a/VCell/include/VCELL/EllipticVolumeEqnBuilder.h b/VCell/include/VCELL/EllipticVolumeEqnBuilder.h index 7e0a68c5e..05471b3e2 100644 --- a/VCell/include/VCELL/EllipticVolumeEqnBuilder.h +++ b/VCell/include/VCELL/EllipticVolumeEqnBuilder.h @@ -23,11 +23,11 @@ struct VolumeNeighbor; class EllipticVolumeEqnBuilder : public SparseMatrixEqnBuilder { public: - EllipticVolumeEqnBuilder(VolumeVariable *species, CartesianMesh *mesh, int numSolveRegions=0, int* solveRegions=0); - ~EllipticVolumeEqnBuilder(); + EllipticVolumeEqnBuilder(VolumeVariable *species, CartesianMesh *mesh, int numSolveRegions=0, int* solveRegions=nullptr); + ~EllipticVolumeEqnBuilder() override; - void initEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); - void buildEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); + void initEquation(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); + void buildEquation(Simulation *sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); void postProcess(); bool isElliptic() { @@ -57,8 +57,8 @@ class EllipticVolumeEqnBuilder : public SparseMatrixEqnBuilder void init(); void computeLHS(int index, double& Aii, int& numCols, int* columnIndices, double* columnValues, bool& bSort); - double computeRHS(int index); - void preProcess(); + double computeRHS(Simulation* sim, int index); + void preProcess(VCellModel *model); bool checkPeriodicCoupledPairsInRegions(int indexm, int indexp); }; diff --git a/VCell/include/VCELL/EqnBuilder.h b/VCell/include/VCELL/EqnBuilder.h index a508290ec..b2fb4d683 100644 --- a/VCell/include/VCELL/EqnBuilder.h +++ b/VCell/include/VCELL/EqnBuilder.h @@ -6,17 +6,20 @@ #define EQNBUILDER_H class Mesh; +class Simulation; class Variable; +class VCellModel; class EqnBuilder { public: EqnBuilder(Variable *var, Mesh *mesh); + virtual ~EqnBuilder() = default; - virtual void initEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize)=0; - virtual void buildEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize)=0; + virtual void initEquation(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize)=0; + virtual void buildEquation(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize)=0; - Mesh* getMesh() { return mesh; }; + Mesh* getMesh() const { return mesh; }; virtual bool isElliptic() { return false; } protected: diff --git a/VCell/include/VCELL/EqnBuilderReactionForward.h b/VCell/include/VCELL/EqnBuilderReactionForward.h index c2b6dee85..967788d65 100644 --- a/VCell/include/VCELL/EqnBuilderReactionForward.h +++ b/VCell/include/VCELL/EqnBuilderReactionForward.h @@ -16,9 +16,9 @@ class EqnBuilderReactionForward : public EqnBuilder public: EqnBuilderReactionForward(VolumeVariable *species, Mesh *mesh, ODESolver *solver); - void initEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) + void initEquation(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) override { } - void buildEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); + void buildEquation(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) override; private: ODESolver* odeSolver; diff --git a/VCell/include/VCELL/FVDataSet.h b/VCell/include/VCELL/FVDataSet.h index 869a4df30..6c8e66a6c 100644 --- a/VCell/include/VCELL/FVDataSet.h +++ b/VCell/include/VCELL/FVDataSet.h @@ -14,7 +14,7 @@ class FVDataSet public: static void read(const char *filename, Simulation *sim); static void write(const char *filename, SimulationExpression *sim, bool bCompress); - static void convolve(Simulation* sim, Variable* var, double* values); + static void convolve(SimulationExpression* sim, Variable* var, double* values); static void readRandomVariables(char* filename, SimulationExpression* sim); }; diff --git a/VCell/include/VCELL/FVSolver.h b/VCell/include/VCELL/FVSolver.h index 18f4a301c..e23c80774 100644 --- a/VCell/include/VCELL/FVSolver.h +++ b/VCell/include/VCELL/FVSolver.h @@ -2,6 +2,8 @@ #define FVSOLVER_H #include + +#include "VCellModel.h" using std::string; using std::istream; @@ -18,61 +20,55 @@ class FastSystemExpression; class SimulationExpression; class Structure; class Membrane; -//class PdeResultSet; class FVSolver { public: - FVSolver(istream& fvinput, int taskID=-1, const char* outdir=0, bool bSimZip=true); + FVSolver(const char* outdir=nullptr); virtual ~FVSolver(); - void createSimTool(istream& ifsInput, int taskID); - void solve(bool bLoadFinal=true, double* paramValues=0); + SimTool* createSimTool(istream& ifsInput, istream& vcgInput, int taskID, bool bSimZip=true); + static void solve(SimTool* sim_tool, bool bLoadFinal=true, double* paramValues=nullptr); - void init(double* paramValues=0); - void step(double* paramValues=0); + static void init(SimTool* sim_tool, double* paramValues=nullptr); + static void step(SimTool* sim_tool, double* paramValues=nullptr); - double getCurrentTime(); - void setEndTime(double endTime); + static double getCurrentTime(SimTool* sim_tool); + static void setEndTime(SimTool* sim_tool, double endTime); - int getNumVariables(); - string getVariableName(int index); - int getVariableLength(string& varName); - double* getValue(string& varName, int arrayID); // arrayID=0 for "old" and 1 for "current" - void setInitialCondition(string& varName, int dataLength, const double* data); + static int getNumVariables(const SimTool* sim_tool); + static string getVariableName(const SimTool* sim_tool, int index); + static int getVariableLength(const SimTool* sim_tool, string& varName); + static double* getValue(const SimTool* sim_tool, string& varName, int arrayID); // arrayID=0 for "old" and 1 for "current" + static void setInitialCondition(SimTool* sim_tool, string& varName, int dataLength, const double* data); - //void reinit(double *paramValues); - //PdeResultSet* getPdeResultSet(); private: - void loadJMSInfo(istream& ifsInput, int taskID); - void loadModel(istream& ifsInput); - void loadSimulation(istream& ifsInput); - void loadSmoldyn(istream& ifsInput); - VCell::Expression* readExpression(istream& ifsInput, string& var_name, string prefix=""); - VarContext* loadEquation(istream& ifsInput, Structure* structure, Variable* var); - void loadJumpCondition(istream& ifsInput, Membrane*, string& var_name); - void loadPseudoConstants(istream& ifsInput, FastSystemExpression* fastSystem); - void loadFastRates(istream& ifsInput, FastSystemExpression* fastSystem); - void loadFastDependencies(istream& ifsInput, FastSystemExpression* fastSystem); - void loadJacobians(istream& ifsInput, FastSystemExpression* fastSystem); - void loadFastSystem(istream& ifsInput, FastSystemExpression* fastSystem); - void loadFeature(istream& ifsInput, Feature* feature); - void loadMembrane(istream& ifsInput, Membrane*); - void loadSimulationParameters(istream& ifsInput); - void loadMesh(istream& ifsInput); - void loadFieldData(istream& ifsInput); - void loadParameters(istream& ifsInput, int numParams); - void loadSerialScanParameters(istream& ifsInput, int numSerialScanParameters); - void loadSerialScanParameterValues(istream& ifsInput, int numSerialScanParamValues); + static void loadJMSInfo(istream& ifsInput, int taskID); + static VCellModel* loadModel(istream& ifsInput); + SimulationExpression* loadSimulation(SimTool* sim_tool, VCellModel* model, istream& ifsInput); + static void loadSmoldyn(SimTool* sim_tool, istream& ifsInput); + static VCell::Expression* readExpression(istream& ifsInput, string& var_name, const string& prefix=""); + static VarContext* loadEquation(istream& ifsInput, Structure* structure, Variable* var); + static void loadJumpCondition(SimulationExpression* simulation, VCellModel* model, istream& ifsInput, Membrane*, string& var_name); + static void loadPseudoConstants(istream& ifsInput, FastSystemExpression* fastSystem); + static void loadFastRates(istream& ifsInput, FastSystemExpression* fastSystem); + static void loadFastDependencies(istream& ifsInput, FastSystemExpression* fastSystem); + static void loadJacobians(istream& ifsInput, FastSystemExpression* fastSystem); + static void loadFastSystem(istream& ifsInput, FastSystemExpression* fastSystem); + static void loadFeature(SimulationExpression* simulation, istream& ifsInput, Feature* feature); + static void loadMembrane(SimulationExpression* simulation, VCellModel* model, istream& ifsInput, Membrane*); + void loadSimulationParameters(SimTool* sim_tool, istream& ifsInput) const; + void loadMeshFromVcg(VCellModel* model, istream& vcgInput); + static void loadFieldData(SimulationExpression* simulation, istream& ifsInput); + static void loadParameters(SimulationExpression* simulation, istream& ifsInput, int numParams); + static void loadSerialScanParameters(SimulationExpression* simulation, istream& ifsInput, int numSerialScanParameters); + static void loadSerialScanParameterValues(SimTool* simTool, SimulationExpression* simulation, istream& ifsInput, int numSerialScanParamValues); const char* outputPath; - SimTool* simTool; - SimulationExpression *simulation; - VCellModel *model; CartesianMesh *mesh; - int loadSolveRegions(istream& instream, int*& solveRegions); + int loadSolveRegions(VCellModel* model, istream& instream, int*& solveRegions) const; - void loadPostProcessingBlock(istream& ifsInput); + static void loadPostProcessingBlock(SimTool* sim_tool, istream& ifsInput); }; #endif diff --git a/VCell/include/VCELL/FastSystem.h b/VCell/include/VCELL/FastSystem.h index 27b2ebb37..722854687 100644 --- a/VCell/include/VCELL/FastSystem.h +++ b/VCell/include/VCELL/FastSystem.h @@ -8,7 +8,7 @@ #include #include -class Simulation; +class SimulationExpression; class Variable; class FastSystem : public AlgebraicSystem @@ -17,11 +17,11 @@ class FastSystem : public AlgebraicSystem FastSystem(int dimension, int numDependents); void setCurrIndex(long index){currIndex = index;} void updateVars(); - inline int getNumDependents() {return numDependents;} + int getNumDependents() const {return numDependents;} void setDependencyMatrix(int i, int j, double value); virtual void updateDependentVars(); void showVars(); - virtual void resolveReferences(Simulation *sim)=0; + virtual void resolveReferences(SimulationExpression *sim)=0; virtual void setCoordinates(double time_sec, WorldCoord& wc){}; protected: diff --git a/VCell/include/VCELL/FastSystemExpression.h b/VCell/include/VCELL/FastSystemExpression.h index 900a2d8a9..562f12895 100644 --- a/VCell/include/VCELL/FastSystemExpression.h +++ b/VCell/include/VCELL/FastSystemExpression.h @@ -2,11 +2,12 @@ #define FASTSYSTEMEXPRESSION_H #include -#include + using std::string; class SimulationExpression; class SimpleSymbolTable; +class SimTool; namespace VCell { class Expression; } @@ -19,22 +20,22 @@ namespace VCell { class FastSystemExpression : public FastSystem { public: FastSystemExpression(int dimension, int numDepend, SimulationExpression* sim); - ~FastSystemExpression(); - - virtual void resolveReferences(Simulation *sim); - void updateDependentVars(); + ~FastSystemExpression() override; + + void resolveReferences(SimulationExpression *sim) override; + void updateDependentVars() override; void setPseudoConstants(string* symbols, VCell::Expression** expressions); void setFastRateExpressions(VCell::Expression** expressions); void setFastDependencyExpressions(string* symbols, VCell::Expression** expressions); void setJacobianExpressions(VCell::Expression** expressions); - void setCoordinates(double time_sec, WorldCoord& wc); + void setCoordinates(double time_sec, WorldCoord& wc) override; - void initVars(); - void setDependentVariables(string* vars); // must be called before other setters - void setIndependentVariables(string* vars); // must be called before other setters - void updateIndepValues(); - void updateMatrix(); + void initVars(SimTool *sim_tool) override; + void setDependentVariables(string* vars) const; // must be called before other setters + void setIndependentVariables(string* vars) const; // must be called before other setters + void updateIndepValues() const; + void updateMatrix() override; private: string* pseudoSymbols; diff --git a/VCell/include/VCELL/Feature.h b/VCell/include/VCELL/Feature.h index ba6ab15c7..8e807785f 100644 --- a/VCell/include/VCELL/Feature.h +++ b/VCell/include/VCELL/Feature.h @@ -11,11 +11,9 @@ using std::vector; class VolumeVarContextExpression; class VolumeRegionVarContextExpression; -//class VolumeParticleContext; -//class ContourParticleContext; class FastSystem; class Feature; -class Simulation; +class SimulationExpression; class VolumeVariable; class MembraneVariable; class VolumeRegionVariable; @@ -25,16 +23,13 @@ class Feature : public Structure { public: Feature(string& name, unsigned char findex, FeatureHandle handle); - ~Feature(); - void resolveReferences(Simulation *sim); + void resolveReferences(SimulationExpression *sim); virtual void initVolumeValues(long volumeIndex); virtual void initVolumeRegionValues(int volumeRegionIndex); FeatureHandle getHandle(); - unsigned char getIndex() { - return index; - } + unsigned char getIndex() const { return index; } VolumeVarContextExpression *getVolumeVarContext(VolumeVariable *var); VolumeRegionVarContextExpression *getVolumeRegionVarContext(VolumeRegionVariable *var); @@ -42,7 +37,7 @@ class Feature : public Structure void addVolumeVarContext(VolumeVarContextExpression *vc); void addVolumeRegionVarContext(VolumeRegionVarContextExpression *vc); - void reinitConstantValues(); + void reinitConstantValues(SimulationExpression* sim); //VolumeParticleContext *getVolumeParticleContext(){return vpc;} //MembraneParticleContext *getMembraneParticleContext(){return mpc;} diff --git a/VCell/include/VCELL/JumpCondition.h b/VCell/include/VCELL/JumpCondition.h index 95453fb3c..db85927a1 100644 --- a/VCell/include/VCELL/JumpCondition.h +++ b/VCell/include/VCELL/JumpCondition.h @@ -13,19 +13,15 @@ class JumpCondition { public: JumpCondition(Membrane*, VCell::Expression*); - ~JumpCondition(void); + ~JumpCondition(); - VCell::Expression *getExpression() { - return expression; - } - Membrane* getMembrane() { - return membrane; - } + VCell::Expression *getExpression() const { return expression; } + Membrane* getMembrane() const { return membrane; } void bindExpression(SymbolTable*); double evaluateExpression(double* values); double evaluateExpression(SimulationExpression*, MembraneElement*); - void reinitConstantValues(); + void reinitConstantValues(SimulationExpression* sim); private: Membrane* membrane; @@ -33,7 +29,7 @@ class JumpCondition double* constantValue; bool bNeedsXYZ; - bool isConstantExpression(); + bool isConstantExpression(SimulationExpression *sim); }; #endif diff --git a/VCell/include/VCELL/Membrane.h b/VCell/include/VCELL/Membrane.h index 96f3d0481..0aa999efe 100644 --- a/VCell/include/VCELL/Membrane.h +++ b/VCell/include/VCELL/Membrane.h @@ -6,7 +6,7 @@ using std::vector; class Feature; -class Simulation; +class SimulationExpression; struct MembraneElement; class MembraneVarContextExpression; class MembraneRegionVarContextExpression; @@ -17,29 +17,23 @@ class Membrane : public Structure { public: Membrane(string& name, Feature* f1, Feature* f2); - ~Membrane(); - + ~Membrane() override; virtual void initMembraneValues(MembraneElement *membraneElement); virtual void initMembraneRegionValues(int membraneRegionIndex); - MembraneVarContextExpression *getMembraneVarContext(string& membraneVarName); MembraneVarContextExpression *getMembraneVarContext(MembraneVariable *var); void addMembraneVarContext(MembraneVarContextExpression *vc); MembraneRegionVarContextExpression *getMembraneRegionVarContext(MembraneRegionVariable *var); void addMembraneRegionVarContext(MembraneRegionVarContextExpression *vc); - void resolveReferences(Simulation *sim); + void resolveReferences(SimulationExpression *sim); - Feature* getFeature1() { - return feature1; - } - Feature* getFeature2() { - return feature2; - } + Feature* getFeature1() const { return feature1; } + Feature* getFeature2() const { return feature2; } bool inBetween(Feature* f1, Feature* f2); - void reinitConstantValues(); + void reinitConstantValues(SimulationExpression* sim); private: Feature* feature1; diff --git a/VCell/include/VCELL/MembraneEqnBuilderDiffusion.h b/VCell/include/VCELL/MembraneEqnBuilderDiffusion.h index 74c8583a9..c6ecc3928 100644 --- a/VCell/include/VCELL/MembraneEqnBuilderDiffusion.h +++ b/VCell/include/VCELL/MembraneEqnBuilderDiffusion.h @@ -19,21 +19,20 @@ class MembraneEqnBuilderDiffusion : public SparseMatrixEqnBuilder { public: MembraneEqnBuilderDiffusion(MembraneVariable *species, Mesh *mesh); - ~MembraneEqnBuilderDiffusion(); + ~MembraneEqnBuilderDiffusion() override; - void initEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); - void buildEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); - - void postProcess(); + void initEquation(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) override; + void buildEquation(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) override; + void postProcess() override; private: vector > periodicPairs; // map of minus and plus pairs of periodic boundary points. - void initEquation_Periodic(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); - void buildEquation_Periodic(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); + void initEquation_Periodic(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); + void buildEquation_Periodic(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); bool bPreProcessed; - void preProcess(); + void preProcess(VCellModel* model); private: double computeDiffusionConstant(int meIndex, int neighborIndex); diff --git a/VCell/include/VCELL/MembraneEqnBuilderForward.h b/VCell/include/VCELL/MembraneEqnBuilderForward.h index 85f6bef19..cc5d648cf 100644 --- a/VCell/include/VCELL/MembraneEqnBuilderForward.h +++ b/VCell/include/VCELL/MembraneEqnBuilderForward.h @@ -19,9 +19,9 @@ class MembraneEqnBuilderForward : public EqnBuilder Mesh *mesh, ODESolver *solver); - void initEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) + void initEquation(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) override {} - void buildEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); + void buildEquation(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) override; private: ODESolver* odeSolver; diff --git a/VCell/include/VCELL/MembraneRegionEqnBuilder.h b/VCell/include/VCELL/MembraneRegionEqnBuilder.h index 69db4d245..eda56af04 100644 --- a/VCell/include/VCELL/MembraneRegionEqnBuilder.h +++ b/VCell/include/VCELL/MembraneRegionEqnBuilder.h @@ -16,10 +16,9 @@ class MembraneRegionEqnBuilder : public EqnBuilder { public: MembraneRegionEqnBuilder(MembraneRegionVariable *var, CartesianMesh *mesh, ODESolver *solver); - void initEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) - { } - - void buildEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); + void initEquation(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) override + { } + void buildEquation(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) override; private: ODESolver* odeSolver; diff --git a/VCell/include/VCELL/MembraneRegionVarContextExpression.h b/VCell/include/VCELL/MembraneRegionVarContextExpression.h index 5ff115415..782d6aae9 100644 --- a/VCell/include/VCELL/MembraneRegionVarContextExpression.h +++ b/VCell/include/VCELL/MembraneRegionVarContextExpression.h @@ -15,15 +15,15 @@ struct MembraneElement; class MembraneRegionVarContextExpression : public VarContext { public: - double getInitialValue(long index); + double getInitialValue(long index) override; double getMembraneReactionRate(MembraneElement *element); double getUniformRate(MembraneRegion *region); MembraneRegionVarContextExpression(Membrane *membrane, MembraneRegionVariable* var); - void resolveReferences(Simulation *sim); + void resolveReferences(SimulationExpression *sim) override; protected: - bool isNullExpressionOK(int expIndex); + bool isNullExpressionOK(int expIndex) const override; }; #endif diff --git a/VCell/include/VCELL/MembraneVarContextExpression.h b/VCell/include/VCELL/MembraneVarContextExpression.h index be398441d..3dd394f42 100644 --- a/VCell/include/VCELL/MembraneVarContextExpression.h +++ b/VCell/include/VCELL/MembraneVarContextExpression.h @@ -15,7 +15,7 @@ class MembraneVarContextExpression : public VarContext public: MembraneVarContextExpression(Membrane *membrane, MembraneVariable* var); - void resolveReferences(Simulation *sim); + void resolveReferences(SimulationExpression *sim) override; double getInitialValue(MembraneElement *element); double getMembraneReactionRate(MembraneElement *element); @@ -36,7 +36,7 @@ class MembraneVarContextExpression : public VarContext double getZpBoundaryFlux(MembraneElement *element); protected: - bool isNullExpressionOK(int expIndex); + bool isNullExpressionOK(int expIndex) const override; }; #endif diff --git a/VCell/include/VCELL/Mesh.h b/VCell/include/VCELL/Mesh.h index e392267e8..b2e7cc52d 100644 --- a/VCell/include/VCELL/Mesh.h +++ b/VCell/include/VCELL/Mesh.h @@ -29,10 +29,6 @@ class Mesh { virtual WorldCoord getVolumeWorldCoord(long volumeIndex) = 0; virtual WorldCoord getMembraneWorldCoord(long membraneIndex) = 0; virtual WorldCoord getMembraneWorldCoord(MembraneElement *element) = 0; - //virtual long getVolumeIndex(WorldCoord coord)=0; - - //virtual double getInsideOld(VolumeVariable *var, MembraneElement *element)=0; - //virtual double getOutsideOld(VolumeVariable *var, MembraneElement *element)=0; virtual void showSummary(FILE *fp) { fprintf(fp, "Mesh::showSummary()...\n"); @@ -45,19 +41,7 @@ class Mesh { MembraneElement *getMembraneElements(); long getNumMembraneElements(); - int getDimension() { - return dimension; - } - - //long getNumContourElements(){return (int)contourElements.size();} - //ContourElement *getContourElement(long index); - - //void addElementToVolumeList(long volumeIndex, ContourElement *element); - //vector getContourElementList(long index){return volumeLists[index];} - //virtual void setVolumeLists()=0; - - //int getNumContourSubdomains(){return (int)pContourSubdomains.size();} - //ContourSubdomain *getContourSubdomain(int i); //returns the ith + int getDimension() const { return dimension; } SparseMatrixPCG* getMembraneCoupling() { if (membraneElementCoupling) { @@ -70,7 +54,6 @@ class Mesh { return membraneElementCoupling; } - virtual int getMembraneNeighborMask(long meindex) = 0; virtual int getMembraneNeighborMask(MembraneElement* element) = 0; virtual double* getMembraneFluxArea(long index) = 0; @@ -80,17 +63,11 @@ class Mesh { virtual ~Mesh(); - //virtual bool sampleContours(); // this should go away just like sampleGeometry() - //void addContourSubdomain(ContourSubdomain *cs); - VolumeElement *pVolumeElement; MembraneElement *pMembraneElement; - //vector contourElements; long numMembrane; long numVolume; int dimension; - //vector pContourSubdomains; - //vector *volumeLists; double captureNeighborhood; SparseMatrixPCG* membraneElementCoupling; virtual void computeMembraneCoupling(void) = 0; diff --git a/VCell/include/VCELL/ODESolver.h b/VCell/include/VCELL/ODESolver.h index 6ad58c7d4..900cca374 100644 --- a/VCell/include/VCELL/ODESolver.h +++ b/VCell/include/VCELL/ODESolver.h @@ -16,11 +16,11 @@ class ODESolver : public Solver int *solveRegions); virtual ~ODESolver(); - virtual void solveEqn(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime); + void solveEqn(SimTool* sim_tool, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime) override; - double *getRates() { return rate; } - long getArraySize() {return arraySize;} - long getGlobalIndex(long arrayIndex) {return Gridmap[arrayIndex];} + double *getRates() const { return rate; } + long getArraySize() const {return arraySize;} + long getGlobalIndex(long arrayIndex) const {return Gridmap[arrayIndex];} protected: Mesh *mesh; diff --git a/VCell/include/VCELL/PDESolver.h b/VCell/include/VCELL/PDESolver.h index a84e3b67f..cdc2c3449 100644 --- a/VCell/include/VCELL/PDESolver.h +++ b/VCell/include/VCELL/PDESolver.h @@ -13,20 +13,17 @@ class MembraneVariable; class EqnBuilder; class Mesh; class CartesianMesh; +class VCellModel; class PDESolver : public Solver { public: PDESolver(Variable *var, bool bTimeDependent); - virtual ~PDESolver(); - virtual void initEqn(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime); - virtual bool isPDESolver() { - return true; - } - bool isTimeDependent() { - return bTimeDependent; - } + void initEqn(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime) override; + + bool isPDESolver() override { return true; } + bool isTimeDependent() const { return bTimeDependent; } private: bool bTimeDependent; }; diff --git a/VCell/include/VCELL/PostProcessingHdf5Writer.h b/VCell/include/VCELL/PostProcessingHdf5Writer.h index 5ef0418c1..aa8df079b 100644 --- a/VCell/include/VCELL/PostProcessingHdf5Writer.h +++ b/VCell/include/VCELL/PostProcessingHdf5Writer.h @@ -23,7 +23,7 @@ class PostProcessingHdf5Writer PostProcessingHdf5Writer(std::string fileName, PostProcessingBlock* postProcessingBlock); virtual ~PostProcessingHdf5Writer(); - void writeOutput(); + void writeOutput(SimTool* sim_tool); static const char* PPGroupName; static const char* TimesDataSetName; diff --git a/VCell/include/VCELL/Scheduler.h b/VCell/include/VCELL/Scheduler.h index ae4820455..7ce752936 100644 --- a/VCell/include/VCELL/Scheduler.h +++ b/VCell/include/VCELL/Scheduler.h @@ -5,16 +5,20 @@ #ifndef SCHEDULER_H #define SCHEDULER_H +class SimTool; +class VCellModel; class Simulation; class Scheduler { public: Scheduler(Simulation *Asim); - virtual void iterate()=0; - virtual void initValues(); - void solveFastSystem(int startVolIndex, int VolSize, int startMemIndex, int MemSize); - bool hasFastSystem() { return bHasFastSystem; } + virtual ~Scheduler() = default; + + virtual void iterate(SimTool* sim_tool)=0; + virtual void initValues(VCellModel* model); + void solveFastSystem(SimTool* sim_tool, int startVolIndex, int VolSize, int startMemIndex, int MemSize); + bool hasFastSystem(VCellModel* model) const { return bHasFastSystem; } protected: Simulation *sim; diff --git a/VCell/include/VCELL/SerialScheduler.h b/VCell/include/VCELL/SerialScheduler.h index 4ef83c09a..0f41a39b5 100644 --- a/VCell/include/VCELL/SerialScheduler.h +++ b/VCell/include/VCELL/SerialScheduler.h @@ -13,7 +13,7 @@ class SerialScheduler : public Scheduler { public: SerialScheduler(Simulation *Asim); - virtual void iterate(); + void iterate(SimTool* sim_tool) override; }; #endif diff --git a/VCell/include/VCELL/SimTool.h b/VCell/include/VCELL/SimTool.h index 0c5cca5e5..8ca33411e 100644 --- a/VCell/include/VCELL/SimTool.h +++ b/VCell/include/VCELL/SimTool.h @@ -7,12 +7,14 @@ /////////////////////////////////////////////////////////// #ifndef SIMTOOL_H #define SIMTOOL_H +#include #include #include #include + #ifndef DIRECTORY_SEPARATOR #if ( defined(WIN32) || defined(WIN64) ) #define DIRECTORY_SEPARATOR '\\' @@ -23,70 +25,61 @@ class VCellModel; class Mesh; -class Simulation; +class SimulationExpression; class Variable; class PostProcessingHdf5Writer; class SimTool { public: - static SimTool* getInstance(); - static void create(); - ~SimTool(); + SimTool(); + virtual ~SimTool(); virtual void start(); virtual void loadFinal(); void setModel(VCellModel* model); - void setSimulation(Simulation* sim); + void setSimulation(SimulationExpression* sim); void setTimeStep(double period); void setSmoldynStepMultiplier(int steps); void setCheckSpatiallyUniform(); - bool isCheckingSpatiallyUniform() { return bCheckSpatiallyUniform; } + bool isCheckingSpatiallyUniform() const { return bCheckSpatiallyUniform; } void setEndTimeSec(double timeSec) { simEndTime = timeSec; } - double getEndTime() { return simEndTime; } + double getEndTime() const { return simEndTime; } void setKeepEvery(int ke) { keepEvery = ke; } void setKeepAtMost(int kam) { keepAtMost = kam; } - void setBaseFilename(char *fname); - char* getBaseFileName() { - return baseFileName; - } - char* getBaseDirName() { - return baseDirName; - } - void setStoreEnable(bool enable) { - bStoreEnable = enable; - } - void setFileCompress(bool compress) { - bSimFileCompress = compress; - } + void setBaseFilename(const std::filesystem::path& fname); + [[nodiscard]] std::filesystem::path getBaseFileName() const { return baseFileName; } + [[nodiscard]] std::filesystem::path getBaseDirName() const { return baseDirName; } + void setStoreEnable(bool enable) { bStoreEnable = enable; } + void setFileCompress(bool compress) { bSimFileCompress = compress; } void requestNoZip(); - Simulation* getSimulation() { return simulation; } - VCellModel* getModel() { return vcellModel; } - bool checkStopRequested(); + SimulationExpression* getSimulation() const { return simulation; } + VCellModel* getModel() const { return vcellModel; } + static bool checkStopRequested(); TimerHandle getTimerHandle(string& timerName); - void startTimer(TimerHandle hnd); - void stopTimer(TimerHandle hnd); - double getElapsedTimeSec(TimerHandle hnd); + void startTimer(TimerHandle hnd) const; + void stopTimer(TimerHandle hnd) const; + double getElapsedTimeSec(TimerHandle hnd) const; virtual void showSummary(FILE *fp); - void setSolver(string& s); - bool isSundialsPdeSolver(); + void setSolver(const string& s); + bool isSundialsPdeSolver() const; void setDiscontinuityTimes(int num, double* times) { numDiscontinuityTimes = num; discontinuityTimes = times; } - int getNumDiscontinuityTimes() { return numDiscontinuityTimes; } - double* getDiscontinuityTimes() { return discontinuityTimes; } + int getNumDiscontinuityTimes() const { return numDiscontinuityTimes; } + double* getDiscontinuityTimes() const { return discontinuityTimes; } - void setSundialsSolverOptions(SundialsSolverOptions& sso) { + void setSundialsSolverOptions(const SundialsSolverOptions& sso) { sundialsSolverOptions = sso; } - const SundialsSolverOptions& getSundialsSolverOptions() { return sundialsSolverOptions; } + const SundialsSolverOptions& getSundialsSolverOptions() const { return sundialsSolverOptions; } void setSpatiallyUniformErrorTolerance(double atol, double rtol) { spatiallyUniformAbsTol = atol; @@ -96,32 +89,30 @@ class SimTool { void setPCGRelativeErrorTolerance(double rtol) { pcgRelTol = rtol; } - double getPCGRelativeErrorTolerance() { + double getPCGRelativeErrorTolerance() const + { return pcgRelTol; } - double getSimStartTime() { return simStartTime; } + double getSimStartTime() const { return simStartTime; } void setSundialsOneStepOutput() { bSundialsOneStepOutput = true; } - bool isSundialsOneStepOutput() { return bSundialsOneStepOutput; } + bool isSundialsOneStepOutput() const { return bSundialsOneStepOutput; } void setSerialParameterScans(int numScans, double** values); void setLoadFinal(bool b) { bLoadFinal = b; } - void checkTaskIdLockFile(); + void checkTaskIdLockFile() const; - void setSmoldynInputFile(string& inputfile); + void setSmoldynInputFile(const string& inputfile); private: - SimTool(); - - FILE* lockForReadWrite(); + FILE* lockForReadWrite() const; bool checkSpatiallyUniform(Variable*); void updateLog(double progress,double time,int iteration); void clearLog(); - int getZipCount(char* zipFileName); - int getZipCount(const std::string* zipFileName); + static int getZipCount(const std::filesystem::path& zipFileName); void start1(); void copyParticleCountsToConcentration(); @@ -129,7 +120,7 @@ class SimTool { bool bSimZip; VCellModel* vcellModel; - Simulation *simulation; + SimulationExpression *simulation; Timer *_timer; bool bSimFileCompress; @@ -140,10 +131,10 @@ class SimTool { int smoldynStepMultiplier; int keepEvery; bool bStoreEnable; - char* baseFileName; + std::filesystem::path baseFileName; int simFileCount; - char* baseSimName; - char* baseDirName; + std::filesystem::path baseSimName; + std::filesystem::path baseDirName; int zipFileCount; string solver; double* discontinuityTimes; diff --git a/VCell/include/VCELL/Simulation.h b/VCell/include/VCELL/Simulation.h index 22dccbb85..20166d1e4 100644 --- a/VCell/include/VCELL/Simulation.h +++ b/VCell/include/VCELL/Simulation.h @@ -5,20 +5,18 @@ #ifndef SIMULATION_H #define SIMULATION_H -#include #include +#include using std::vector; using std::string; -//#define PARTICLE_ALL -1 +class Simulation; +class VCellModel; -//class VolumeParticleContext; -//class MembraneParticleContext; -//class ContourParticleContext; -//class Particle; class VolumeVariable; class Variable; -class Mesh; +class CartesianMesh; +class SimTool; class Solver; class Scheduler; class PostProcessingBlock; @@ -26,22 +24,24 @@ class PostProcessingBlock; class Simulation { public: - Simulation(Mesh *mesh); - ~Simulation(); + explicit Simulation(CartesianMesh *mesh); + virtual ~Simulation(); - virtual void resolveReferences(); - void initSimulation(); // initializes to t=0 - void iterate(); // computes 1 time step + virtual void resolveReferences(SimTool* sim_tool); + virtual void initSimulation(SimTool* sim_tool) = 0; // initializes to t=0 + void iterate(SimTool* sim_tool); // computes 1 time step virtual void update(); // copies new to old values - double getTime_sec(); + double getTime_sec(SimTool* sim_tool); void setCurrIteration(int curriter) { currIteration = curriter; } - int getCurrIteration() { + int getCurrIteration() const + { return currIteration; } - double getDT_sec() { + double getDT_sec() const + { return _dT_sec; } void setDT_sec(double dT) { @@ -58,28 +58,28 @@ class Simulation Variable *getVariableFromName(string& name); Variable *getVariableFromName(char* name); Solver *getSolverFromVariable(Variable *var); - Mesh *getMesh() { + CartesianMesh *getMesh() const + { return _mesh; } Solver* getSolver(int index); - //void addParticle(Particle *particle); - //long getNumParticles() { - // return (int)globalParticleList.size(); - //} - int getNumVariables() { - return (int)varList.size(); + int getNumVariables() const + { + return static_cast(varList.size()); } - int getNumSolvers() { - return (int)solverList.size(); + int getNumSolvers() const + { + return static_cast(solverList.size()); } void addVariable(Variable *var); void addSolver(Solver *solver); - void setSimStartTime(double st); + void setSimStartTime(SimTool* sim_tool, double st); virtual void createPostProcessingBlock()=0; - PostProcessingBlock* getPostProcessingBlock() { + PostProcessingBlock* getPostProcessingBlock() const + { return postProcessingBlock; } @@ -91,7 +91,7 @@ class Simulation vector solverList; vector varList; //vector globalParticleList; - Mesh *_mesh; + CartesianMesh *_mesh; bool _advanced; bool _initEquations; diff --git a/VCell/include/VCELL/SimulationExpression.h b/VCell/include/VCELL/SimulationExpression.h index 04283a36b..d8b6014c9 100644 --- a/VCell/include/VCELL/SimulationExpression.h +++ b/VCell/include/VCELL/SimulationExpression.h @@ -27,20 +27,21 @@ class MembraneParticleVariable; class SimulationExpression : public Simulation { public: - SimulationExpression(Mesh *mesh); - ~SimulationExpression(); + explicit SimulationExpression(CartesianMesh *mesh); + ~SimulationExpression() override; - void resolveReferences(); // create symbol table - void update(); // copies new to old values - void advanceTimeOn(); - void advanceTimeOff(); + void initSimulation(SimTool* sim_tool) override; + void resolveReferences(SimTool* sim_tool) override; // create symbol table + void update(SimTool* sim_tool); // copies new to old values + void advanceTimeOn(SimTool* sim_tool); + void advanceTimeOff(SimTool* sim_tool); - void writeData(const char *filename, bool bCompress); + void writeData(const char *filename, bool bCompress) override; void addFieldData(FieldData* fd) { fieldDataList.push_back(fd); } - int getNumFields() { return (int)fieldDataList.size(); } + int getNumFields() const { return static_cast(fieldDataList.size()); } string* getFieldSymbols(); void populateFieldValues(double* darray, int index); @@ -48,36 +49,41 @@ class SimulationExpression : public Simulation randomVarList.push_back(rv); } RandomVariable* getRandomVariableFromName(char* rvname); - int getNumRandomVariables() { - return (int)randomVarList.size(); + int getNumRandomVariables() const + { + return static_cast(randomVarList.size()); } - RandomVariable* getRandomVariable(int index) { + RandomVariable* getRandomVariable(int index) const + { return randomVarList[index]; } void populateRandomValues(double* darray, int index); - int* getIndices() { return indices; }; - SymbolTable* getSymbolTable() { + int* getIndices() const { return indices; }; + SymbolTable* getSymbolTable() const + { return symbolTable; }; void setCurrentCoordinate(WorldCoord& wc); - bool isVolumeVariableDefinedInRegion(int volVarIndex, int regionIndex) { - if (volVariableRegionMap[volVarIndex] == 0) { + bool isVolumeVariableDefinedInRegion(int volVarIndex, int regionIndex) const + { + if (volVariableRegionMap[volVarIndex] == nullptr) { return true; } return volVariableRegionMap[volVarIndex][regionIndex]; } void addParameter(string& param); - void setParameterValues(double* paramValues); - int getNumParameters() { - return (int)paramList.size(); + void setParameterValues(SimTool* sim_tool, double* paramValues); + int getNumParameters() const + { + return static_cast(paramList.size()); } void populateParameterValues(double* darray); // right now bSolveRegion is only applicable for volume variables - void addVariable(Variable *var, bool* bSolveRegions=0); + void addVariable(Variable *var, bool* bSolveRegions=nullptr); void addVolumeVariable(VolumeVariable *var, bool* bSolveRegions); void addVolumeParticleVariable(VolumeParticleVariable *var); void addMembraneVariable(MembraneVariable *var); @@ -85,57 +91,61 @@ class SimulationExpression : public Simulation void addVolumeRegionVariable(VolumeRegionVariable *var); void addMembraneRegionVariable(MembraneRegionVariable *var); - double* getDiscontinuityTimes() { return discontinuityTimes; } + double* getDiscontinuityTimes() const { return discontinuityTimes; } void setDiscontinuityTimes(double* stopTimes) { discontinuityTimes=stopTimes; } - int getNumVolVariables() { return volVarSize; } - VolumeVariable* getVolVariable(int i) { return volVarList[i]; } + int getNumVolVariables() const { return volVarSize; } + VolumeVariable* getVolVariable(int i) const { return volVarList[i]; } - int getNumMemVariables() { return memVarSize; } - MembraneVariable* getMemVariable(int i) { return memVarList[i]; } + int getNumMemVariables() const { return memVarSize; } + MembraneVariable* getMemVariable(int i) const { return memVarList[i]; } - int getNumVolRegionVariables() { return volRegionVarSize; } - VolumeRegionVariable* getVolRegionVariable(int i) { return volRegionVarList[i]; } + int getNumVolRegionVariables() const { return volRegionVarSize; } + VolumeRegionVariable* getVolRegionVariable(int i) const { return volRegionVarList[i]; } - int getNumMemRegionVariables() { return memRegionVarSize; } - MembraneRegionVariable* getMemRegionVariable(int i) { return memRegionVarList[i]; } + int getNumMemRegionVariables() const { return memRegionVarSize; } + MembraneRegionVariable* getMemRegionVariable(int i) const { return memRegionVarList[i]; } - int getNumVolParticleVariables() { + int getNumVolParticleVariables() const + { return volParticleVarSize; } - int getNumMemParticleVariables() { + int getNumMemParticleVariables() const + { return memParticleVarSize; } - int getNumMemPde() { return numMemPde; } - int getNumVolPde() { return numVolPde; } + int getNumMemPde() const { return numMemPde; } + int getNumVolPde() const { return numVolPde; } void setHasTimeDependentDiffusionAdvection() { bHasTimeDependentDiffusionAdvection = true; } - bool hasTimeDependentDiffusionAdvection() { return bHasTimeDependentDiffusionAdvection; } + bool hasTimeDependentDiffusionAdvection() const { return bHasTimeDependentDiffusionAdvection; } void setPSFFieldDataIndex(int idx) { psfFieldDataIndex = idx; } - FieldData* getPSFFieldData() { + FieldData* getPSFFieldData() const + { if (psfFieldDataIndex >= 0) { return fieldDataList[psfFieldDataIndex]; } - - return 0; + return nullptr; } bool isParameter(string& symbol); // can be serial scan parameter or opt parameter bool isVariable(string& symbol); - int getNumRegionSizeVariables() { + int getNumRegionSizeVariables() const + { return numRegionSizeVars; } - RegionSizeVariable* getRegionSizeVariable(int index) { + RegionSizeVariable* getRegionSizeVariable(int index) const + { return regionSizeVarList[index]; } void populateRegionSizeVariableValues(double* darray, bool bVolumeRegion, int regionIndex); - void createPostProcessingBlock(); + void createPostProcessingBlock() override; void populateFieldValuesNew(double* darray, int index); void populateRegionSizeVariableValuesNew(double* darray, bool bVolumeRegion, int regionIndex); @@ -143,25 +153,25 @@ class SimulationExpression : public Simulation void populateRandomValuesNew(double* darray, int index); void populateParticleVariableValuesNew(double* array, bool bVolume, int index); - int symbolIndexOfT() { return symbolIndexOffset_T; } - int symbolIndexOfXyz() { return symbolIndexOffset_Xyz; } - int symbolIndexOfVolVar() { return symbolIndexOffset_VolVar; } - int symbolIndexOfMemVar() { return symbolIndexOffset_MemVar; } - int symbolIndexOfVolRegionVar() { return symbolIndexOffset_VolRegionVar; } - int symbolIndexOfMemRegionVar() { return symbolIndexOffset_MemRegionVar; } - int symbolIndexOfVolParticleVar() { return symbolIndexOffset_VolParticleVar; } - int symbolIndexOfMemParticleVar() { return symbolIndexOffset_MemParticleVar; } - int symbolIndexOfRegionSizeVariable() { return symbolIndexOffset_RegionSizeVariable; } - int symbolIndexOfFieldData() { return symbolIndexOffset_FieldData; } - int symbolIndexOfRandomVar() { return symbolIndexOffset_RandomVar; } - int symbolIndexOfParameters() { return symbolIndexOffset_Parameters; } - int numOfSymbols() { return numSymbols;} + int symbolIndexOfT() const { return symbolIndexOffset_T; } + int symbolIndexOfXyz() const { return symbolIndexOffset_Xyz; } + int symbolIndexOfVolVar() const { return symbolIndexOffset_VolVar; } + int symbolIndexOfMemVar() const { return symbolIndexOffset_MemVar; } + int symbolIndexOfVolRegionVar() const { return symbolIndexOffset_VolRegionVar; } + int symbolIndexOfMemRegionVar() const { return symbolIndexOffset_MemRegionVar; } + int symbolIndexOfVolParticleVar() const { return symbolIndexOffset_VolParticleVar; } + int symbolIndexOfMemParticleVar() const { return symbolIndexOffset_MemParticleVar; } + int symbolIndexOfRegionSizeVariable() const { return symbolIndexOffset_RegionSizeVariable; } + int symbolIndexOfFieldData() const { return symbolIndexOffset_FieldData; } + int symbolIndexOfRandomVar() const { return symbolIndexOffset_RandomVar; } + int symbolIndexOfParameters() const { return symbolIndexOffset_Parameters; } + int numOfSymbols() const { return numSymbols;} private: SymbolTable* symbolTable; int* indices; - void createSymbolTable(); + void createSymbolTable(SimTool* sim_tool); ScalarValueProxy* valueProxyTime; ScalarValueProxy* valueProxyX; @@ -200,7 +210,7 @@ class SimulationExpression : public Simulation int psfFieldDataIndex; - void reinitConstantValues(); + void reinitConstantValues(SimTool* sim_tool); int numSymbols; int symbolIndexOffset_T; diff --git a/VCell/include/VCELL/Solver.h b/VCell/include/VCELL/Solver.h index 0c00fee45..b40611832 100644 --- a/VCell/include/VCELL/Solver.h +++ b/VCell/include/VCELL/Solver.h @@ -11,19 +11,23 @@ class MembraneVariable; class EqnBuilder; class Mesh; class CartesianMesh; +class SimTool; +class Simulation; +class VCellModel; class Solver { public: - Solver(Variable *variable); + explicit Solver(Variable *variable); + virtual ~Solver() = default; - virtual void initEqn(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime); - virtual void buildEqn(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime); - virtual void solveEqn(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime)=0; + virtual void initEqn(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime); + virtual void buildEqn(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime); + virtual void solveEqn(SimTool* sim_tool, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime)=0; - Variable *getVar() { return var; } + Variable *getVar() const { return var; } void setEqnBuilder(EqnBuilder *builder) { eqnBuilder = builder; } - EqnBuilder *getEqnBuilder() { return eqnBuilder; } + EqnBuilder *getEqnBuilder() const { return eqnBuilder; } virtual bool isPDESolver() { return false; } protected: diff --git a/src/SolverMain.h b/VCell/include/VCELL/SolverMain.h similarity index 56% rename from src/SolverMain.h rename to VCell/include/VCELL/SolverMain.h index 5b69637f4..fb6cb2565 100644 --- a/src/SolverMain.h +++ b/VCell/include/VCELL/SolverMain.h @@ -8,6 +8,6 @@ #include std::string version(); -int solve(const std::string& inputFilename, const std::string& outputDir); +int solve(const std::string& fvInputFilename, const std::string& vcgInputFilename, const std::string& outputDir); #endif //SOLVERMAIN_H diff --git a/VCell/include/VCELL/SparseLinearSolver.h b/VCell/include/VCELL/SparseLinearSolver.h index 7b8e24299..9fc87fd41 100644 --- a/VCell/include/VCELL/SparseLinearSolver.h +++ b/VCell/include/VCELL/SparseLinearSolver.h @@ -10,9 +10,9 @@ class SparseLinearSolver : public PDESolver { public: SparseLinearSolver(Variable *Var, SparseMatrixEqnBuilder* eqnbuilder, double rtol, bool AbTimeDependent); - ~SparseLinearSolver(); + ~SparseLinearSolver() override; - virtual void solveEqn(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime); + void solveEqn(SimTool* sim_tool, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime) override; private: bool enableRetry; @@ -20,7 +20,7 @@ class SparseLinearSolver : public PDESolver double pcgRelErr; protected: - int* PCGSolve(bool bRecomputeIncompleteFactorization); + int* PCGSolve(SimTool* sim_tool, bool bRecomputeIncompleteFactorization); SparseMatrixEqnBuilder* smEqnBuilder; double* pcg_workspace; long nWork; diff --git a/VCell/include/VCELL/SparseVolumeEqnBuilder.h b/VCell/include/VCELL/SparseVolumeEqnBuilder.h index 01edab545..ceccbf9ff 100644 --- a/VCell/include/VCELL/SparseVolumeEqnBuilder.h +++ b/VCell/include/VCELL/SparseVolumeEqnBuilder.h @@ -23,12 +23,12 @@ struct CoupledNeighbors; class SparseVolumeEqnBuilder : public SparseMatrixEqnBuilder { public: - SparseVolumeEqnBuilder(VolumeVariable *species, CartesianMesh *mesh, bool bNoConvection, int numSolveRegions=0, int* solveRegions=0); - ~SparseVolumeEqnBuilder(); + SparseVolumeEqnBuilder(VolumeVariable *species, CartesianMesh *mesh, bool bNoConvection, int numSolveRegions=0, int* solveRegions=nullptr); + ~SparseVolumeEqnBuilder() override; - void initEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); - void buildEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); - void postProcess(); + void initEquation(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) override; + void buildEquation(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) override; + void postProcess() override; private: bool bSymmetricStorage; // no convection, A would be symmetric @@ -52,8 +52,8 @@ class SparseVolumeEqnBuilder : public SparseMatrixEqnBuilder void init(); void computeLHS(int index, double* lambdas, double& Aii, int& numCols, int* columnIndices, double* columnValues, bool& bSort); - double computeRHS(int index, double deltaTime, double* lambdas, double bInit); - void preProcess(); + double computeRHS(Simulation* sim, int index, double deltaTime, double* lambdas, double bInit); + void preProcess(VCellModel* model); bool checkPeriodicCoupledPairsInRegions(int indexm, int indexp); }; diff --git a/VCell/include/VCELL/Structure.h b/VCell/include/VCELL/Structure.h index 002ad8484..55fad1ddb 100644 --- a/VCell/include/VCELL/Structure.h +++ b/VCell/include/VCELL/Structure.h @@ -13,8 +13,8 @@ class FastSystem; class Structure { public: - Structure(string& Aname); - ~Structure(void); + explicit Structure(string& Aname); + virtual ~Structure() = default; virtual BoundaryType getXmBoundaryType() { return boundaryType[0]; } virtual BoundaryType getXpBoundaryType() { return boundaryType[1]; } @@ -31,17 +31,19 @@ class Structure virtual void setZpBoundaryType(BoundaryType bt) { boundaryType[5] = bt; } const string& getName() { return name; } - int getNumRegions() { - return (int)regionList.size(); + int getNumRegions() const + { + return static_cast(regionList.size()); } void addRegion(Region* r); - Region* getRegion(int i) { + Region* getRegion(int i) const + { return regionList.at(i); } int getNumElements(); void setFastSystem(FastSystem* arg_fastSystem) {fastSystem = arg_fastSystem; } - FastSystem *getFastSystem(){ return fastSystem; } + FastSystem *getFastSystem() const { return fastSystem; } protected: string name; diff --git a/VCell/include/VCELL/SundialsPdeScheduler.h b/VCell/include/VCELL/SundialsPdeScheduler.h index 5ae94f338..8402289b8 100644 --- a/VCell/include/VCELL/SundialsPdeScheduler.h +++ b/VCell/include/VCELL/SundialsPdeScheduler.h @@ -22,11 +22,11 @@ class Feature; class SundialsPdeScheduler : public Scheduler { public: - SundialsPdeScheduler(Simulation *sim, const SundialsSolverOptions& sso, int numDisTimes, double* disTimes, bool bDefaultOutput); - ~SundialsPdeScheduler(); + SundialsPdeScheduler(SimulationExpression *sim, const SundialsSolverOptions& sso, int numDisTimes, double* disTimes, bool bDefaultOutput); + ~SundialsPdeScheduler() override; - void iterate(); - double getCurrentTime() { return currentTime; } + void iterate(SimTool* sim_tool) override; + double getCurrentTime() const { return currentTime; } void setSimStartTime(double st) { currentTime = st; } @@ -39,15 +39,15 @@ class SundialsPdeScheduler : public Scheduler static int RHS_callback(realtype t, N_Vector y, N_Vector r, void *fdata); int CVodeRHS(double t, double* yinput, double* rhs); - static int pcSetup_callback(realtype t, N_Vector y, N_Vector fy, booleantype jok, booleantype *jcurPtr, realtype gamma, - void *P_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); - int pcSetup(realtype t, N_Vector y, N_Vector fy, booleantype jok, booleantype *jcurPtr, realtype gamma); - /* + static int pcSetup_callback(realtype t, N_Vector y, N_Vector fy, booleantype jok, booleantype *jcurPtr, realtype gamma, + void *P_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); + int pcSetup(realtype t, N_Vector y, N_Vector fy, booleantype jok, booleantype *jcurPtr, realtype gamma); + /* Definition typedef int (*CVSpilsPrecSetupFn)(realtype t, N Vector y, N Vector fy, booleantype jok, booleantype *jcurPtr, realtype gamma, void *p_data, N_Vector tmp1, N_Vector tmp2, - N_Vector tmp3); + N_Vector tmp3); Purpose This function evaluates and/or preprocesses Jacobian-related data needed by the pre- conditioner. @@ -94,9 +94,9 @@ class SundialsPdeScheduler : public Scheduler If the user's CVSpilsPrecSetupFn function uses difference quotient approximations, it may need to access quantities not in the call list. These include the current stepsize, the error weights, etc. To obtain these, use the CVodeGet* functions described in $5.5.7.1. - The unit roundoff can be accessed as UNIT ROUNDOFF defined in sundials_types.h. - */ - static int pcSolve_callback(realtype t, N_Vector y, N_Vector fy, N_Vector r, N_Vector z, realtype gamma, realtype delta, + The unit roundoff can be accessed as UNIT ROUNDOFF defined in sundials_types.h. + */ + static int pcSolve_callback(realtype t, N_Vector y, N_Vector fy, N_Vector r, N_Vector z, realtype gamma, realtype delta, int lr, void *P_data, N_Vector tmp); int pcSolve(realtype t, N_Vector y, N_Vector fy, N_Vector r, N_Vector z, realtype gamma, realtype delta, int lr); /* @@ -112,19 +112,19 @@ class SundialsPdeScheduler : public Scheduler fy is the current value of the vector f(t; y). r is the right-hand side vector of the linear system. z is the output vector computed. - gamma is the scalar ° appearing in the Newton matrix M = I - gammaJ. + gamma is the scalar � appearing in the Newton matrix M = I - gammaJ. delta is an input tolerance to be used if an iterative method is employed in the solu- tion. In that case, the residual vector Res = r - Pz of the system should be - made less than delta in weighted l2 norm, i.e., pPi(Resi ¢ ewti)2 < delta. + made less than delta in weighted l2 norm, i.e., pPi(Resi � ewti)2 < delta. To obtain the N Vector ewt, call CVodeGetErrWeights (see x5.5.7.1). - lr is an input °ag indicating whether the preconditioner solve function is to use + lr is an input �ag indicating whether the preconditioner solve function is to use the left preconditioner (lr = 1) or the right preconditioner (lr = 2); p_data is a pointer to user data | the same as the p data parameter passed to the function CVSp*SetPreconditioner. tmp is a pointer to memory allocated for a variable of type N Vector which can be used for work space. Return value - The value to be returned by the preconditioner solve function is a °ag indicating whether + The value to be returned by the preconditioner solve function is a �ag indicating whether it was successful. This value should be 0 if successful, positive for a recoverable error (in which case the step will be retried), negative for an unrecoverable error (in which case the integration is halted). @@ -178,8 +178,8 @@ class SundialsPdeScheduler : public Scheduler void applyMembraneFluxOperator(double t, double* yinput, double* rhs); void applyVolumeRegionReactionOperator(double t, double* yinput, double* rhs); void applyMembraneRegionReactionOperator(double t, double* yinput, double* rhs); - void initSundialsSolver(); - void solve(); + void initSundialsSolver(VCellModel* model); + void solve(SimTool* sim_tool); double *statePointValues, **neighborStatePointValues; void updateVolumeStatePointValues(int volIndex, double t, double* yinput, double* values); diff --git a/VCell/include/VCELL/VCellModel.h b/VCell/include/VCELL/VCellModel.h index 9c77c1f54..ad9d9c790 100644 --- a/VCell/include/VCELL/VCellModel.h +++ b/VCell/include/VCELL/VCellModel.h @@ -7,14 +7,12 @@ #include #include -#include using std::vector; using std::string; class Feature; class Membrane; -//class Contour; -class Simulation; +class SimulationExpression; struct MembraneElement; class VCellModel @@ -23,27 +21,25 @@ class VCellModel VCellModel(); ~VCellModel(); - //Contour *getContour(int type); - //int getNumContours(); - //void addContour(Contour *contour); - - int getNumFeatures() { - return (int)featureList.size(); + int getNumFeatures() const + { + return static_cast(featureList.size()); } Feature* addFeature(string& name, FeatureHandle handle); Feature* getFeatureFromHandle(FeatureHandle handle); Feature* getFeatureFromName(string& name); Feature* getFeatureFromIndex(int index); - int getNumMembranes() { - return (int)membraneList.size(); + int getNumMembranes() const + { + return static_cast(membraneList.size()); } Membrane* addMembrane(string& name, string& feature1_name, string& feature2_name); Membrane* getMembraneFromIndex(int index); Membrane* getMembraneFromName(string& mem_name); Membrane* getMembrane(Feature* f1, Feature* f2); - void resolveReferences(); + void resolveReferences(SimulationExpression* sim); bool hasFastSystem(); diff --git a/VCell/include/VCELL/VarContext.h b/VCell/include/VCELL/VarContext.h index 8b9e03e71..c7bcc9d23 100644 --- a/VCell/include/VCELL/VarContext.h +++ b/VCell/include/VCELL/VarContext.h @@ -7,6 +7,8 @@ #include #include + +#include "Expression.h" using std::string; using std::vector; @@ -30,7 +32,7 @@ static string String_Expression_Index[] = {"INITIAL_VALUE_EXP", "DIFF_RATE_EXP", class Variable; class VolumeVariable; class Structure; -class Simulation; +class SimulationExpression; class Mesh; class EqnBuilder; class SimulationExpression; @@ -44,62 +46,61 @@ namespace VCell { class VarContext { public: - ~VarContext(); + virtual ~VarContext(); - Variable *getVar() { return species; } + Variable *getVar() const { return species; } // // ALWAYS CALL ParentClass::resolveReferences() first // - virtual void resolveReferences(Simulation *sim); + virtual void resolveReferences(SimulationExpression *sim); virtual double getInitialValue(long index); virtual bool hasExact() { return false; } virtual bool hasRemainders() { return false; } - Structure *getStructure() { return structure; } + Structure *getStructure() const { return structure; } void setExpression(VCell::Expression* newexp, int expIndex); // exclusively for sundials pde - double evaluateJumpCondition(MembraneElement*, double* values); - double evaluateExpression(long expIndex, double* values); - double evaluateConstantExpression(long expIndex); + double evaluateJumpCondition(MembraneElement*, double* values) const; + double evaluateExpression(long expIndex, double* values) const; + double evaluateConstantExpression(long expIndex) const; - void addJumpCondition(Membrane* membrane, VCell::Expression* exp); - JumpCondition* getJumpCondition(); - void reinitConstantValues(); + void addJumpCondition(Membrane* membrane, VCell::Expression* exp); + void reinitConstantValues(SimulationExpression* sim); protected: VarContext(Structure *s, Variable* var); Variable *species; Structure *structure; - Simulation *sim; + SimulationExpression *sim; void bindAll(SimulationExpression* simulation); - double evaluateJumpCondition(MembraneElement*); - double evaluateExpression(long volIndex, long expIndex); // for volume - double evaluateMembraneRegionExpression(long memRegionIdex, long expIndex); // for membrane region - double evaluateVolumeRegionExpression(long volRegionIdex, long expIndex); // for volume regin - double evaluateExpression(MembraneElement* element, long expIndex); // for membrane - - virtual bool isNullExpressionOK(int expIndex) { return false; } - bool isConstantExpression(long expIndex); - bool isXYZOnlyExpression(long expIndex); - bool isNullExpression(int expIndex) { - return expressions[expIndex] == NULL; + double evaluateJumpCondition(MembraneElement*) const; + double evaluateExpression(long volIndex, long expIndex) const; // for volume + double evaluateMembraneRegionExpression(long memRegionIdex, long expIndex) const; // for membrane region + double evaluateVolumeRegionExpression(long volRegionIdex, long expIndex) const; // for volume regin + double evaluateExpression(MembraneElement* element, long expIndex) const; // for membrane + + virtual bool isNullExpressionOK(int expIndex) const = 0; + bool isConstantExpression(long expIndex) const; + bool isXYZOnlyExpression(long expIndex) const; + bool isNullExpression(int expIndex) const + { + return expressions[expIndex] == nullptr; } private: - VCell::Expression** expressions; - double** constantValues; - //bool* needsXYZ; - + vector expressions; + vector constantValues; + vector dependencyMask; vector jumpConditionList; - unsigned char* dependencyMask; - void computeDependencyMask(int expIndex); + + void computeDependencyMask(SimulationExpression* sim, int expIndex); }; #endif diff --git a/VCell/include/VCELL/VolumeRegionEqnBuilder.h b/VCell/include/VCELL/VolumeRegionEqnBuilder.h index 4cadbf904..e6306daf6 100644 --- a/VCell/include/VCELL/VolumeRegionEqnBuilder.h +++ b/VCell/include/VCELL/VolumeRegionEqnBuilder.h @@ -16,9 +16,9 @@ class VolumeRegionEqnBuilder : public EqnBuilder { public: VolumeRegionEqnBuilder(VolumeRegionVariable *species, CartesianMesh *mesh, ODESolver *solver); - void initEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) - {} - void buildEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize); + void initEquation(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) override + {} + void buildEquation(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) override; private: ODESolver* odeSolver; diff --git a/VCell/include/VCELL/VolumeRegionVarContextExpression.h b/VCell/include/VCELL/VolumeRegionVarContextExpression.h index b4263342e..eeacdeb89 100644 --- a/VCell/include/VCELL/VolumeRegionVarContextExpression.h +++ b/VCell/include/VCELL/VolumeRegionVarContextExpression.h @@ -14,17 +14,17 @@ struct MembraneElement; class VolumeRegionVarContextExpression : public VarContext { public: - void resolveReferences(Simulation *sim); + void resolveReferences(SimulationExpression *sim) override; - double getInitialValue(long index); - double getReactionRate(long volumeIndex); - double getUniformRate(VolumeRegion *region); - double getFlux(MembraneElement *element); + double getInitialValue(long index) override; + double getReactionRate(long volumeIndex) const; + double getUniformRate(VolumeRegion *region) const; + double getFlux(MembraneElement *element) const; VolumeRegionVarContextExpression(Feature *feature, VolumeRegionVariable* var); protected: - bool isNullExpressionOK(int expIndex); + bool isNullExpressionOK(int expIndex) const override; }; #endif diff --git a/VCell/include/VCELL/VolumeVarContextExpression.h b/VCell/include/VCELL/VolumeVarContextExpression.h index eca0cdad6..5be124f98 100644 --- a/VCell/include/VCELL/VolumeVarContextExpression.h +++ b/VCell/include/VCELL/VolumeVarContextExpression.h @@ -15,43 +15,43 @@ class VolumeVarContextExpression : public VarContext public: VolumeVarContextExpression(Feature *feature, VolumeVariable* var); - void resolveReferences(Simulation *sim); - - double getInitialValue(long index); - double getDiffusionRate(long index); - double getReactionRate(long volumeIndex); - double getFlux(MembraneElement *element); - - double getXmBoundaryValue(long index); - double getXpBoundaryValue(long index); - double getYmBoundaryValue(long index); - double getYpBoundaryValue(long index); - double getZmBoundaryValue(long index); - double getZpBoundaryValue(long index); - - double getXmBoundaryFlux(long index); - double getXpBoundaryFlux(long index); - double getYmBoundaryFlux(long index); - double getYpBoundaryFlux(long index); - double getZmBoundaryFlux(long index); - double getZpBoundaryFlux(long index); + void resolveReferences(SimulationExpression *sim) override; + + double getInitialValue(long index) override; + double getDiffusionRate(long index) const; + double getReactionRate(long volumeIndex) const; + double getFlux(MembraneElement *element) const; + + double getXmBoundaryValue(long index) const; + double getXpBoundaryValue(long index) const; + double getYmBoundaryValue(long index) const; + double getYpBoundaryValue(long index) const; + double getZmBoundaryValue(long index) const; + double getZpBoundaryValue(long index) const; + + double getXmBoundaryFlux(long index) const; + double getXpBoundaryFlux(long index) const; + double getYmBoundaryFlux(long index) const; + double getYpBoundaryFlux(long index) const; + double getZmBoundaryFlux(long index) const; + double getZpBoundaryFlux(long index) const; - double getXBoundaryPeriodicConstant(); - double getYBoundaryPeriodicConstant(); - double getZBoundaryPeriodicConstant(); + double getXBoundaryPeriodicConstant() const; + double getYBoundaryPeriodicConstant() const; + double getZBoundaryPeriodicConstant() const; - double getConvectionVelocity_X(long index); - double getConvectionVelocity_Y(long index); - double getConvectionVelocity_Z(long index); + double getConvectionVelocity_X(long index) const; + double getConvectionVelocity_Y(long index) const; + double getConvectionVelocity_Z(long index) const; - bool hasConstantDiffusion(); - bool hasConstantCoefficients(int dimension); - bool hasXYZOnlyDiffusion(); + bool hasConstantDiffusion() const; + bool hasConstantCoefficients(int dimension) const; + bool hasXYZOnlyDiffusion() const; - bool hasGradient(int dir); + bool hasGradient(int dir) const; private: - bool isNullExpressionOK(int expIndex); + bool isNullExpressionOK(int expIndex) const override; }; #endif diff --git a/VCell/src/CartesianMesh.cpp b/VCell/src/CartesianMesh.cpp index 399b2c22f..3f4f09c7e 100644 --- a/VCell/src/CartesianMesh.cpp +++ b/VCell/src/CartesianMesh.cpp @@ -17,21 +17,18 @@ using std::min; using std::set; #include -#include -#include +#include +#include #include #include #include #include #include -#include #include #include #include #include #include -#include -#include #include #include #include @@ -82,7 +79,22 @@ extern "C" //static const int NumNaturalNeighbors_Membrane[] = {0, 0, 2, 4}; -CartesianMesh::CartesianMesh(double AcaptureNeighborhood) : Mesh(AcaptureNeighborhood){ +CartesianMesh::CartesianMesh(double AcaptureNeighborhood) : Mesh(AcaptureNeighborhood), domainSizeX(0), domainSizeY(0), + domainSizeZ(0), + domainOriginX(0), + domainOriginY(0), + domainOriginZ(0), numX(0), + numY(0), numZ(0), + numXY(0), + scaleX_um(0), + scaleY_um(0), + scaleZ_um(0), + areaX_squm(0), + areaY_squm(0), + areaZ_squm(0), + volume_cu(0), + membraneInfo() +{ } //void CartesianMesh::setVolumeLists() @@ -160,19 +172,19 @@ CartesianMesh::CartesianMesh(double AcaptureNeighborhood) : Mesh(AcaptureNeighbo // } //} -void CartesianMesh::initialize(istream& ifs) +void CartesianMesh::initialize(VCellModel* model, istream& ifs) { - if (pVolumeElement!=NULL) { + if (pVolumeElement!=nullptr) { return; } - readGeometryFile(ifs); + readGeometryFile(model, ifs); initScale(); //sampleContours(); //setVolumeLists(); - printf("numVolume=%d\n",numVolume); + printf("numVolume=%ld\n",numVolume); setBoundaryConditions(); findMembraneNeighbors(); @@ -217,7 +229,7 @@ assert (alt == v); return (unsigned char)v; } */ -void CartesianMesh::readGeometryFile(istream& ifs) { +void CartesianMesh::readGeometryFile(VCellModel *model, istream& ifs) { stringstream ss; string line; @@ -246,6 +258,8 @@ void CartesianMesh::readGeometryFile(istream& ifs) { case 3: ss >> name >> domainSizeX >> domainSizeY >> domainSizeZ; break; + default: + throw std::runtime_error("CartesianMesh::readGeometryFile() : invalid dimension"); } //origin getline(ifs, line); @@ -270,7 +284,10 @@ void CartesianMesh::readGeometryFile(istream& ifs) { membraneInfo.nDirections = 4; membraneInfo.oppositeDirection= 2; break; + default: + throw std::runtime_error("CartesianMesh::readGeometryFile() : invalid dimension"); } + //volumeRegions getline(ifs, line); int numVolumeRegions; @@ -278,7 +295,6 @@ void CartesianMesh::readGeometryFile(istream& ifs) { ss.str(line); ss >> name >> numVolumeRegions; - VCellModel* model = SimTool::getInstance()->getModel(); for (int i = 0; i < numVolumeRegions; i ++) { int fhi; double volume; @@ -286,9 +302,9 @@ void CartesianMesh::readGeometryFile(istream& ifs) { ss.clear(); ss.str(line); ss >> name >> volume >> fhi; - FeatureHandle fh = (FeatureHandle)fhi; + auto fh = (FeatureHandle)fhi; Feature *feature = model->getFeatureFromHandle(fh); - VolumeRegion* vr = new VolumeRegion(i, name, this, feature); + auto* vr = new VolumeRegion(i, name, this, feature); vr->setSize(volume); feature->addRegion(vr); @@ -344,10 +360,10 @@ void CartesianMesh::readGeometryFile(istream& ifs) { //volumeSamples compressed, changed from byte to short int twiceNumVolume = 2 * numVolume; - unsigned char* bytes_from_compressed = new unsigned char[twiceNumVolume + 1000]; + auto* bytes_from_compressed = new unsigned char[twiceNumVolume + 1000]; memset(bytes_from_compressed, 0, (twiceNumVolume + 1000) * sizeof(unsigned char)); - unsigned char* compressed_hex = new unsigned char[twiceNumVolume + 1000]; + auto* compressed_hex = new unsigned char[twiceNumVolume + 1000]; memset(compressed_hex, 0, (twiceNumVolume + 1000) * sizeof(unsigned char)); ifs.getline((char*)compressed_hex, twiceNumVolume + 1000); const std::streamsize line_len = ifs.gcount() - 1; //don't count null termination character @@ -364,13 +380,13 @@ void CartesianMesh::readGeometryFile(istream& ifs) { bytes_from_compressed[j] = fromHex(compressed_hex + i); } - unsigned char* inflated_bytes = new unsigned char[twiceNumVolume + 1]; + auto* inflated_bytes = new unsigned char[twiceNumVolume + 1]; memset(inflated_bytes, 0, (twiceNumVolume + 1) * sizeof(unsigned char)); unsigned long inflated_len = twiceNumVolume; int retVal = uncompress(inflated_bytes, &inflated_len, bytes_from_compressed, compressed_len/2); - unsigned short* volsamples = new unsigned short[numVolume]; + auto* volsamples = new unsigned short[numVolume]; if (inflated_len == numVolume) { for (unsigned long i = 0; i < inflated_len; i ++) { volsamples[i] = inflated_bytes[i]; @@ -1263,7 +1279,7 @@ void CartesianMesh::findMembranePointInCurve(int N, long index, int neighbor_di } } if (!found) { - throw "CartesianMesh::findMembranePointInCurve(), new_neighbor can never be < 0"; + throw std::runtime_error("CartesianMesh::findMembranePointInCurve(), new_neighbor can never be < 0"); } return findMembranePointInCurve(N, neighborIndex, new_neighbor_dir, leftOverN, returnNeighbor); } @@ -1511,7 +1527,7 @@ void CartesianMesh::computeNormal(MembraneElement& meptr, const UnitVector3* nor case 3: break; default: - throw "CartesianMesh::computeNormal(), dimension should be 2 or 3."; + throw std::runtime_error("CartesianMesh::computeNormal(), dimension should be 2 or 3."); } // overwriting the normal computed by neighborhood if unnormalized length is much smaller than 1. diff --git a/VCell/src/EllipticVolumeEqnBuilder.cpp b/VCell/src/EllipticVolumeEqnBuilder.cpp index 00435aa4b..d097823e1 100644 --- a/VCell/src/EllipticVolumeEqnBuilder.cpp +++ b/VCell/src/EllipticVolumeEqnBuilder.cpp @@ -392,7 +392,7 @@ void EllipticVolumeEqnBuilder::computeLHS(int index, double& Aii, int& numCols, // Left Hand Side // //------------------------------------------------------------------ -void EllipticVolumeEqnBuilder::initEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) +void EllipticVolumeEqnBuilder::initEquation(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) { /** * we decided to scale both sides with deltaT/VOLUME @@ -406,7 +406,7 @@ void EllipticVolumeEqnBuilder::initEquation(double deltaTime, int volumeIndexSta * B_i = U_old * volumeScale_i + R * deltaT + boundary conditions **/ if (!bPreProcessed) { - preProcess(); + preProcess(model); } ASSERTION(solver->getVar() == var); @@ -494,10 +494,9 @@ void EllipticVolumeEqnBuilder::initEquation(double deltaTime, int volumeIndexSta delete[] columnValues; } -double EllipticVolumeEqnBuilder::computeRHS(int index) { +double EllipticVolumeEqnBuilder::computeRHS(Simulation* sim, int index) { string varname = var->getName(); double b = 0; - Simulation* sim = SimTool::getInstance()->getSimulation(); VolumeElement *pVolumeElement = mesh->getVolumeElements(); MembraneElement *pMembraneElement = mesh->getMembraneElements(); @@ -640,16 +639,16 @@ double EllipticVolumeEqnBuilder::computeRHS(int index) { // Right Hand side // //------------------------------------------------------------------ -void EllipticVolumeEqnBuilder::buildEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) +void EllipticVolumeEqnBuilder::buildEquation(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) { if (bSolveWholeMesh) { for (int index = volumeIndexStart; index < volumeIndexStart + volumeIndexSize; index ++){ - B[index] = computeRHS(index); + B[index] = computeRHS(sim, index); } } else { for (int localIndex = 0; localIndex < getSize() ; localIndex ++) { int globalIndex = LocalToGlobalMap[localIndex]; - B[localIndex] = computeRHS(globalIndex); + B[localIndex] = computeRHS(sim, globalIndex); } } // to make the matrix symmetric @@ -684,7 +683,7 @@ bool EllipticVolumeEqnBuilder::checkPeriodicCoupledPairsInRegions(int indexm, in return false; // both of them are not in the solved regions, so we don't have to add them to the list } -void EllipticVolumeEqnBuilder::preProcess() { +void EllipticVolumeEqnBuilder::preProcess(VCellModel *model) { if (bPreProcessed) { return; } @@ -693,8 +692,8 @@ void EllipticVolumeEqnBuilder::preProcess() { // check if there is periodic boundary condition in the model bool bPeriodic = false; - for (int i = 0; i < SimTool::getInstance()->getModel()->getNumFeatures(); i ++) { - Feature* feature = SimTool::getInstance()->getModel()->getFeatureFromIndex(i); + for (int i = 0; i < model->getNumFeatures(); i ++) { + Feature* feature = model->getFeatureFromIndex(i); if (feature->getXmBoundaryType() == BOUNDARY_PERIODIC || feature->getYmBoundaryType() == BOUNDARY_PERIODIC || feature->getZmBoundaryType() == BOUNDARY_PERIODIC) { diff --git a/VCell/src/EqnBuilderReactionForward.cpp b/VCell/src/EqnBuilderReactionForward.cpp index c328bdd92..3e9b2e89f 100644 --- a/VCell/src/EqnBuilderReactionForward.cpp +++ b/VCell/src/EqnBuilderReactionForward.cpp @@ -18,14 +18,13 @@ EqnBuilderReactionForward::EqnBuilderReactionForward(VolumeVariable *Avar, Mesh odeSolver = Asolver; } -void EqnBuilderReactionForward::buildEquation(double deltaTime, +void EqnBuilderReactionForward::buildEquation(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) { Feature *feature; VolumeVarContextExpression *varContext; - Simulation *sim = SimTool::getInstance()->getSimulation(); long arraySize = odeSolver->getArraySize(); if(arraySize==0){ ASSERTION((volumeIndexStart>=0) && ((volumeIndexStart+volumeIndexSize)<=mesh->getNumVolumeElements())); diff --git a/VCell/src/FVDataSet.cpp b/VCell/src/FVDataSet.cpp index 457f59955..79fab77af 100644 --- a/VCell/src/FVDataSet.cpp +++ b/VCell/src/FVDataSet.cpp @@ -25,7 +25,7 @@ using std::endl; #define CONVOLVE_SUFFIX "_Convolved" -FieldData* getPSFFieldData(); +FieldData* getPSFFieldData(const SimulationExpression* sim); void FVDataSet::readRandomVariables(char *filename, SimulationExpression *sim) { @@ -161,9 +161,9 @@ void FVDataSet::read(const char *filename, Simulation *sim) fclose(fp); } -void FVDataSet::convolve(Simulation* sim, Variable* var, double* values) { +void FVDataSet::convolve(SimulationExpression* sim, Variable* var, double* values) { - FieldData* psfFieldData = getPSFFieldData(); + FieldData* psfFieldData = getPSFFieldData(sim); if (psfFieldData == 0) { throw "psf field data is not defined"; } @@ -293,7 +293,7 @@ void FVDataSet::write(const char *filename, SimulationExpression *sim, bool bCom cout << "DataSet::write() - no variables defined" << endl; } - FieldData* psfFieldData = getPSFFieldData(); + FieldData* psfFieldData = getPSFFieldData(sim); int numBlocks = psfFieldData == 0 ? numVars : numVars*2; // region size variable diff --git a/VCell/src/FVSolver.cpp b/VCell/src/FVSolver.cpp index 0a5da4901..22340064a 100644 --- a/VCell/src/FVSolver.cpp +++ b/VCell/src/FVSolver.cpp @@ -1,8 +1,8 @@ #include -#include +#include #include -#include +#include #include #include using std::ifstream; @@ -55,8 +55,8 @@ using VCell::Expression; #include #include -FieldData *getPSFFieldData() { - return ((SimulationExpression*)SimTool::getInstance()->getSimulation())->getPSFFieldData(); +FieldData *getPSFFieldData(const SimulationExpression* sim) { + return sim->getPSFFieldData(); } /* @@ -87,7 +87,7 @@ void FVSolver::loadJMSInfo(istream& ifsInput, int taskID) { nextToken = ""; lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "JMS_PARAM_END") { @@ -136,9 +136,9 @@ FEATURE ec 1 value value value value MEMBRANE cyt_ec_membrane cyt ec value value value value MODEL_END */ -void FVSolver::loadModel(istream& ifsInput) { +VCellModel* FVSolver::loadModel(istream& ifsInput) { //cout << "loading model " << endl; - model = new VCellModel(); + auto* model = new VCellModel(); string nextToken, line; string feature_name; int handle; @@ -150,14 +150,14 @@ void FVSolver::loadModel(istream& ifsInput) { nextToken = ""; lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "MODEL_END") { break; } - Structure* structure = 0; + Structure* structure = nullptr; if (nextToken == "FEATURE") { numFeatures ++; lineInput >> feature_name >> handle; @@ -177,7 +177,7 @@ void FVSolver::loadModel(istream& ifsInput) { break; } lineInput >> btstr; - if (btstr.length() == 0) { + if (btstr.empty()) { break; } BoundaryType bt = BOUNDARY_VALUE; @@ -190,7 +190,7 @@ void FVSolver::loadModel(istream& ifsInput) { } else { stringstream ss1; ss1 << "loadModel(), wrong boundary type " << btstr; - throw ss1.str(); + throw runtime_error(ss1.str()); } switch (i) { case 0: // XM @@ -214,26 +214,28 @@ void FVSolver::loadModel(istream& ifsInput) { } } } + return model; } -int FVSolver::loadSolveRegions(istream& lineInput, int*& solveRegions) { +int FVSolver::loadSolveRegions(VCellModel* model, istream& lineInput, int*& solveRegions) const +{ int numVolumeRegions = mesh->getNumVolumeRegions(); int regionCount = 0; while (true) { - string feature_name = ""; + string feature_name; lineInput >> feature_name; - if (feature_name == "") { + if (feature_name.empty()) { break; } for (int i = 0; i < numVolumeRegions; i++){ VolumeRegion *volRegion = mesh->getVolumeRegion(i); - Feature* feature = SimTool::getInstance()->getModel()->getFeatureFromName(feature_name); - if (feature == NULL) { + Feature* feature = model->getFeatureFromName(feature_name); + if (feature == nullptr) { stringstream ss; ss << "Feature '" << feature_name << "' doesn't exist!"; - throw ss.str(); + throw runtime_error(ss.str()); } - if (solveRegions == 0) { + if (solveRegions == nullptr) { solveRegions = new int[numVolumeRegions]; } if (volRegion->getFeature()->getHandle() == feature->getHandle()){ @@ -255,9 +257,9 @@ VOLUME_PDE Na false false true MEMBRANE_REGION Voltage_membrane VARIABLE_END */ -void FVSolver::loadSimulation(istream& ifsInput) { +SimulationExpression* FVSolver::loadSimulation(SimTool* simTool, VCellModel *model, istream& ifsInput) { //cout << "loading simulation" << endl; - simulation = new SimulationExpression(mesh); + auto* simulation = new SimulationExpression(mesh); string nextToken, line; long sizeX = mesh->getNumVolumeX(); long sizeY = mesh->getNumVolumeY(); @@ -271,7 +273,7 @@ void FVSolver::loadSimulation(istream& ifsInput) { nextToken = ""; lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "VARIABLE_END") { @@ -281,12 +283,12 @@ void FVSolver::loadSimulation(istream& ifsInput) { if (nextToken == "VOLUME_RANDOM") { string name; ifsInput >> name; - RandomVariable* rv = new RandomVariable(name, VAR_VOLUME, mesh->getNumVolumeElements()); + auto* rv = new RandomVariable(name, VAR_VOLUME, mesh->getNumVolumeElements()); simulation->addRandomVariable(rv); } else if (nextToken == "MEMBRANE_RANDOM") { string name; ifsInput >> name; - RandomVariable* rv = new RandomVariable(name, VAR_MEMBRANE, mesh->getNumMembraneElements()); + auto* rv = new RandomVariable(name, VAR_MEMBRANE, mesh->getNumMembraneElements()); simulation->addRandomVariable(rv); } else if (nextToken == "VOLUME_PDE" || nextToken == "VOLUME_PDE_STEADY") { bool bSteady = false; @@ -304,7 +306,7 @@ void FVSolver::loadSimulation(istream& ifsInput) { bool bSolveVariable = true; int numSolveRegions = 0; // flag specifying to solve for all regions - int *solveRegions = 0; + int *solveRegions = nullptr; bool* vrmap = new bool[numVolumeRegions]; for (int i = 0; i < numVolumeRegions; i ++) { @@ -313,7 +315,7 @@ void FVSolver::loadSimulation(istream& ifsInput) { if (solve_whole_mesh_flag == "false") { memset(vrmap, 0, numVolumeRegions * sizeof(bool)); // if not solve everywhere; - numSolveRegions = loadSolveRegions(lineInput, solveRegions); + numSolveRegions = loadSolveRegions(model, lineInput, solveRegions); if (numSolveRegions == 0) { bSolveVariable = false; @@ -334,9 +336,9 @@ void FVSolver::loadSimulation(istream& ifsInput) { bool bGradient = gradient_flag == "true"; Feature* feature = model->getFeatureFromName(structure_name); - VolumeVariable* volumeVar = new VolumeVariable(variable_name, feature, sizeX, sizeY, sizeZ, true, !bNoConvection, bGradient); + auto* volumeVar = new VolumeVariable(variable_name, feature, sizeX, sizeY, sizeZ, true, !bNoConvection, bGradient); if (bSolveVariable && !simTool->isSundialsPdeSolver()) { - SparseMatrixEqnBuilder* builder = 0; + SparseMatrixEqnBuilder* builder = nullptr; if (bSteady) { builder = new EllipticVolumeEqnBuilder(volumeVar,mesh, numSolveRegions, solveRegions); } else { @@ -353,7 +355,7 @@ void FVSolver::loadSimulation(istream& ifsInput) { assert(solve_whole_mesh_flag == "true" || solve_whole_mesh_flag == "false"); int numSolveRegions = 0; // flag specifying to solve for all regions - int *solveRegions = 0; + int *solveRegions = nullptr; bool* vrmap = new bool[numVolumeRegions]; for (int i = 0; i < numVolumeRegions; i ++) { @@ -363,7 +365,7 @@ void FVSolver::loadSimulation(istream& ifsInput) { bool bSolveVariable = true; if (solve_whole_mesh_flag == "false") { memset(vrmap, 0, numVolumeRegions * sizeof(bool)); // if not solve everywhere; - numSolveRegions = loadSolveRegions(lineInput, solveRegions); + numSolveRegions = loadSolveRegions(model, lineInput, solveRegions); if (numSolveRegions == 0) { bSolveVariable = false; @@ -375,9 +377,9 @@ void FVSolver::loadSimulation(istream& ifsInput) { } Feature* feature = model->getFeatureFromName(structure_name); - VolumeVariable* volumeVar = new VolumeVariable(variable_name, feature, sizeX, sizeY, sizeZ, false); + auto* volumeVar = new VolumeVariable(variable_name, feature, sizeX, sizeY, sizeZ, false); if (bSolveVariable && !simTool->isSundialsPdeSolver()) { - ODESolver* odeSolver = new ODESolver(volumeVar,mesh,numSolveRegions,solveRegions); + auto* odeSolver = new ODESolver(volumeVar,mesh,numSolveRegions,solveRegions); EqnBuilder* builder = new EqnBuilderReactionForward(volumeVar,mesh,odeSolver); odeSolver->setEqnBuilder(builder); simulation->addSolver(odeSolver); @@ -386,22 +388,22 @@ void FVSolver::loadSimulation(istream& ifsInput) { } else if (nextToken == "VOLUME_PARTICLE") { lineInput >> variable_name >> structure_name; Feature* feature = model->getFeatureFromName(structure_name); - VolumeParticleVariable* volumeParticleVar = new VolumeParticleVariable(variable_name, feature, sizeX, sizeY, sizeZ); + auto* volumeParticleVar = new VolumeParticleVariable(variable_name, feature, sizeX, sizeY, sizeZ); simulation->addVolumeParticleVariable(volumeParticleVar); } else if (nextToken == "MEMBRANE_PARTICLE") { lineInput >> variable_name >> structure_name; Membrane* membrane = model->getMembraneFromName(structure_name); - MembraneParticleVariable* membraneParticleVar = new MembraneParticleVariable(variable_name, membrane, mesh->getNumMembraneElements()); + auto* membraneParticleVar = new MembraneParticleVariable(variable_name, membrane, mesh->getNumMembraneElements()); simulation->addMembraneParticleVariable(membraneParticleVar); } else if (nextToken == "MEMBRANE_ODE") { lineInput >> variable_name >> structure_name; int numSolveRegions = 0; // flag specifying to solve for all regions - int *solveRegions = NULL; + int *solveRegions = nullptr; Membrane* membrane = model->getMembraneFromName(structure_name); MembraneVariable* membraneVar = new MembraneVariable(variable_name, membrane, mesh->getNumMembraneElements(), false); if (!simTool->isSundialsPdeSolver()) { - ODESolver* odeSolver = new ODESolver(membraneVar,mesh,numSolveRegions,solveRegions); + auto* odeSolver = new ODESolver(membraneVar,mesh,numSolveRegions,solveRegions); EqnBuilder* builder = new MembraneEqnBuilderForward(membraneVar,mesh,odeSolver); odeSolver->setEqnBuilder(builder); simulation->addSolver(odeSolver); @@ -418,22 +420,22 @@ void FVSolver::loadSimulation(istream& ifsInput) { } Membrane* membrane = model->getMembraneFromName(structure_name); - MembraneVariable* membraneVar = new MembraneVariable(variable_name, membrane, mesh->getNumMembraneElements(), true); + auto* membraneVar = new MembraneVariable(variable_name, membrane, mesh->getNumMembraneElements(), true); if (!simTool->isSundialsPdeSolver()) { SparseMatrixEqnBuilder* smbuilder = new MembraneEqnBuilderDiffusion(membraneVar,mesh); - SparseLinearSolver* slSolver = new SparseLinearSolver(membraneVar,smbuilder,simTool->getPCGRelativeErrorTolerance(),bTimeDependent); + auto* slSolver = new SparseLinearSolver(membraneVar,smbuilder,simTool->getPCGRelativeErrorTolerance(),bTimeDependent); simulation->addSolver(slSolver); } simulation->addMembraneVariable(membraneVar); } else if (nextToken == "VOLUME_REGION") { lineInput >> variable_name >> structure_name; int numSolveRegions = 0; // flag specifying to solve for all regions - int *solveRegions = NULL; + int *solveRegions = nullptr; Feature* feature = model->getFeatureFromName(structure_name); - VolumeRegionVariable* volumeRegionVar = new VolumeRegionVariable(variable_name, feature, mesh->getNumVolumeRegions()); + auto* volumeRegionVar = new VolumeRegionVariable(variable_name, feature, mesh->getNumVolumeRegions()); if (!simTool->isSundialsPdeSolver()) { - ODESolver* odeSolver = new ODESolver(volumeRegionVar,mesh,numSolveRegions,solveRegions); + auto* odeSolver = new ODESolver(volumeRegionVar,mesh,numSolveRegions,solveRegions); EqnBuilder* builder = new VolumeRegionEqnBuilder(volumeRegionVar,mesh,odeSolver); odeSolver->setEqnBuilder(builder); simulation->addSolver(odeSolver); @@ -442,12 +444,12 @@ void FVSolver::loadSimulation(istream& ifsInput) { } else if (nextToken == "MEMBRANE_REGION") { lineInput >> variable_name >> structure_name; int numSolveRegions = 0; // flag specifying to solve for all regions - int *solveRegions = NULL; + int *solveRegions = nullptr; Membrane* membrane = model->getMembraneFromName(structure_name); - MembraneRegionVariable* memRegionVariable = new MembraneRegionVariable(variable_name, membrane, mesh->getNumMembraneRegions()); + auto* memRegionVariable = new MembraneRegionVariable(variable_name, membrane, mesh->getNumMembraneRegions()); if (!simTool->isSundialsPdeSolver()) { - ODESolver* odeSolver = new ODESolver(memRegionVariable,mesh,numSolveRegions,solveRegions); + auto* odeSolver = new ODESolver(memRegionVariable,mesh,numSolveRegions,solveRegions); EqnBuilder* builder = new MembraneRegionEqnBuilder(memRegionVariable,mesh,odeSolver); odeSolver->setEqnBuilder(builder); simulation->addSolver(odeSolver); @@ -455,9 +457,10 @@ void FVSolver::loadSimulation(istream& ifsInput) { simulation->addMembraneRegionVariable(memRegionVariable); } } + return simulation; } -Expression* FVSolver::readExpression(istream& lineInput, string& var_name, string prefix) { +Expression* FVSolver::readExpression(istream& lineInput, string& var_name, const string& prefix) { string expStr; getline(lineInput, expStr); expStr = prefix + expStr; @@ -465,7 +468,7 @@ Expression* FVSolver::readExpression(istream& lineInput, string& var_name, strin if (expStr[expStr.size()-1] != ';') { stringstream msg; msg << "Expression for [" << var_name << "] is not terminated by ';'"; - throw msg.str(); + throw runtime_error(msg.str()); } return new Expression(expStr); } @@ -487,7 +490,7 @@ VarContext* FVSolver::loadEquation(istream& ifsInput, Structure* structure, Vari //cout << "loading volume var context " << var_name << endl; - VarContext* varContext = NULL; + VarContext* varContext = nullptr; if (var->getVarType() == VAR_VOLUME) { varContext = new VolumeVarContextExpression((Feature*)structure, (VolumeVariable*)var); } else if (var->getVarType() == VAR_VOLUME_REGION) { @@ -499,7 +502,7 @@ VarContext* FVSolver::loadEquation(istream& ifsInput, Structure* structure, Vari } else { stringstream ss; ss << "loadEquation: variable type not supported yet: " << var->getName(); - throw ss.str(); + throw runtime_error(ss.str()); } string nextToken, line; @@ -510,7 +513,7 @@ VarContext* FVSolver::loadEquation(istream& ifsInput, Structure* structure, Vari nextToken = ""; lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "EQUATION_END") { @@ -568,7 +571,7 @@ FLUX Nucleus (150.0 * (Ca_Cytosol_membrane - Ca_Nucleus_membrane)); FLUX Cytosol - (150.0 * (Ca_Cytosol_membrane - Ca_Nucleus_membrane)); JUMP_CONDITION_END */ -void FVSolver::loadJumpCondition(istream& ifsInput, Membrane* membrane, string& var_name) { +void FVSolver::loadJumpCondition(SimulationExpression* simulation, VCellModel* model, istream& ifsInput, Membrane* membrane, string& var_name) { //cout << "loading jump condition " << var_name << endl; string nextToken, line; @@ -579,7 +582,7 @@ void FVSolver::loadJumpCondition(istream& ifsInput, Membrane* membrane, string& nextToken = ""; lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "JUMP_CONDITION_END") { @@ -592,22 +595,22 @@ void FVSolver::loadJumpCondition(istream& ifsInput, Membrane* membrane, string& Feature* f = model->getFeatureFromName(featurename); string var_name = var->getName(); Expression* exp = readExpression(lineInput, var_name); - VarContext *varContext = 0; + VarContext *varContext = nullptr; if (var->getVarType() == VAR_VOLUME) { varContext = f->getVolumeVarContext((VolumeVariable*)var); } else if (var->getVarType() == VAR_VOLUME_REGION) { varContext = f->getVolumeRegionVarContext((VolumeRegionVariable*)var); } else { - throw "Only volume variables and volume region variables have jump conditions"; + throw runtime_error("Only volume variables and volume region variables have jump conditions"); } - if (varContext == NULL) { + if (varContext == nullptr) { stringstream ss; ss << "variable '" << var->getName() << "' is not defined in " << featurename; - throw ss.str(); + throw runtime_error(ss.str()); } varContext->addJumpCondition(membrane, exp); } else { - throw "Expecting FLUX in JumpCondition."; + throw runtime_error("Expecting FLUX in JumpCondition."); } } } @@ -623,8 +626,8 @@ void FVSolver::loadPseudoConstants(istream& ifsInput, FastSystemExpression* fast string nextToken, line; int count = 0; int numDep = fastSystem->getNumDependents(); - string* vars = new string[numDep]; - Expression **expressions = new Expression*[numDep]; + auto* vars = new string[numDep]; + auto **expressions = new Expression*[numDep]; while (!ifsInput.eof()) { getline(ifsInput, line); @@ -632,7 +635,7 @@ void FVSolver::loadPseudoConstants(istream& ifsInput, FastSystemExpression* fast nextToken = ""; lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "PSEUDO_CONSTANT_END") { @@ -660,7 +663,7 @@ void FVSolver::loadFastRates(istream& ifsInput, FastSystemExpression* fastSystem string nextToken, line; int count = 0; int numIndep = fastSystem->getDimension(); - Expression **expressions = new Expression*[numIndep]; + auto **expressions = new Expression*[numIndep]; while (!ifsInput.eof()) { getline(ifsInput, line); @@ -668,7 +671,7 @@ void FVSolver::loadFastRates(istream& ifsInput, FastSystemExpression* fastSystem nextToken = ""; lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "FAST_RATE_END") { @@ -680,7 +683,7 @@ void FVSolver::loadFastRates(istream& ifsInput, FastSystemExpression* fastSystem count ++; } if (count != numIndep) { - throw "In the fast system the number of fast rates should be the same as that of independent variables"; + throw runtime_error("In the fast system the number of fast rates should be the same as that of independent variables"); } fastSystem->setFastRateExpressions(expressions); } @@ -696,8 +699,8 @@ void FVSolver::loadFastDependencies(istream& ifsInput, FastSystemExpression* fas string nextToken, line; int count = 0; int numDep = fastSystem->getNumDependents(); - string* vars = new string[numDep]; - Expression **expressions = new Expression*[numDep]; + auto* vars = new string[numDep]; + auto **expressions = new Expression*[numDep]; while (!ifsInput.eof()) { getline(ifsInput, line); @@ -705,7 +708,7 @@ void FVSolver::loadFastDependencies(istream& ifsInput, FastSystemExpression* fas nextToken = ""; lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "FAST_DEPENDENCY_END") { @@ -719,7 +722,7 @@ void FVSolver::loadFastDependencies(istream& ifsInput, FastSystemExpression* fas fastSystem->setFastDependencyExpressions(vars, expressions); delete[] vars; if (count != numDep) { - throw "In the fast system the number of fast dependencies should be the same as that of dependent variables"; + throw runtime_error("In the fast system the number of fast dependencies should be the same as that of dependent variables"); } } @@ -736,7 +739,7 @@ void FVSolver::loadJacobians(istream& ifsInput, FastSystemExpression* fastSystem string nextToken, line; int count = 0; int numIndep = fastSystem->getDimension(); - Expression **expressions = new Expression*[numIndep * numIndep]; + auto **expressions = new Expression*[numIndep * numIndep]; while (!ifsInput.eof()) { getline(ifsInput, line); @@ -744,7 +747,7 @@ void FVSolver::loadJacobians(istream& ifsInput, FastSystemExpression* fastSystem nextToken = ""; lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "JACOBIAN_END") { @@ -756,7 +759,7 @@ void FVSolver::loadJacobians(istream& ifsInput, FastSystemExpression* fastSystem count ++; } if (count != numIndep * numIndep) { - throw "In the fast system the number of Jacobian should dim*dim"; + throw runtime_error("In the fast system the number of Jacobian should dim*dim"); } fastSystem->setJacobianExpressions(expressions); } @@ -803,7 +806,7 @@ void FVSolver::loadFastSystem(istream& ifsInput, FastSystemExpression* fastSyste nextToken = ""; lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "FAST_SYSTEM_END") { @@ -811,14 +814,14 @@ void FVSolver::loadFastSystem(istream& ifsInput, FastSystemExpression* fastSyste } if (nextToken == "DEPENDENT_VARIALBES") { - string* vars = new string[numDep]; + auto* vars = new string[numDep]; for (int i = 0; i < numDep; i ++) { lineInput >> vars[i]; } fastSystem->setDependentVariables(vars); delete[] vars; } else if (nextToken == "INDEPENDENT_VARIALBES") { - string* vars = new string[numIndep]; + auto* vars = new string[numIndep]; for (int i = 0; i < numIndep; i ++) { lineInput >> vars[i]; } @@ -873,7 +876,7 @@ EQUATION_END COMPARTMENT_END */ -void FVSolver::loadFeature(istream& ifsInput, Feature* feature) { +void FVSolver::loadFeature(SimulationExpression* simulation, istream& ifsInput, Feature* feature) { //cout << "loading feature " << feature->getName() << endl; string nextToken, line; @@ -883,7 +886,7 @@ void FVSolver::loadFeature(istream& ifsInput, Feature* feature) { nextToken = ""; lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "COMPARTMENT_END") { @@ -903,7 +906,7 @@ void FVSolver::loadFeature(istream& ifsInput, Feature* feature) { } else if (nextToken == "FAST_SYSTEM_BEGIN") { int dimension, num_of_dependents; lineInput >> dimension >> num_of_dependents; - FastSystemExpression* fastSystem = new FastSystemExpression(dimension, num_of_dependents, simulation); + auto* fastSystem = new FastSystemExpression(dimension, num_of_dependents, simulation); loadFastSystem(ifsInput, fastSystem); feature->setFastSystem(fastSystem); } @@ -948,7 +951,7 @@ JUMP_CONDITION_END MEMBRANE_END */ -void FVSolver::loadMembrane(istream& ifsInput, Membrane* membrane) { +void FVSolver::loadMembrane(SimulationExpression* simulation, VCellModel* model, istream& ifsInput, Membrane* membrane) { //cout << "loading membrane " << var_name << endl; string nextToken, line; @@ -958,7 +961,7 @@ void FVSolver::loadMembrane(istream& ifsInput, Membrane* membrane) { nextToken = ""; lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "MEMBRANE_END") { @@ -978,13 +981,13 @@ void FVSolver::loadMembrane(istream& ifsInput, Membrane* membrane) { } else if (nextToken == "FAST_SYSTEM_BEGIN") { int dimension, num_of_dependents; lineInput >> dimension >> num_of_dependents; - FastSystemExpression* fastSystem = new FastSystemExpression(dimension, num_of_dependents, simulation); + auto* fastSystem = new FastSystemExpression(dimension, num_of_dependents, simulation); loadFastSystem(ifsInput, fastSystem); membrane->setFastSystem(fastSystem); } else if (nextToken == "JUMP_CONDITION_BEGIN") { string var_name; lineInput >> var_name; - loadJumpCondition(ifsInput, membrane, var_name); + loadJumpCondition(simulation, model, ifsInput, membrane, var_name); } } } @@ -1008,7 +1011,8 @@ TIME_STEP 0.01 KEEP_EVERY 10 SIMULATION_PARAM_END */ -void FVSolver::loadSimulationParameters(istream& ifsInput) { +void FVSolver::loadSimulationParameters(SimTool* simTool, istream& ifsInput) const +{ string nextToken, line; while (!ifsInput.eof()) { @@ -1017,7 +1021,7 @@ void FVSolver::loadSimulationParameters(istream& ifsInput) { nextToken = ""; lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "SIMULATION_PARAM_END") { @@ -1025,7 +1029,7 @@ void FVSolver::loadSimulationParameters(istream& ifsInput) { } if (nextToken == "SOLVER") { - string solver=""; + string solver; lineInput >> solver; simTool->setSolver(solver); if (solver == FV_SOLVER) { @@ -1064,7 +1068,7 @@ void FVSolver::loadSimulationParameters(istream& ifsInput) { int numDisTimes = 0; lineInput >> numDisTimes; if (numDisTimes > 0) { - double* discontinuityTimes = 0; + double* discontinuityTimes = nullptr; discontinuityTimes = new double[numDisTimes]; for (int i = 0; i < numDisTimes; i ++) { lineInput >> discontinuityTimes[i]; @@ -1075,11 +1079,11 @@ void FVSolver::loadSimulationParameters(istream& ifsInput) { string basefilename; getline(lineInput, basefilename); trimString(basefilename); - if (outputPath == 0) { + if (outputPath == nullptr) { simTool->setBaseFilename((char*)basefilename.c_str()); } else { const char* baseSimName = strrchr(basefilename.c_str(), DIRECTORY_SEPARATOR); - if (baseSimName == NULL) { + if (baseSimName == nullptr) { baseSimName = basefilename.c_str(); } else { baseSimName += 1; @@ -1127,9 +1131,7 @@ void FVSolver::loadSimulationParameters(istream& ifsInput) { lineInput >> bStoreEnable; simTool->setStoreEnable(bStoreEnable!=0); } else { - stringstream ss; - ss << "loadSimulationParameters(), encountered unknown token " << nextToken << endl; - throw ss.str(); + throw runtime_error("loadSimulationParameters(), encountered unknown token " + nextToken); } } @@ -1142,65 +1144,14 @@ MESH_BEGIN VCG_FILE \\cfs01.vcell.uchc.edu\raid\Vcell\users\fgao\SimID_36230826_0_.vcg MESH_END */ -void FVSolver::loadMesh(istream& ifsInput) { - if (SimTool::getInstance()->getModel() == 0) { - throw "Model has to be initialized before mesh initialization"; - } - - string meshfile = ""; - string nextToken, line; - string vcgText = ""; - - while (!ifsInput.eof()) { - getline(ifsInput, line); - istringstream lineInput(line); - - nextToken = ""; - lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { - continue; - } - if (nextToken == "MESH_END") { - break; - } - - if (nextToken == "VCG_FILE") { - getline(lineInput, meshfile); - trimString(meshfile); - struct stat buf; - if (stat(meshfile.c_str(), &buf)) { - stringstream ss; - ss << "Mesh file(.vcg) [" << meshfile <<"] doesn't exist"; - throw ss.str(); - } - } else { // VCG In file - vcgText += nextToken; - getline(lineInput, nextToken); - vcgText += nextToken + "\n"; - } +void FVSolver::loadMeshFromVcg(VCellModel* model, istream& vcgInput) { + if (model == nullptr) { + throw std::runtime_error("Model has to be initialized before mesh initialization"); } SimulationMessaging::getInstVar()->setWorkerEvent(new WorkerEvent(JOB_STARTING, "initializing mesh")); mesh = new CartesianMesh(); - if (meshfile.size() != 0) { - ifstream ifs(meshfile.c_str()); - if (!ifs.is_open()){ - stringstream ss; - ss << "Can't open geometry file '" << meshfile << "'"; - throw ss.str(); - } - cout << "Reading mesh from vcg file '" << meshfile << "'" << endl; - mesh->initialize(ifs); - ifs.close(); - } else { - if (vcgText.size() == 0) { - throw "no mesh specified"; - } - cout << "Reading mesh from text..." << endl; - //cout << vcgText << endl; - istringstream iss(vcgText); - mesh->initialize(iss); - } + mesh->initialize(model, vcgInput); SimulationMessaging::getInstVar()->setWorkerEvent(new WorkerEvent(JOB_STARTING, "mesh initialized")); } @@ -1211,9 +1162,9 @@ FIELD_DATA_BEGIN 0 _VCell_FieldData_0 FRAP_binding_ALPHA rfB 0.1 \\\\SAN2\\raid\\Vcell\\users\\fgao\\SimID_22489731_0_FRAP_binding_ALPHA_rfB_0_1.fdat FIELD_DATA_END */ -void FVSolver::loadFieldData(istream& ifsInput) { - if (simulation == 0) { - throw "Simulation has to be initialized before loading field data"; +void FVSolver::loadFieldData(SimulationExpression* simulation, istream& ifsInput) { + if (simulation == nullptr) { + throw runtime_error("Simulation has to be initialized before loading field data"); } string nextToken, line; @@ -1227,7 +1178,7 @@ void FVSolver::loadFieldData(istream& ifsInput) { nextToken = ""; lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "FIELD_DATA_END") { @@ -1245,8 +1196,8 @@ void FVSolver::loadFieldData(istream& ifsInput) { lineInput >> fdVarType >> fdID >> fdName >> fdVarName >> fdTime; getline(lineInput, fdFile); trimString(fdFile); - if (fdVarType == "" || fdID == "" || fdName == "" || fdVarName == "" || fdFile == "" || fdIndex < 0 || fdTime < 0) { - throw "loadFieldData(), wrong input"; + if (fdVarType.empty() || fdID.empty() || fdName.empty() || fdVarName.empty() || fdFile.empty() || fdIndex < 0 || fdTime < 0) { + throw runtime_error("loadFieldData(), wrong input"); } VariableType varType = VAR_UNKNOWN; if (fdVarType == "Membrane") { @@ -1254,7 +1205,7 @@ void FVSolver::loadFieldData(istream& ifsInput) { } else if (fdVarType == "Volume") { varType = VAR_VOLUME; } else { - throw "field data is only supported for volume and membrane variables"; + throw runtime_error("field data is only supported for volume and membrane variables"); } simulation->addFieldData(new FieldData(fdIndex, varType, fdID, fdName, fdVarName, fdTime, fdFile)); } @@ -1269,9 +1220,9 @@ U0 U1 PARAMETER_END */ -void FVSolver::loadParameters(istream& ifsInput, int numParameters) { - if (simulation == 0) { - throw "Simulation has to be initialized before loading field data"; +void FVSolver::loadParameters(SimulationExpression* simulation, istream& ifsInput, int numParameters) { + if (simulation == nullptr) { + throw runtime_error("Simulation has to be initialized before loading field data"); } string nextToken, line; @@ -1281,7 +1232,7 @@ void FVSolver::loadParameters(istream& ifsInput, int numParameters) { istringstream lineInput(line); lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "PARAMETER_END") { @@ -1294,9 +1245,9 @@ void FVSolver::loadParameters(istream& ifsInput, int numParameters) { assert(nread == numParameters); } -void FVSolver::loadSerialScanParameters(istream& ifsInput, int numSerialScanParameters) { - if (simulation == 0) { - throw "Simulation has to be initialized before loading serial scan parameters"; +void FVSolver::loadSerialScanParameters(SimulationExpression* simulation, istream& ifsInput, int numSerialScanParameters) { + if (simulation == nullptr) { + throw runtime_error("Simulation has to be initialized before loading serial scan parameters"); } string nextToken, line; @@ -1306,7 +1257,7 @@ void FVSolver::loadSerialScanParameters(istream& ifsInput, int numSerialScanPara istringstream lineInput(line); lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "SERIAL_SCAN_PARAMETER_END") { @@ -1326,14 +1277,14 @@ PARAMETER_SCAN_BEGIN 2 2.0 PARAMETER_SCAN_END */ -void FVSolver::loadSerialScanParameterValues(istream& ifsInput, int numSerialScanParameterValues) { - if (simulation == 0) { - throw "Simulation has to be initialized before loading serial scan parameter values"; +void FVSolver::loadSerialScanParameterValues(SimTool* simTool, SimulationExpression* simulation, istream& ifsInput, int numSerialScanParameterValues) { + if (simulation == nullptr) { + throw runtime_error("Simulation has to be initialized before loading serial scan parameter values"); } string nextToken, line; int numSerialScanParameters = simulation->getNumParameters(); - double** serialScanParameterValues = new double*[numSerialScanParameterValues]; + auto** serialScanParameterValues = new double*[numSerialScanParameterValues]; for (int i = 0; i < numSerialScanParameterValues; i ++) { serialScanParameterValues[i] = new double[numSerialScanParameters]; @@ -1352,7 +1303,7 @@ void FVSolver::loadSerialScanParameterValues(istream& ifsInput, int numSerialSca lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "SERIAL_SCAN_PARAMETER_VALUE_END") { @@ -1448,10 +1399,9 @@ JUMP_CONDITION_END MEMBRANE_END */ -void FVSolver::createSimTool(istream& ifsInput, int taskID) +SimTool* FVSolver::createSimTool(istream& ifsInput, istream& vcgInput, int taskID, bool bSimZip) { - SimTool::create(); - simTool = SimTool::getInstance(); + SimTool* simTool = new SimTool(); if (taskID < 0) { // no messaging SimulationMessaging::create(); @@ -1478,66 +1428,83 @@ void FVSolver::createSimTool(istream& ifsInput, int taskID) SimulationMessaging::getInstVar()->setWorkerEvent(new WorkerEvent(JOB_STARTING, "preprocessing started")); } else if (nextToken == "SIMULATION_PARAM_BEGIN") { - loadSimulationParameters(ifsInput); + loadSimulationParameters(simTool, ifsInput); simTool->checkTaskIdLockFile();//check if taskid compatible with taskid lockfile } else if (nextToken == "POST_PROCESSING_BLOCK_BEGIN") { - loadPostProcessingBlock(ifsInput); + loadPostProcessingBlock(simTool, ifsInput); } else if (nextToken == "MODEL_BEGIN") { - loadModel(ifsInput); - if (model == NULL) { - throw "Model has 0 features"; + VCellModel* model = loadModel(ifsInput); + if (model == nullptr) { + throw runtime_error("Model has 0 features"); } simTool->setModel(model); } else if (nextToken == "MESH_BEGIN") { - loadMesh(ifsInput); + // consume lines until MESH_END + while (!ifsInput.eof()) + { + getline(ifsInput, line); + istringstream lineInput(line); + + nextToken = ""; + lineInput >> nextToken; + if (nextToken.size() == 0 || nextToken[0] == '#') { + continue; + } + if (nextToken == "MESH_END") { + break; + } + } + loadMeshFromVcg(simTool->getModel(), vcgInput); } else if (nextToken == "VARIABLE_BEGIN") { - loadSimulation(ifsInput); + SimulationExpression* simulation = loadSimulation(simTool, simTool->getModel(), ifsInput); simTool->setSimulation(simulation); } else if (nextToken == "SMOLDYN_BEGIN") { - loadSmoldyn(ifsInput); + loadSmoldyn(simTool, ifsInput); } else if (nextToken == "PARAMETER_BEGIN") { int numParams = 0; lineInput >> numParams; - loadParameters(ifsInput, numParams); + loadParameters(simTool->getSimulation(), ifsInput, numParams); } else if (nextToken == "SERIAL_SCAN_PARAMETER_BEGIN") { int numSerialScanParams = 0; lineInput >> numSerialScanParams; - loadSerialScanParameters(ifsInput, numSerialScanParams); + loadSerialScanParameters(simTool->getSimulation(), ifsInput, numSerialScanParams); } else if (nextToken == "SERIAL_SCAN_PARAMETER_VALUE_BEGIN") { int numSerialScanParamValues = 0; lineInput >> numSerialScanParamValues; - loadSerialScanParameterValues(ifsInput, numSerialScanParamValues); + loadSerialScanParameterValues(simTool, simTool->getSimulation(), ifsInput, numSerialScanParamValues); } else if (nextToken == "FIELD_DATA_BEGIN") { - loadFieldData(ifsInput); + loadFieldData(simTool->getSimulation(), ifsInput); } else if (nextToken == "COMPARTMENT_BEGIN") { string feature_name; lineInput >> feature_name; - Feature* feature = model->getFeatureFromName(feature_name); - if (feature != NULL) { - loadFeature(ifsInput, feature); + Feature* feature = simTool->getModel()->getFeatureFromName(feature_name); + if (feature != nullptr) { + loadFeature(simTool->getSimulation(), ifsInput, feature); } else { - throw "createSimTool(), Invalid compartment when loading feature!"; + throw runtime_error("createSimTool(), Invalid compartment when loading feature!"); } } else if (nextToken == "MEMBRANE_BEGIN") { string mem_name, feature1_name, feature2_name; lineInput >> mem_name >> feature1_name >> feature2_name; - Membrane* membrane = model->getMembraneFromName(mem_name); - if (membrane != 0) { - loadMembrane(ifsInput, membrane); + Membrane* membrane = simTool->getModel()->getMembraneFromName(mem_name); + if (membrane != nullptr) { + loadMembrane(simTool->getSimulation(), simTool->getModel(), ifsInput, membrane); } else { - throw "createSimTool(), Invalid compartment when loading membrane!"; + throw runtime_error("createSimTool(), Invalid compartment when loading membrane!"); } } else { - stringstream ss; - ss << "createSimTool(), encountered unknown token " << nextToken << endl; - throw ss.str(); + throw runtime_error("createSimTool(), encountered unknown token " + nextToken); } } + if (!bSimZip) { + simTool->requestNoZip(); + } + return simTool; } -void FVSolver::loadSmoldyn(istream& ifsInput) { - if (simulation == 0) { - throw "Simulation has to be initialized before loading smoldyn"; +void FVSolver::loadSmoldyn(SimTool* sim_tool, istream& ifsInput) { + if (sim_tool->getSimulation() == nullptr) { + throw runtime_error("Simulation has to be initialized before loading smoldyn"); } string nextToken, line; @@ -1547,7 +1514,7 @@ void FVSolver::loadSmoldyn(istream& ifsInput) { istringstream lineInput(line); lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "SMOLDYN_END") { @@ -1557,35 +1524,31 @@ void FVSolver::loadSmoldyn(istream& ifsInput) { string inputfile; getline(lineInput, inputfile); trimString(inputfile); - simTool->setSmoldynInputFile(inputfile); + sim_tool->setSmoldynInputFile(inputfile); } } } -FVSolver::FVSolver(istream& fvinput, int taskID, const char* outdir, bool bSimZip) { - simTool = 0; - simulation = 0; - model = 0; - mesh = 0; +FVSolver::FVSolver(const char* outdir) { + // simTool = nullptr; + // simulation = nullptr; + // model = nullptr; + mesh = nullptr; outputPath = outdir; - createSimTool(fvinput, taskID); - if (!bSimZip) { - SimTool::getInstance()->requestNoZip(); - } } FVSolver::~FVSolver() { delete SimulationMessaging::getInstVar(); - delete simulation; - delete model; + // delete simulation; + // delete model; delete mesh; - delete simTool; + // delete simTool; } -void FVSolver::solve(bool bLoadFinal, double* paramValues) +void FVSolver::solve(SimTool* simTool, bool bLoadFinal, double* paramValues) { - if (paramValues != 0) { - simulation->setParameterValues(paramValues); + if (paramValues != nullptr) { + simTool->getSimulation()->setParameterValues(simTool, paramValues); } simTool->setLoadFinal(bLoadFinal); @@ -1594,63 +1557,59 @@ void FVSolver::solve(bool bLoadFinal, double* paramValues) simTool->start(); } -void FVSolver::init(double* paramValues){ +void FVSolver::init(SimTool* sim_tool, double* paramValues){ // setting initial conditions. - simulation->setParameterValues(paramValues); - simulation->initSimulation(); + sim_tool->getSimulation()->setParameterValues(sim_tool, paramValues); + sim_tool->getSimulation()->initSimulation(sim_tool); } -void FVSolver::step(double* paramValues) +void FVSolver::step(SimTool* sim_tool, double* paramValues) { //simulation->setParameterValues(paramValues); - simulation->iterate(); - simulation->update(); + sim_tool->getSimulation()->iterate(sim_tool); + sim_tool->getSimulation()->update(sim_tool); } -string FVSolver::getVariableName(int index){ - return simulation->getVariable(index)->getName(); +string FVSolver::getVariableName(const SimTool* sim_tool, int index){ + return sim_tool->getSimulation()->getVariable(index)->getName(); } -int FVSolver::getNumVariables(){ - return simulation->getNumVariables(); +int FVSolver::getNumVariables(const SimTool* sim_tool){ + return sim_tool->getSimulation()->getNumVariables(); } -double* FVSolver::getValue(string& var, int arrayID) { +double* FVSolver::getValue(const SimTool* sim_tool, string& var, int arrayID) { if (arrayID==0){ - return simulation->getVariableFromName(var)->getOld(); + return sim_tool->getSimulation()->getVariableFromName(var)->getOld(); } else if (arrayID==1){ - return simulation->getVariableFromName(var)->getCurr(); + return sim_tool->getSimulation()->getVariableFromName(var)->getCurr(); } else { - throw "arrayID out of bounds"; + throw runtime_error("arrayID out of bounds"); } } -int FVSolver::getVariableLength(string& var) { - return simulation->getVariableFromName(var)->getSize(); +int FVSolver::getVariableLength(const SimTool* sim_tool, string& var) { + return sim_tool->getSimulation()->getVariableFromName(var)->getSize(); } -void FVSolver::setInitialCondition(string& varName, int dataLength, const double* data) { - Variable* var = simulation->getVariableFromName(varName); - if (var == 0) { - char errMsg[512]; - sprintf(errMsg, "FVSolver::setInitialCondition() : variable %s doesn't exist", varName.c_str()); - throw errMsg; +void FVSolver::setInitialCondition(SimTool* sim_tool, string& varName, int dataLength, const double* data) { + Variable* var = sim_tool->getSimulation()->getVariableFromName(varName); + if (var == nullptr) { + throw runtime_error("FVSolver::setInitialCondition() : variable " + varName + " doesn't exist"); } if (var->getSize() != dataLength) { - char errMsg[512]; - sprintf(errMsg, "FVSolver::setInitialCondition() : variable %s doesn't match in size, %d, %d", varName.c_str(), var->getSize(), dataLength); - throw errMsg; + throw runtime_error("FVSolver::setInitialCondition() : variable " + varName + " doesn't match in size, " + to_string(var->getSize()) + ", " + to_string(dataLength)); } memcpy(var->getCurr(), data, dataLength * sizeof(double)); var->update(); } -double FVSolver::getCurrentTime(){ - return simulation->getTime_sec(); +double FVSolver::getCurrentTime(SimTool* sim_tool){ + return sim_tool->getSimulation()->getTime_sec(sim_tool); } -void FVSolver::setEndTime(double endTime){ - simTool->setEndTimeSec(endTime); +void FVSolver::setEndTime(SimTool* sim_tool, double endTime){ + sim_tool->setEndTimeSec(endTime); } /** @@ -1659,14 +1618,14 @@ POST_PROCESSING_BLOCK_BEGIN PROJECTION_DATA_GENERATOR postDex cell x sum (5.0 * Dex_cell); POST_PROCESSING_BLOCK_END */ -void FVSolver::loadPostProcessingBlock(istream& ifsInput){ - if (simulation == 0) { - throw "Simulation has to be initialized before loading field data"; +void FVSolver::loadPostProcessingBlock(SimTool* sim_tool, istream& ifsInput){ + if (sim_tool->getSimulation() == nullptr) { + throw runtime_error("Simulation has to be initialized before loading field data"); } // create post processing block; - simulation->createPostProcessingBlock(); - PostProcessingBlock* postProcessingBlock = simulation->getPostProcessingBlock(); + sim_tool->getSimulation()->createPostProcessingBlock(); + PostProcessingBlock* postProcessingBlock = sim_tool->getSimulation()->getPostProcessingBlock(); // add var statistics data generator always postProcessingBlock->addDataGenerator(new VariableStatisticsDataGenerator()); @@ -1677,7 +1636,7 @@ void FVSolver::loadPostProcessingBlock(istream& ifsInput){ nextToken = ""; lineInput >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "POST_PROCESSING_BLOCK_END") { @@ -1687,15 +1646,15 @@ void FVSolver::loadPostProcessingBlock(istream& ifsInput){ if (nextToken == "PROJECTION_DATA_GENERATOR") { string name, domain_name, op, axis; lineInput >> name >> domain_name >> axis >> op; - Feature* feature = model->getFeatureFromName(domain_name); + Feature* feature = sim_tool->getModel()->getFeatureFromName(domain_name); Expression* function = readExpression(lineInput, name); - ProjectionDataGenerator* dataGenerator = new ProjectionDataGenerator(name, feature, axis, op, function); + auto* dataGenerator = new ProjectionDataGenerator(name, feature, axis, op, function); postProcessingBlock->addDataGenerator(dataGenerator); } else if (nextToken == "GAUSSIAN_CONVOLUTION_DATA_GENERATOR") { string name, domain_name; double sigmaXY, sigmaZ; lineInput >> name >> domain_name >> sigmaXY >> sigmaZ; - Feature* feature = model->getFeatureFromName(domain_name); + Feature* feature = sim_tool->getModel()->getFeatureFromName(domain_name); string volFunctionKeyword = "GAUSSIAN_CONVOLUTION_VOL_FUNCTION"; string memFunctionKeyword = "GAUSSIAN_CONVOLUTION_MEM_FUNCTION"; size_t indexOfVolumeFunctionKeyword = line.find(volFunctionKeyword); @@ -1704,21 +1663,21 @@ void FVSolver::loadPostProcessingBlock(istream& ifsInput){ size_t indexOfVolumeExp = indexOfVolumeFunctionKeyword+volFunctionKeyword.size(); string volumeFunctionString = line.substr(indexOfVolumeExp, indexOfMembraneFunctionKeyword-indexOfVolumeExp); trimString(volumeFunctionString); - Expression* volumeFunction = new Expression(volumeFunctionString); + auto* volumeFunction = new Expression(volumeFunctionString); size_t indexOfMembraneExp = indexOfMembraneFunctionKeyword+memFunctionKeyword.size(); string membraneFunctionString = line.substr(indexOfMembraneExp, string::npos); trimString(membraneFunctionString); - Expression* membraneFunction = new Expression(membraneFunctionString); - GaussianConvolutionDataGenerator* dataGenerator = new GaussianConvolutionDataGenerator(name, feature, sigmaXY, sigmaZ, volumeFunction, membraneFunction); + auto* membraneFunction = new Expression(membraneFunctionString); + auto* dataGenerator = new GaussianConvolutionDataGenerator(name, feature, sigmaXY, sigmaZ, volumeFunction, membraneFunction); postProcessingBlock->addDataGenerator(dataGenerator); string garbage; getline(lineInput, garbage); } else if (nextToken == "ROI_DATA_GENERATOR_BEGIN") { - int* volumePoints = 0; - int* membranePoints = 0; + int* volumePoints = nullptr; + int* membranePoints = nullptr; int numVolumePoints = 0; int numMembranePoints = 0; FieldData* sampleImage; @@ -1734,7 +1693,7 @@ void FVSolver::loadPostProcessingBlock(istream& ifsInput){ nextToken = ""; lineInput_roi >> nextToken; - if (nextToken.size() == 0 || nextToken[0] == '#') { + if (nextToken.empty() || nextToken[0] == '#') { continue; } if (nextToken == "ROI_DATA_GENERATOR_END") { @@ -1759,7 +1718,7 @@ void FVSolver::loadPostProcessingBlock(istream& ifsInput){ } else if (nextToken == "StoreEnabled") { lineInput_roi >> storeEnabledStr; if (storeEnabledStr == "false") { - simTool->setStoreEnable(false); + sim_tool->setStoreEnable(false); } } else if (nextToken == "SampleImage") { lineInput_roi >> numImageRegions >> zSlice; @@ -1776,12 +1735,10 @@ void FVSolver::loadPostProcessingBlock(istream& ifsInput){ sampleImage = new FieldData(0, VAR_VOLUME, "", fieldName, varName, time, fieldfilename); } }// end while loop for roi data generator - RoiDataGenerator* dataGenerator = new RoiDataGenerator(roiDGName, volumePoints, numVolumePoints, membranePoints, numMembranePoints, sampleImage, numImageRegions, zSlice); + auto* dataGenerator = new RoiDataGenerator(roiDGName, volumePoints, numVolumePoints, membranePoints, numMembranePoints, sampleImage, numImageRegions, zSlice); postProcessingBlock->addDataGenerator(dataGenerator); } else { - stringstream ss; - ss << "loadPostProcessingBlock(), encountered unknown token " << nextToken << endl; - throw ss.str(); + throw runtime_error("loadPostProcessingBlock(), encountered unknown token " + nextToken); } } } diff --git a/VCell/src/FastSystemExpression.cpp b/VCell/src/FastSystemExpression.cpp index b97e1b91a..72ba030fd 100644 --- a/VCell/src/FastSystemExpression.cpp +++ b/VCell/src/FastSystemExpression.cpp @@ -149,14 +149,14 @@ void FastSystemExpression::setCoordinates(double time_sec, WorldCoord& wc) { } -void FastSystemExpression::initVars() +void FastSystemExpression::initVars(SimTool* sim_tool) { // independent variables int numSymbols = 4 + dimension + numDependents; double* values = new double[numSymbols]; WorldCoord wc = simulation->getMesh()->getVolumeWorldCoord(currIndex); - values[0] = simulation->getTime_sec(); + values[0] = simulation->getTime_sec(sim_tool); values[1] = wc.x; values[2] = wc.y; values[3] = wc.z; @@ -188,7 +188,8 @@ void FastSystemExpression::updateDependentVars() { } } -void FastSystemExpression::setDependentVariables(string* vars) { +void FastSystemExpression::setDependentVariables(string* vars) const +{ for (int i = 0; i < numDependents; i ++) { pDependentVars[i] = simulation->getVariableFromName(vars[i]); if (pDependentVars[i] == NULL){ @@ -199,7 +200,8 @@ void FastSystemExpression::setDependentVariables(string* vars) { } } -void FastSystemExpression::setIndependentVariables(string* vars) { +void FastSystemExpression::setIndependentVariables(string* vars) const +{ for (int i = 0; i < dimension; i ++) { pVars[i] = simulation->getVariableFromName(vars[i]); if (pVars[i] == NULL){ @@ -211,10 +213,10 @@ void FastSystemExpression::setIndependentVariables(string* vars) { } SimpleSymbolTable* FastSystemExpression::getFastSymbolTable() { - if (fastSymbolTable == NULL) { + if (fastSymbolTable == nullptr) { for (int i = 0; i < dimension; i ++) { - if (pVars[i] == NULL) { - throw "No independent variables defined"; + if (pVars[i] == nullptr) { + throw std::runtime_error("No independent variables defined"); } } @@ -260,14 +262,15 @@ SimpleSymbolTable* FastSystemExpression::getFastSymbolTable() { return fastSymbolTable; } -void FastSystemExpression::updateIndepValues() { +void FastSystemExpression::updateIndepValues() const +{ int indepOffset = 4 + simulation->getNumFields() + simulation->getNumRandomVariables(); for (int i = 0; i < dimension; i ++) { fastValues[indepOffset + i] = getX(i); } } -void FastSystemExpression::resolveReferences(Simulation *sim) { +void FastSystemExpression::resolveReferences(SimulationExpression *sim) { bindAllExpressions(); } diff --git a/VCell/src/Feature.cpp b/VCell/src/Feature.cpp index 6120cebeb..94d937d87 100644 --- a/VCell/src/Feature.cpp +++ b/VCell/src/Feature.cpp @@ -28,16 +28,12 @@ Feature::Feature(string& name, unsigned char findex, FeatureHandle Ahandle) : St fastSystem = NULL; } -Feature::~Feature() -{ -} - //double Feature::getMaxIterationTime() //{ // return 0.0; //} -void Feature::resolveReferences(Simulation *sim) +void Feature::resolveReferences(SimulationExpression *sim) { for (int i = 0; i < (int)volumeVarContextList.size(); i ++) { VolumeVarContextExpression* volumeVarContext = volumeVarContextList[i]; @@ -120,12 +116,12 @@ VolumeRegionVarContextExpression* Feature::getVolumeRegionVarContext(VolumeRegio return 0; } -void Feature::reinitConstantValues() { +void Feature::reinitConstantValues(SimulationExpression* sim) { for (int i = 0; i < (int)volumeVarContextList.size(); i ++) { - volumeVarContextList[i]->reinitConstantValues(); + volumeVarContextList[i]->reinitConstantValues(sim); } for (int i = 0; i < (int)volumeRegionVarContextList.size(); i ++) { - volumeRegionVarContextList[i]->reinitConstantValues(); + volumeRegionVarContextList[i]->reinitConstantValues(sim); } } diff --git a/VCell/src/FiniteVolume.cpp b/VCell/src/FiniteVolume.cpp index 830a58d90..b8bfb28b7 100644 --- a/VCell/src/FiniteVolume.cpp +++ b/VCell/src/FiniteVolume.cpp @@ -35,7 +35,7 @@ void printUsage() { #ifdef USE_MESSAGING cout << "Arguments : [-d output] [-nz] [-tid taskID] fvInputFile" << endl; #else - cout << "Arguments : [-d output] [-nz] fvInputFile" << endl; + cout << "Arguments : [-d output] [-nz] fvInputFile vcgInputFile" << endl; #endif } @@ -49,15 +49,17 @@ int main(int argc, char *argv[]) int returnCode = 0; string errorMsg = "Exception : "; - char* outputPath = 0; - char* fvInputFile = 0; + char* outputPath = nullptr; + char* fvInputFile = nullptr; + char* vcgInputFile = nullptr; ifstream ifsInput; - FVSolver* fvSolver = NULL; + ifstream vcgInput; + FVSolver* fvSolver = nullptr; bool bSimZip = true; try { int taskID = -1; - if (argc < 2) { + if (argc < 3) { cout << "Missing arguments!" << endl; printUsage(); exit(1); @@ -94,8 +96,10 @@ int main(int argc, char *argv[]) printUsage(); exit(1); #endif - } else { + } else if (fvInputFile == nullptr){ fvInputFile = argv[i]; + } else { + vcgInputFile = argv[i]; } } struct stat buf; @@ -106,25 +110,39 @@ int main(int argc, char *argv[]) // strip " in case that file name has " around int fl = strlen(fvInputFile); - if (fvInputFile[0] == '"' && fvInputFile[fl-1] == '"') { - fvInputFile[fl-1] = 0; - fvInputFile ++; - } + if (fvInputFile[0] == '"' && fvInputFile[fl-1] == '"') { + fvInputFile[fl-1] = 0; + fvInputFile ++; + } ifsInput.open(fvInputFile); if (!ifsInput.is_open()) { cout << "File doesn't exist: " << fvInputFile << endl; exit(102); } - fvSolver = new FVSolver(ifsInput, taskID, outputPath, bSimZip); + // strip " in case that file name has " around + fl = strlen(vcgInputFile); + if (vcgInputFile[0] == '"' && vcgInputFile[fl-1] == '"') { + vcgInputFile[fl-1] = 0; + vcgInputFile ++; + } + vcgInput.open(vcgInputFile); + if (!vcgInput.is_open()) { + cout << "File doesn't exist: " << vcgInputFile << endl; + exit(102); + } + + fvSolver = new FVSolver(outputPath); + SimTool* sim_tool = fvSolver->createSimTool(ifsInput, vcgInput, taskID, bSimZip); ifsInput.close(); + vcgInput.close(); - if(fvSolver->getNumVariables() == 0){ + if(fvSolver->getNumVariables(sim_tool) == 0){ //sims with no reactions and no diffusing species cause exit logic to 'wait' forever //never sending a job failed or job finished message and never cleaning up threads throw invalid_argument("FiniteVolume error: Must have at least 1 variable or reaction to solve"); } - fvSolver->solve(); + fvSolver->solve(sim_tool); } catch (const char *exStr){ errorMsg += exStr; diff --git a/VCell/src/JumpCondition.cpp b/VCell/src/JumpCondition.cpp index 56d7abe35..6386664e0 100644 --- a/VCell/src/JumpCondition.cpp +++ b/VCell/src/JumpCondition.cpp @@ -64,7 +64,7 @@ double JumpCondition::evaluateExpression(double* values) { return expression->evaluateVector(values); } -bool JumpCondition::isConstantExpression() { +bool JumpCondition::isConstantExpression(SimulationExpression *sim) { // pure constant if (constantValue != 0) { return true; @@ -81,19 +81,19 @@ bool JumpCondition::isConstantExpression() { vector symbols; expression->getSymbols(symbols); for (int i = 0; i < (int)symbols.size(); i ++) { - if (!((SimulationExpression*)SimTool::getInstance()->getSimulation())->isParameter(symbols[i])) { + if (!sim->isParameter(symbols[i])) { return false; } } return true; } -void JumpCondition::reinitConstantValues() { - if (expression == 0 || !isConstantExpression()) { +void JumpCondition::reinitConstantValues(SimulationExpression *sim) { + if (expression == nullptr || !isConstantExpression(sim)) { return; } double d = expression->evaluateProxy(); - if (constantValue == 0) { + if (constantValue == nullptr) { constantValue = new double[1]; } constantValue[0] = d; diff --git a/VCell/src/Membrane.cpp b/VCell/src/Membrane.cpp index df4a781f8..1bd2da4c2 100644 --- a/VCell/src/Membrane.cpp +++ b/VCell/src/Membrane.cpp @@ -88,7 +88,7 @@ bool Membrane::inBetween(Feature* f1, Feature* f2) { return (f1 == feature1 && f2 == feature2 || f1 == feature2 && f2 == feature1); } -void Membrane::resolveReferences(Simulation *sim) +void Membrane::resolveReferences(SimulationExpression *sim) { for (int i = 0; i < (int)membraneVarContextList.size(); i ++) { MembraneVarContextExpression *membraneVarContext = membraneVarContextList[i]; @@ -105,11 +105,11 @@ void Membrane::resolveReferences(Simulation *sim) } } -void Membrane::reinitConstantValues() { +void Membrane::reinitConstantValues(SimulationExpression* sim) { for (int i = 0; i < (int)membraneVarContextList.size(); i ++) { - membraneVarContextList[i]->reinitConstantValues(); + membraneVarContextList[i]->reinitConstantValues(sim); } for (int i = 0; i < (int)membraneRegionVarContextList.size(); i ++) { - membraneRegionVarContextList[i]->reinitConstantValues(); + membraneRegionVarContextList[i]->reinitConstantValues(sim); } } diff --git a/VCell/src/MembraneEqnBuilderDiffusion.cpp b/VCell/src/MembraneEqnBuilderDiffusion.cpp index ec4c44e69..a877f1511 100644 --- a/VCell/src/MembraneEqnBuilderDiffusion.cpp +++ b/VCell/src/MembraneEqnBuilderDiffusion.cpp @@ -42,15 +42,15 @@ MembraneEqnBuilderDiffusion::~MembraneEqnBuilderDiffusion() { // Initializes Matrix only // //------------------------------------------------------------------ -void MembraneEqnBuilderDiffusion::initEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, +void MembraneEqnBuilderDiffusion::initEquation(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) { if (!bPreProcessed) { - preProcess(); + preProcess(model); } if (periodicPairs.size() > 0) { - initEquation_Periodic(deltaTime, volumeIndexStart, volumeIndexSize, membraneIndexStart, membraneIndexSize); + initEquation_Periodic(model, deltaTime, volumeIndexStart, volumeIndexSize, membraneIndexStart, membraneIndexSize); return; } @@ -113,11 +113,11 @@ void MembraneEqnBuilderDiffusion::initEquation(double deltaTime, int volumeIndex // updates B vector only with reaction rates and boundary conditions // //------------------------------------------------------------------ -void MembraneEqnBuilderDiffusion::buildEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, +void MembraneEqnBuilderDiffusion::buildEquation(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) { if (periodicPairs.size() > 0) { - buildEquation_Periodic(deltaTime, volumeIndexStart, volumeIndexSize, + buildEquation_Periodic(sim, deltaTime, volumeIndexStart, volumeIndexSize, membraneIndexStart, membraneIndexSize); return; } @@ -125,7 +125,6 @@ void MembraneEqnBuilderDiffusion::buildEquation(double deltaTime, int volumeInde VolumeElement *pVolumeElement = mesh->getVolumeElements(); MembraneElement *pMembraneElement = mesh->getMembraneElements(); - Simulation *sim = SimTool::getInstance()->getSimulation(); SparseMatrixPCG* membraneElementCoupling = ((CartesianMesh*)mesh)->getMembraneCoupling(); MembraneElement* membraneElement = pMembraneElement; @@ -260,10 +259,10 @@ double MembraneEqnBuilderDiffusion::computeDiffusionConstant(int meIndex, int ne return (Di + Dj < epsilon)?(0.0):(2 * Di * Dj/(Di + Dj)); } -void MembraneEqnBuilderDiffusion::initEquation_Periodic(double deltaTime, int volumeIndexStart, int volumeIndexSize, +void MembraneEqnBuilderDiffusion::initEquation_Periodic(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) { if (!bPreProcessed) { - preProcess(); + preProcess(model); } VolumeElement *pVolumeElement = mesh->getVolumeElements(); @@ -374,13 +373,12 @@ void MembraneEqnBuilderDiffusion::initEquation_Periodic(double deltaTime, int vo } // end i } -void MembraneEqnBuilderDiffusion::buildEquation_Periodic(double deltaTime, int volumeIndexStart, int volumeIndexSize, +void MembraneEqnBuilderDiffusion::buildEquation_Periodic(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) { VolumeElement *pVolumeElement = mesh->getVolumeElements(); MembraneElement *pMembraneElement = mesh->getMembraneElements(); ASSERTION(pMembraneElement); - Simulation *sim = SimTool::getInstance()->getSimulation(); SparseMatrixPCG* membraneElementCoupling = ((CartesianMesh*)mesh)->getMembraneCoupling(); MembraneElement* membraneElement = pMembraneElement; @@ -521,7 +519,7 @@ void MembraneEqnBuilderDiffusion::buildEquation_Periodic(double deltaTime, int v #define PERIODIC_GEOMETRY_ERROR_MESSAGE "Geometry is not compatible with periodic boundary conditions. Couldn't find corresponding membrane element." #define PERIODIC_NORMAL_ERROR_MESSAGE "Non periodic surface at periodic membrane (normals don't agree)."; -void MembraneEqnBuilderDiffusion::preProcess() { +void MembraneEqnBuilderDiffusion::preProcess(VCellModel* model) { if (bPreProcessed) { return; } @@ -531,7 +529,6 @@ void MembraneEqnBuilderDiffusion::preProcess() { // check if there is periodic boundary condition in the model bool bHasPeriodicBC = false; - VCellModel* model = SimTool::getInstance()->getModel(); for (int i = 0; i < model->getNumMembranes(); i ++) { Membrane *membrane = model->getMembraneFromIndex(i); if (membrane->getXmBoundaryType() == BOUNDARY_PERIODIC diff --git a/VCell/src/MembraneEqnBuilderForward.cpp b/VCell/src/MembraneEqnBuilderForward.cpp index 7631ba7e5..e16d613d6 100644 --- a/VCell/src/MembraneEqnBuilderForward.cpp +++ b/VCell/src/MembraneEqnBuilderForward.cpp @@ -16,10 +16,8 @@ MembraneEqnBuilderForward::MembraneEqnBuilderForward(MembraneVariable *Avar, Mes odeSolver = Asolver; } -void MembraneEqnBuilderForward::buildEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) +void MembraneEqnBuilderForward::buildEquation(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) { - Simulation *sim = SimTool::getInstance()->getSimulation(); - ASSERTION((membraneIndexStart>=0) && ((membraneIndexStart+membraneIndexSize)<=mesh->getNumMembraneElements())); MembraneElement *pMembraneElement = mesh->getMembraneElements() + membraneIndexStart; diff --git a/VCell/src/MembraneRegionEqnBuilder.cpp b/VCell/src/MembraneRegionEqnBuilder.cpp index 6c9829900..8b2adb5e9 100644 --- a/VCell/src/MembraneRegionEqnBuilder.cpp +++ b/VCell/src/MembraneRegionEqnBuilder.cpp @@ -18,7 +18,7 @@ MembraneRegionEqnBuilder::MembraneRegionEqnBuilder(MembraneRegionVariable *Avar, odeSolver = Asolver; } -void MembraneRegionEqnBuilder::buildEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) +void MembraneRegionEqnBuilder::buildEquation(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) { int size = ((CartesianMesh*)mesh)->getNumMembraneRegions(); ASSERTION(size==var->getSize()); diff --git a/VCell/src/MembraneRegionVarContextExpression.cpp b/VCell/src/MembraneRegionVarContextExpression.cpp index 70285dad7..d258bc94a 100644 --- a/VCell/src/MembraneRegionVarContextExpression.cpp +++ b/VCell/src/MembraneRegionVarContextExpression.cpp @@ -17,9 +17,9 @@ double MembraneRegionVarContextExpression::getMembraneReactionRate(MembraneEleme return evaluateExpression(element, REACT_RATE_EXP); } -void MembraneRegionVarContextExpression::resolveReferences(Simulation* sim) { +void MembraneRegionVarContextExpression::resolveReferences(SimulationExpression* sim) { VarContext::resolveReferences(sim); - bindAll((SimulationExpression*)sim); + bindAll(sim); } double MembraneRegionVarContextExpression::getInitialValue(long regionIndex) { @@ -30,7 +30,7 @@ double MembraneRegionVarContextExpression::getUniformRate(MembraneRegion *region return evaluateMembraneRegionExpression(region->getIndex(), UNIFORM_RATE_EXP); } -bool MembraneRegionVarContextExpression::isNullExpressionOK(int expIndex) { +bool MembraneRegionVarContextExpression::isNullExpressionOK(int expIndex) const { if (expIndex == INITIAL_VALUE_EXP || expIndex == REACT_RATE_EXP || expIndex == UNIFORM_RATE_EXP) { return false; } diff --git a/VCell/src/MembraneVarContextExpression.cpp b/VCell/src/MembraneVarContextExpression.cpp index 7e25b6bad..41a84df62 100644 --- a/VCell/src/MembraneVarContextExpression.cpp +++ b/VCell/src/MembraneVarContextExpression.cpp @@ -4,22 +4,19 @@ */ #include #include -#include -#include -#include #include #include #include -#include +#include MembraneVarContextExpression::MembraneVarContextExpression(Membrane *membrane, MembraneVariable* var) : VarContext(membrane, var) { } -void MembraneVarContextExpression::resolveReferences(Simulation* sim) { +void MembraneVarContextExpression::resolveReferences(SimulationExpression* sim) { VarContext::resolveReferences(sim); - bindAll((SimulationExpression*)sim); + bindAll(sim); } double MembraneVarContextExpression::getInitialValue(MembraneElement *element){ @@ -71,11 +68,12 @@ double MembraneVarContextExpression::getZpBoundaryFlux(MembraneElement *element) return evaluateExpression(element, BOUNDARY_ZP_EXP); } -bool MembraneVarContextExpression::isNullExpressionOK(int expIndex) { +bool MembraneVarContextExpression::isNullExpressionOK(int expIndex) const +{ if (expIndex == INITIAL_VALUE_EXP || expIndex == REACT_RATE_EXP) { return false; } - Solver* solver = SimTool::getInstance()->getSimulation()->getSolverFromVariable(species); + Solver* solver = sim->getSolverFromVariable(species); if (solver != 0 && solver->isPDESolver()) { if (expIndex == DIFF_RATE_EXP) { return false; diff --git a/VCell/src/ODESolver.cpp b/VCell/src/ODESolver.cpp index 060e5a742..8ee1c28f9 100644 --- a/VCell/src/ODESolver.cpp +++ b/VCell/src/ODESolver.cpp @@ -38,7 +38,7 @@ ODESolver::~ODESolver() if (rate) delete[] rate; } -void ODESolver::solveEqn(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime) +void ODESolver::solveEqn(SimTool* sim_tool, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime) { if(arraySize==0){ ASSERTION((volumeIndexStart>=0) && diff --git a/VCell/src/PDESolver.cpp b/VCell/src/PDESolver.cpp index 4fa8ddb80..33e5cd86c 100644 --- a/VCell/src/PDESolver.cpp +++ b/VCell/src/PDESolver.cpp @@ -11,16 +11,13 @@ PDESolver::PDESolver(Variable *Var, bool AbTimeDependent) : Solver(Var) bTimeDependent = AbTimeDependent; } -PDESolver::~PDESolver() -{ -} -void PDESolver::initEqn(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime) { +void PDESolver::initEqn(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime) { if (!bFirstTime && !isTimeDependent()) { return; } ASSERTION(eqnBuilder); - eqnBuilder->initEquation(deltaTime, volumeIndexStart, volumeIndexSize, membraneIndexStart, membraneIndexSize); + eqnBuilder->initEquation(model, deltaTime, volumeIndexStart, volumeIndexSize, membraneIndexStart, membraneIndexSize); } diff --git a/VCell/src/PostProcessingHdf5Writer.cpp b/VCell/src/PostProcessingHdf5Writer.cpp index 2167ea131..56dbca81b 100644 --- a/VCell/src/PostProcessingHdf5Writer.cpp +++ b/VCell/src/PostProcessingHdf5Writer.cpp @@ -177,7 +177,7 @@ void PostProcessingHdf5Writer::createGroups() { } } -void PostProcessingHdf5Writer::writeOutput() { +void PostProcessingHdf5Writer::writeOutput(SimTool* sim_tool) { try { int timesRank = 1; hsize_t timesDims = 1; @@ -187,7 +187,7 @@ void PostProcessingHdf5Writer::writeOutput() { createGroups(); // write current time - double currTime = postProcessingBlock->simulation->getTime_sec(); + double currTime = postProcessingBlock->simulation->getTime_sec(sim_tool); hsize_t size = timeList.size() + 1; timesDataSet->extend(&size); hsize_t dim = 1; diff --git a/VCell/src/Scheduler.cpp b/VCell/src/Scheduler.cpp index 9b83d7f20..778b195cc 100644 --- a/VCell/src/Scheduler.cpp +++ b/VCell/src/Scheduler.cpp @@ -22,7 +22,7 @@ Scheduler::Scheduler(Simulation *Asim) bHasFastSystem = false; } -void Scheduler::initValues() +void Scheduler::initValues(VCellModel* model) { bFirstTime = true; // @@ -80,10 +80,10 @@ void Scheduler::initValues() } } - bHasFastSystem = SimTool::getInstance()->getModel()->hasFastSystem(); + bHasFastSystem = model->hasFastSystem(); } -void Scheduler::solveFastSystem(int volStart, int volSize, int memStart, int memSize) +void Scheduler::solveFastSystem(SimTool* sim_tool, int volStart, int volSize, int memStart, int memSize) { Feature *feature = NULL; FastSystem *fs = NULL; @@ -104,8 +104,8 @@ void Scheduler::solveFastSystem(int volStart, int volSize, int memStart, int mem if (fs = feature->getFastSystem()){ fs->setCurrIndex(i); WorldCoord wc = mesh->getVolumeWorldCoord(i); - fs->setCoordinates(sim->getTime_sec(), wc); - fs->initVars(); + fs->setCoordinates(sim->getTime_sec(sim_tool), wc); + fs->initVars(sim_tool); fs->solveSystem(); fs->updateVars(); } @@ -121,8 +121,8 @@ void Scheduler::solveFastSystem(int volStart, int volSize, int memStart, int mem if (fs = membrane->getFastSystem()){ fs->setCurrIndex(i); WorldCoord wc = mesh->getMembraneWorldCoord(i); - fs->setCoordinates(sim->getTime_sec(), wc); - fs->initVars(); + fs->setCoordinates(sim->getTime_sec(sim_tool), wc); + fs->initVars(sim_tool); fs->solveSystem(); fs->updateVars(); } diff --git a/VCell/src/SerialScheduler.cpp b/VCell/src/SerialScheduler.cpp index 24bf27b9b..ac661ddbd 100644 --- a/VCell/src/SerialScheduler.cpp +++ b/VCell/src/SerialScheduler.cpp @@ -6,7 +6,6 @@ #include #include #include -#include #include #include @@ -15,10 +14,9 @@ SerialScheduler::SerialScheduler(Simulation *Asim) { } -void SerialScheduler::iterate() +void SerialScheduler::iterate(SimTool* sim_tool) { - - VCellModel *model = SimTool::getInstance()->getModel(); + /* Contour *contour = NULL; int numContours = model->getNumContours(); @@ -71,25 +69,25 @@ void SerialScheduler::iterate() for (int i = 0; i < sim->getNumSolvers(); i ++) { solver = sim->getSolver(i); string timername = solver->getVar()->getName() + " Build"; - TimerHandle tHndBuild = SimTool::getInstance()->getTimerHandle(timername); + TimerHandle tHndBuild = sim_tool->getTimerHandle(timername); timername = solver->getVar()->getName() + " Solve"; - TimerHandle tHndSolve = SimTool::getInstance()->getTimerHandle(timername); + TimerHandle tHndSolve = sim_tool->getTimerHandle(timername); // // initialize equations first time around // - solver->initEqn(sim->getDT_sec(),0,volumeSize,0,membraneSize, bFirstTime); + solver->initEqn(sim_tool->getModel(), sim->getDT_sec(),0,volumeSize,0,membraneSize, bFirstTime); - SimTool::getInstance()->startTimer(tHndBuild); - solver->buildEqn(sim->getDT_sec(),0,volumeSize,0,membraneSize, bFirstTime); - SimTool::getInstance()->stopTimer(tHndBuild); + sim_tool->startTimer(tHndBuild); + solver->buildEqn(sim, sim->getDT_sec(),0,volumeSize,0,membraneSize, bFirstTime); + sim_tool->stopTimer(tHndBuild); - SimTool::getInstance()->startTimer(tHndSolve); - solver->solveEqn(sim->getDT_sec(),0,volumeSize,0,membraneSize, bFirstTime); - SimTool::getInstance()->stopTimer(tHndSolve); + sim_tool->startTimer(tHndSolve); + solver->solveEqn(sim_tool, sim->getDT_sec(),0,volumeSize,0,membraneSize, bFirstTime); + sim_tool->stopTimer(tHndSolve); } - if(hasFastSystem()){ + if(hasFastSystem(sim_tool->getModel())){ Mesh *mesh = sim->getMesh(); - solveFastSystem(0, mesh->getNumVolumeElements(), 0, mesh->getNumMembraneElements()); + solveFastSystem(sim_tool, 0, mesh->getNumVolumeElements(), 0, mesh->getNumMembraneElements()); } bFirstTime=false; } diff --git a/VCell/src/SimTool.cpp b/VCell/src/SimTool.cpp index 069951644..acfd81c6c 100644 --- a/VCell/src/SimTool.cpp +++ b/VCell/src/SimTool.cpp @@ -16,13 +16,11 @@ using std::endl; #include #include #include -#include #include #include #include #include #include -#include #include #include #include @@ -75,17 +73,15 @@ using std::max; #endif */ -SimTool* SimTool::instance = 0; - static int NUM_TOKENS_PER_LINE = 4; -static const int numRetries = 2; -static const int retryWaitSeconds = 5; +static constexpr int numRetries = 2; +static constexpr int retryWaitSeconds = 5; SimTool::SimTool() :bSimZip(true ), - vcellModel(0), - simulation(0), - _timer(0), + vcellModel(nullptr), + simulation(nullptr), + _timer(nullptr), bSimFileCompress(false), simEndTime(0), simStartTime(0), @@ -94,13 +90,10 @@ SimTool::SimTool() smoldynStepMultiplier(1), keepEvery(100), bStoreEnable(true), - baseFileName(0), simFileCount(0), - baseSimName(0), - baseDirName(0 ), zipFileCount(0), solver(FV_SOLVER ), - discontinuityTimes(0), + discontinuityTimes(nullptr), numDiscontinuityTimes(0), bLoadFinal(true ), sundialsSolverOptions( ), @@ -111,22 +104,15 @@ SimTool::SimTool() bSundialsOneStepOutput(false), keepAtMost(5000), - serialScanParameterValues(0), + serialScanParameterValues(nullptr), numSerialParameterScans(0), - postProcessingHdf5Writer(0), - smoldynSim(0 ), - smoldynInputFile("" ) + postProcessingHdf5Writer(nullptr), + smoldynSim(nullptr) { } SimTool::~SimTool() { - if (baseSimName != baseDirName) { - delete[] baseSimName; - } - delete[] baseDirName; - delete[] baseFileName; - delete[] discontinuityTimes; for (int i = 0; i < numSerialParameterScans; i ++) { @@ -138,15 +124,15 @@ SimTool::~SimTool() } void SimTool::setModel(VCellModel* model) { - if (model == 0) { - throw "SimTool::setModel(), model can't be null"; + if (model == nullptr) { + throw runtime_error("SimTool::setModel(), model can't be null"); } vcellModel = model; } -void SimTool::setSimulation(Simulation* sim) { - if (sim == 0) { - throw "SimTool::setSimulation(), simulation can't be null"; +void SimTool::setSimulation(SimulationExpression* sim) { + if (sim == nullptr) { + throw runtime_error("SimTool::setSimulation(), simulation can't be null"); } simulation = sim; simulation->setDT_sec(simDeltaTime); @@ -158,26 +144,11 @@ void SimTool::setTimeStep(double period) { void SimTool::setSmoldynStepMultiplier(int steps) { if(steps <= 0) { - throw "SimTool::setSmoldynStepMultiplier(), smoldynStepMultiplier must be a positive integer"; + throw runtime_error("SimTool::setSmoldynStepMultiplier(), smoldynStepMultiplier must be a positive integer"); } smoldynStepMultiplier = steps; } -void SimTool::create() { - if (instance == 0) { - instance = new SimTool(); - } else { - throw "SimTool (singleton) has been created"; - } -} - -SimTool* SimTool::getInstance() { - if (instance == 0) { - throw "SimTool (singleton) has not been created"; - } - return instance; -} - void SimTool::requestNoZip() { bSimZip = false; } @@ -188,23 +159,23 @@ bool SimTool::checkStopRequested() { TimerHandle SimTool::getTimerHandle(string& identifier) { - if (_timer==NULL){ + if (_timer==nullptr){ _timer = new Timer(); } return _timer->registerID(identifier); } -void SimTool::startTimer(TimerHandle hnd) +void SimTool::startTimer(TimerHandle hnd) const { _timer->start(hnd); } -void SimTool::stopTimer(TimerHandle hnd) +void SimTool::stopTimer(TimerHandle hnd) const { _timer->stop(hnd); } -double SimTool::getElapsedTimeSec(TimerHandle hnd) +double SimTool::getElapsedTimeSec(TimerHandle hnd) const { double time=0.0; _timer->getElapsedTimeSec(hnd, time); @@ -215,7 +186,7 @@ void SimTool::showSummary(FILE *fp) { fprintf(fp,"--------------- sim summary ----------------\n"); simulation->getMesh()->showSummary(fp); - fprintf(fp,"\ttime = %lg sec\n",simulation->getTime_sec()); + fprintf(fp,"\ttime = %lg sec\n",simulation->getTime_sec(this)); fprintf(fp,"\tdT = %lg sec\n",simulation->getDT_sec()); fprintf(fp,"--------------------------------------------\n"); if (_timer){ @@ -225,28 +196,13 @@ void SimTool::showSummary(FILE *fp) } } -void SimTool::setBaseFilename(char *fname) { - if (fname == 0 || strlen(fname) == 0) { - throw "invalid base file name for data set"; - } - baseFileName = new char[strlen(fname) + 1]; - memset(baseFileName, 0, strlen(fname) + 1); - memcpy(baseFileName, fname, strlen(fname) * sizeof(char)); - - // extract directory - baseDirName = new char[strlen(baseFileName) + 1]; - baseSimName = NULL; - - strcpy(baseDirName, baseFileName); - char* p = strrchr(baseDirName, DIRECTORY_SEPARATOR); - if (p == NULL) { - baseSimName = baseDirName; - baseDirName = 0; - } else { - baseSimName = new char[strlen(p+1) + 1]; - strcpy(baseSimName, p + 1); - *(p + 1)= 0; +void SimTool::setBaseFilename(const std::filesystem::path& fname) { + if (fname.string().empty()) { + throw runtime_error("invalid base file name for data set"); } + baseFileName = fname; + baseDirName = baseFileName.parent_path(); + baseSimName = baseFileName.filename(); } static void retryWait(int seconds) { @@ -258,11 +214,11 @@ static void retryWait(int seconds) { } static FILE* openFileWithRetry(const char* fileName, const char* mode) { - FILE *fp = NULL; + FILE *fp = nullptr; for (int retry = 0; retry < numRetries; retry ++) { fp = fopen(fileName, mode); - if (fp != NULL) { + if (fp != nullptr) { break; } if (retry < numRetries - 1) { @@ -273,7 +229,7 @@ static FILE* openFileWithRetry(const char* fileName, const char* mode) { return fp; } -static bool zipUnzipWithRetry(bool bZip, const char* zipFileName, const char* simFileName, std::string* errmsg) { +static bool zipUnzipWithRetry(bool bZip, const filesystem::path& zipFileName, const filesystem::path& simFileName, std::string* errmsg) { bool bSuccess = false; for (int attemptNum = 0; attemptNum <= numRetries; attemptNum++) { try { @@ -301,7 +257,7 @@ static bool zipUnzipWithRetry(bool bZip, const char* zipFileName, const char* si } if (!bSuccess) { errmsg->clear(); - errmsg->append("Writing zip file <").append(zipFileName).append("> failed."); + errmsg->append("Writing zip file <").append(zipFileName.string()).append("> failed."); } return bSuccess; } @@ -313,7 +269,7 @@ void SimTool::loadFinal() return; } - if (smoldynSim != NULL) { + if (smoldynSim != nullptr) { clearLog(); return; } @@ -324,21 +280,20 @@ void SimTool::loadFinal() bool bStartOver = true; - std::string logFileName{baseFileName}; + std::filesystem::path logFileName = baseFileName; logFileName.append(LOG_FILE_EXT); - std::string zipFileName{baseFileName}; - const std::string dataFileName; + std::filesystem::path zipFileName = baseFileName; + const std::filesystem::path dataFileName; FILE* tidFP = lockForReadWrite(); - FILE* logFP = fopen(logFileName.c_str(), "r"); + FILE* logFP = fopen(logFileName.string().c_str(), "r"); - if (logFP != NULL) { + if (logFP != nullptr) { bStartOver = false; // log file exists, there is old data - struct stat buf; zipFileName.append("00").append(ZIP_FILE_EXT); - if (stat(zipFileName.c_str(), &buf)) { + if (exists(zipFileName)) { bSimZip = false; NUM_TOKENS_PER_LINE = 3; } else { @@ -347,7 +302,7 @@ void SimTool::loadFinal() } simStartTime = -1; - int tempIteration = -1, tempFileCount = 0, tempZipCount = 0; + int tempIteration = -1, tempFileCount = 0; int numTimes = 0; while (!feof(logFP)){ @@ -359,18 +314,18 @@ void SimTool::loadFinal() // // look for line with a valid filename (includes basename) // - if (strstr(logBuffer, baseSimName)){ + if (strstr(logBuffer, baseSimName.string().c_str())){ // // parse iteration number and time // - int numTokens = 0; + int numTokens; if (bSimZip) { - numTokens = sscanf(logBuffer, "%d %s %s %lg", &tempIteration, dataFileName.c_str(), zipFileName.c_str(), &simStartTime); + numTokens = sscanf(logBuffer, "%d %s %s %lg", &tempIteration, dataFileName.string().c_str(), zipFileName.string().c_str(), &simStartTime); } else { - numTokens = sscanf(logBuffer, "%d %s %lg", &tempIteration, dataFileName.c_str(), &simStartTime); + numTokens = sscanf(logBuffer, "%d %s %lg", &tempIteration, dataFileName.string().c_str(), &simStartTime); } if (numTokens != NUM_TOKENS_PER_LINE){ - printf("SimTool::load(), error reading log file %s, reading iteration\n", logFileName.c_str()); + printf("SimTool::load(), error reading log file %s, reading iteration\n", logFileName.string().c_str()); printf("error in line %d = '%s'\n", tempFileCount, logBuffer); bStartOver = true; break; @@ -389,55 +344,54 @@ void SimTool::loadFinal() if (!bStartOver) { if (bSimZip) { // check if zip file exists - std::string zipFileAbsoluteName; - const bool needs_path_completion = - strchr(zipFileName.c_str(), DIRECTORY_SEPARATOR) == 0 && baseDirName != 0; + std::filesystem::path zipFileAbsoluteName; + const bool needs_path_completion = zipFileName.is_absolute() && is_directory(baseDirName); if (needs_path_completion) { - zipFileAbsoluteName = std::string(baseDirName) + zipFileName; + zipFileAbsoluteName = baseDirName / zipFileName; } else { - zipFileAbsoluteName = std::string(zipFileName); + zipFileAbsoluteName = zipFileName; } - if (stat(zipFileAbsoluteName.c_str(), &buf)) { + if (exists(zipFileAbsoluteName)) { cout << "SimTool::loadFinal(), unable to open zip file <" << zipFileAbsoluteName << ">" << endl; bStartOver = true; } else { // unzip the file (without directory) into exdir, currently we // unzip the file to the current working directory std::string errmsg; - zipUnzipWithRetry(false, zipFileAbsoluteName.c_str(), dataFileName.c_str(), &errmsg); + zipUnzipWithRetry(false, zipFileAbsoluteName.string().c_str(), dataFileName.string().c_str(), &errmsg); } } if (!bStartOver) { // otherwise check if sim file exists - if (stat(dataFileName.c_str(), &buf)) { + if (exists(dataFileName)) { cout << "SimTool::loadFinal(), unable to open sim file <" << dataFileName << ">" << endl; bStartOver = true; } else { try { - if (postProcessingHdf5Writer != NULL) { + if (postProcessingHdf5Writer != nullptr) { bStartOver = !postProcessingHdf5Writer->loadFinal(numTimes); } if (!bStartOver) { - FVDataSet::read(dataFileName.c_str(), simulation); + FVDataSet::read(dataFileName.string().c_str(), simulation); simulation->setCurrIteration(tempIteration); // set start time on sundials if (isSundialsPdeSolver()) { - simulation->setSimStartTime(simStartTime); + simulation->setSimStartTime(this, simStartTime); } simFileCount = tempFileCount; if (bSimZip) { - remove(dataFileName.c_str()); - zipFileCount = getZipCount(&zipFileName); + remove(dataFileName.string().c_str()); + zipFileCount = getZipCount(zipFileName); // wrong zip file Name if (zipFileCount < 0) { // should never happen bStartOver = true; } else { // check if this zip file is already big enough - if (stat(zipFileName.c_str(), &buf) == 0) { - if (buf.st_size > ZIP_FILE_LIMIT) { + if (exists(zipFileName)) { + if (file_size(zipFileName) > ZIP_FILE_LIMIT) { zipFileCount ++; } } @@ -457,41 +411,41 @@ void SimTool::loadFinal() } // if (logFP != NULL) // close tid file - if (tidFP != 0) { + if (tidFP != nullptr) { fclose(tidFP); } if (bStartOver) { clearLog(); } } -void SimTool::checkTaskIdLockFile(){ +void SimTool::checkTaskIdLockFile() const +{ FILE* fp = lockForReadWrite(); - if (fp != 0 ) { + if (fp != nullptr) { fclose(fp); } } -FILE* SimTool::lockForReadWrite() { +FILE* SimTool::lockForReadWrite() const +{ const int myTaskID = SimulationMessaging::getInstVar()->getTaskID(); - if (myTaskID < 0) return 0; + if (myTaskID < 0) return nullptr; - std::string tidFileName{baseFileName}; + std::string tidFileName = baseFileName.string(); tidFileName.append(TID_FILE_EXT); bool bExist = false; - struct stat buf; + struct stat buf{}; if (stat(tidFileName.c_str(), &buf) == 0) { // if exists bExist = true; } FILE* fp = openFileWithRetry(tidFileName.c_str(), bExist ? "r+" : "w+"); - if (fp == 0){ - std::string errMsg; - errMsg.append("SimTool::lockForReadWrite() - error opening .tid file <").append(tidFileName).append(">"); - throw errMsg.c_str(); + if (fp == nullptr){ + throw runtime_error("SimTool::lockForReadWrite() - error opening .tid file <" + tidFileName + ">"); } if (bExist) { @@ -530,20 +484,20 @@ void SimTool::updateLog(double progress, double time, int iteration) std::string errorMsg; FILE* tidFP = lockForReadWrite(); - struct stat buf; + struct stat buf{}; #if ( defined(WIN32) || defined(WIN64) ) // Windows wstring TempPath; wchar_t wcharPath[128]; // This needs to be fixed when on windows. if (GetTempPathW(128, wcharPath)){ TempPath = wcharPath; }else{ - throw "failed to obtain system temp directory"; + throw runtime_error("failed to obtain system temp directory"); } std::string tempDir( TempPath.begin(), TempPath.end() ); #else std::string tempDir = "/tmp/"; char const *folder = getenv("TMPDIR"); - if (folder != 0) + if (folder != nullptr) tempDir = folder; if (tempDir.at(tempDir.length()-1) != '/') tempDir += '/'; @@ -566,9 +520,9 @@ void SimTool::updateLog(double progress, double time, int iteration) std::stringstream simFileNameStream; // write sim files to local if (bSimZip && bUseTempDir) { - simFileNameStream << tempDir << getFileName(baseSimName) << std::setfill('0') << std::setw(4) << simFileCount << SIM_FILE_EXT; + simFileNameStream << tempDir << getFileName(baseSimName.string()) << std::setfill('0') << std::setw(4) << simFileCount << SIM_FILE_EXT; } else { - simFileNameStream << baseFileName << std::setfill('0') << std::setw(4) << simFileCount << SIM_FILE_EXT; + simFileNameStream << baseFileName.string() << std::setfill('0') << std::setw(4) << simFileCount << SIM_FILE_EXT; } simFileName = simFileNameStream.str(); std::cout << "sim file name is " << simFileName << std::endl; @@ -578,12 +532,12 @@ void SimTool::updateLog(double progress, double time, int iteration) simulation->writeData(simFileName.c_str(), bSimFileCompress); - std::string logFileName{baseFileName}; + std::string logFileName = baseFileName.string(); logFileName.append(LOG_FILE_EXT); logFP = openFileWithRetry(logFileName.c_str(), "a"); - if (logFP == 0) { + if (logFP == nullptr) { errorMsg.append("SimTool::updateLog() - error opening log file <").append(logFileName).append(">"); bSuccess = false; } else { @@ -592,7 +546,7 @@ void SimTool::updateLog(double progress, double time, int iteration) if (bSimZip) { //int retcode = 0; std::stringstream zipFileNameStream; - zipFileNameStream << baseFileName << std::setfill('0') << std::setw(2) << zipFileCount << ZIP_FILE_EXT; + zipFileNameStream << baseFileName.string() << std::setfill('0') << std::setw(2) << zipFileCount << ZIP_FILE_EXT; std::string zipFileName{zipFileNameStream.str()}; if (stat(particleFileName.c_str(), &buf) == 0) { // has particle // retcode = zip32(2, zipFileName, simFileName, particleFileName); @@ -606,17 +560,17 @@ void SimTool::updateLog(double progress, double time, int iteration) // write the log file if (bSuccess) { // write hdf5 post processing before writing log entry - if (postProcessingHdf5Writer != NULL) { - postProcessingHdf5Writer->writeOutput(); + if (postProcessingHdf5Writer != nullptr) { + postProcessingHdf5Writer->writeOutput(this); } std::stringstream zipFileNameWithoutPathStream; - zipFileNameWithoutPathStream << baseSimName << std::setfill('0') << std::setw(2) << + zipFileNameWithoutPathStream << baseSimName.string() << std::setfill('0') << std::setw(2) << zipFileCount << ZIP_FILE_EXT; std::stringstream simFileNameWithoutPathStream; - simFileNameWithoutPathStream << baseSimName << std::setfill('0') << std::setw(4) << + simFileNameWithoutPathStream << baseSimName.string() << std::setfill('0') << std::setw(4) << simFileCount << SIM_FILE_EXT; fprintf(logFP,"%4d %s %s %.15lg\n", iteration, simFileNameWithoutPathStream.str().c_str(), zipFileNameWithoutPathStream.str().c_str(), time); @@ -629,7 +583,7 @@ void SimTool::updateLog(double progress, double time, int iteration) } } else { // old format, no zip std::stringstream simFileNameWithoutPathStream; - simFileNameWithoutPathStream << baseSimName << std::setfill('0') << std::setw(4) << + simFileNameWithoutPathStream << baseSimName.string() << std::setfill('0') << std::setw(4) << simFileCount << SIM_FILE_EXT; fprintf(logFP,"%4d %s %.15lg\n", iteration, simFileNameWithoutPathStream.str().c_str(), time); @@ -638,18 +592,18 @@ void SimTool::updateLog(double progress, double time, int iteration) // close log file fclose(logFP); // close tid file - if (tidFP != 0) { + if (tidFP != nullptr) { fclose(tidFP); } if (bSuccess) { simFileCount++; } else { - throw errorMsg.c_str(); + throw runtime_error(errorMsg); } } else{ // write hdf5 post processing before writing log entry - if (postProcessingHdf5Writer != NULL) { - postProcessingHdf5Writer->writeOutput(); + if (postProcessingHdf5Writer != nullptr) { + postProcessingHdf5Writer->writeOutput(this); } } @@ -657,9 +611,9 @@ void SimTool::updateLog(double progress, double time, int iteration) SimulationMessaging::getInstVar()->setWorkerEvent(new WorkerEvent(JOB_DATA, progress, time)); } -int SimTool::getZipCount(char* zipFileName) { - char* p = strstr(zipFileName, ZIP_FILE_EXT); - if (p == NULL) { +int SimTool::getZipCount(const filesystem::path& zipFileName) { + const char* p = strstr(zipFileName.string().c_str(), ZIP_FILE_EXT); + if (p == nullptr) { return -1; } @@ -669,21 +623,14 @@ int SimTool::getZipCount(char* zipFileName) { return atoi(str); } -int SimTool::getZipCount(const std::string* zipFileName) { - // We need a char buffer because underlying ststr() call needs non-const char-ptr - char buffer[zipFileName->size()]; // C99 is great - strcpy(buffer, zipFileName->c_str()); - return this->getZipCount(buffer); -} - void SimTool::clearLog() { simStartTime = 0; simFileCount = 0; zipFileCount = 0; - if (SimTool::getInstance()->isSundialsPdeSolver()) { - simulation->setSimStartTime(0); + if (isSundialsPdeSolver()) { + simulation->setSimStartTime(this, 0); } FILE *fp; @@ -691,25 +638,25 @@ void SimTool::clearLog() std::string buffer; // remove mesh file - buffer.append(baseFileName).append(MESH_FILE_EXT); + buffer.append(baseFileName.string()).append(MESH_FILE_EXT); remove(buffer.c_str()); buffer.clear(); - buffer.append(baseFileName).append(MESHMETRICS_FILE_EXT); + buffer.append(baseFileName.string()).append(MESHMETRICS_FILE_EXT); remove(buffer.c_str()); buffer.clear(); - buffer.append(baseFileName).append(ZIP_FILE_EXT); + buffer.append(baseFileName.string()).append(ZIP_FILE_EXT); remove(buffer.c_str()); buffer.clear(); - buffer.append(baseFileName).append("00").append(ZIP_FILE_EXT); + buffer.append(baseFileName.string()).append("00").append(ZIP_FILE_EXT); remove(buffer.c_str()); buffer.clear(); - logFileName.append(baseFileName).append(LOG_FILE_EXT); + logFileName.append(baseFileName.string()).append(LOG_FILE_EXT); - if ((fp=fopen(logFileName.c_str(), "r")) == NULL){ + if ((fp=fopen(logFileName.c_str(), "r")) == nullptr){ printf("error opening log file <%s>\n", logFileName.c_str()); return; } @@ -720,7 +667,7 @@ void SimTool::clearLog() double time; while (true) { - int numTokens = 0; + int numTokens; if (bSimZip) { numTokens = fscanf(fp,"%d %s %s %lg\n", &iteration, simFileName.c_str(), zipFileName.c_str(), &time); } else { @@ -737,16 +684,16 @@ void SimTool::clearLog() *dotSim = '\0'; - buffer.append(baseFileName).append(SIM_FILE_EXT); + buffer.append(baseFileName.string()).append(SIM_FILE_EXT); remove(buffer.c_str()); buffer.clear(); - buffer.append(baseFileName).append(SIM_FILE_EXT).append(PARTICLE_FILE_EXT); + buffer.append(baseFileName.string()).append(SIM_FILE_EXT).append(PARTICLE_FILE_EXT); remove(buffer.c_str()); buffer.clear(); if (!bSimZip) continue; - const int count = getZipCount(&zipFileName); + const int count = getZipCount(zipFileName); if (oldCount != count && count >= 0) { remove(zipFileName.c_str()); oldCount = count; @@ -758,15 +705,16 @@ void SimTool::clearLog() remove(logFileName.c_str()); } -bool SimTool::isSundialsPdeSolver() { +bool SimTool::isSundialsPdeSolver() const +{ return solver == SUNDIALS_PDE_SOLVER; } -void SimTool::setSolver(string& s) { - if (s.length() > 0 && s != FV_SOLVER && s != SUNDIALS_PDE_SOLVER) { +void SimTool::setSolver(const string& s) { + if (!s.empty() && s != FV_SOLVER && s != SUNDIALS_PDE_SOLVER) { stringstream ss; ss << "unknown solver : " << s; - throw ss.str(); + throw runtime_error(ss.str()); } solver = s; } @@ -775,7 +723,7 @@ void SimTool::start() { if (simulation->getNumVariables() == 0) { return; } - simulation->resolveReferences(); + simulation->resolveReferences(this); if (numSerialParameterScans == 0 ) { // Is Not paramScan, the .fvinput had no SERIAL_SCAN_PARAMETER_... section start1(); @@ -783,47 +731,47 @@ void SimTool::start() { // Is paramScan, the .fvinput had 'SERIAL_SCAN_PARAMETER_...' section // Is also a remote server (green-button) run , the .fvinput has 'JMS_PARAM_BEGIN' section // Each paramScan is run as a separate job with serial jobIndex in the 'JMS_PARAM_BEGIN' section - SimulationExpression* sim = (SimulationExpression*)simulation; - sim->setParameterValues(serialScanParameterValues[SimulationMessaging::getInstVar()->getJobIndex()]); + SimulationExpression* sim = simulation; + sim->setParameterValues(this, serialScanParameterValues[SimulationMessaging::getInstVar()->getJobIndex()]); start1(); } else { // Is paramScan, the .fvinput had SERIAL_SCAN_PARAMETER_... section // Is also a local (blue-button) run where each paramScan is run in the following loop - SimulationExpression* sim = (SimulationExpression*)simulation; + SimulationExpression* sim = simulation; for (int scan = 0; scan < numSerialParameterScans; scan ++) { if (scan > 0) { - string bfn(baseFileName); + string bfn(baseFileName.string()); char oldIndex[10], newIndex[10]; sprintf(oldIndex, "_%d_\0", scan - 1); sprintf(newIndex, "_%d_\0", scan); int p = (int)bfn.rfind(oldIndex); bfn.replace(p, strlen(oldIndex), newIndex); - setBaseFilename((char*)bfn.c_str()); + setBaseFilename(filesystem::path(bfn)); } - sim->setParameterValues(serialScanParameterValues[scan]); + sim->setParameterValues(this, serialScanParameterValues[scan]); start1(); } } } -void SimTool::setSmoldynInputFile(string& inputfile) { +void SimTool::setSmoldynInputFile(const string& inputfile) { smoldynInputFile = inputfile; } void SimTool::start1() { // create post processing hdf5 writer - if (simulation->getPostProcessingBlock() != NULL) { - if (postProcessingHdf5Writer == NULL) { + if (simulation->getPostProcessingBlock() != nullptr) { + if (postProcessingHdf5Writer == nullptr) { std::string h5PPFileName; - h5PPFileName.append(baseFileName).append(HDF5_FILE_EXT); + h5PPFileName.append(baseFileName.string()).append(HDF5_FILE_EXT); postProcessingHdf5Writer = new PostProcessingHdf5Writer(h5PPFileName, simulation->getPostProcessingBlock()); } } - simulation->initSimulation(); + simulation->initSimulation(this); - if (smoldynInputFile != "") { + if (!smoldynInputFile.empty()) { smoldynSim = vcellhybrid::smoldynInit(this, smoldynInputFile);//smoldynInit will write output therefore computeHistogram is called by VCellSmoldynOutput.write(). copyParticleCountsToConcentration(); // since smoldyn only initializes variable current value, @@ -840,12 +788,12 @@ void SimTool::start1() { return; } - if (bStoreEnable && (baseFileName == NULL || strlen(baseFileName) == 0)) { - throw "Invalid base file name for dataset"; + if (bStoreEnable && !baseFileName.has_filename()) { + throw runtime_error("Invalid base file name "+baseFileName.string()+"for dataset"); } std::string message; - message.append("simulation [").append(baseSimName).append("] started"); + message.append("simulation [").append(baseSimName.string()).append("] started"); SimulationMessaging::getInstVar()->setWorkerEvent(new WorkerEvent(JOB_STARTING, message.c_str())); // @@ -859,23 +807,23 @@ void SimTool::start1() { // simulation starts from scratch, needs to // write .mesh and .meshmetrics file if (bStoreEnable){ - FILE *fp = NULL; + FILE *fp; std::string filename; - filename.append(baseFileName).append(MESH_FILE_EXT); - if ((fp=fopen(filename.c_str(),"w")) == NULL){ + filename.append(baseFileName.string()).append(MESH_FILE_EXT); + if ((fp=fopen(filename.c_str(),"w")) == nullptr){ std::string errorMessage; errorMessage.append("cannot open mesh file ").append(filename).append(" for writing"); - throw errorMessage.c_str(); + throw runtime_error(errorMessage); } simulation->getMesh()->write(fp); fclose(fp); - filename.erase(strlen(baseFileName), 5); + filename.erase(strlen(baseFileName.string().c_str()), 5); filename.append(MESHMETRICS_FILE_EXT); - if ((fp=fopen(filename.c_str(),"w")) == NULL){ + if ((fp=fopen(filename.c_str(),"w")) == nullptr){ std::string errorMessage; errorMessage.append("cannot open mesh metrics file ").append(filename).append(" for writing"); - throw errorMessage.c_str(); + throw runtime_error(errorMessage); } simulation->getMesh()->writeMeshMetrics(fp); fclose(fp); @@ -892,12 +840,12 @@ void SimTool::start1() { if (bCheckSpatiallyUniform) { increment = 0.001; } - double epsilon = 1e-12; clock_t oldTime = clock(); // to control the output of progress, send progress every 2 seconds. int counter = 0; while (true) { - if (simulation->getTime_sec() - simEndTime > epsilon // currentTime past endTime - || fabs(simEndTime - simulation->getTime_sec()) < epsilon) { // reached the end time + double epsilon = 1e-12; + if (simulation->getTime_sec(this) - simEndTime > epsilon // currentTime past endTime + || fabs(simEndTime - simulation->getTime_sec(this)) < epsilon) { // reached the end time break; } if (isSundialsPdeSolver() && bSundialsOneStepOutput) { @@ -911,8 +859,8 @@ void SimTool::start1() { } counter++; - simulation->iterate(); - if (smoldynSim != NULL && counter == smoldynStepMultiplier) { + simulation->iterate(this); + if (smoldynSim != nullptr && counter == smoldynStepMultiplier) { vcellhybrid::smoldynOneStep(smoldynSim);//smoldynOneStep includes computeHistogram after each time step. copyParticleCountsToConcentration(); counter = 0; @@ -922,19 +870,19 @@ void SimTool::start1() { return; } - simulation->update(); + simulation->update(this); - if (simulation->getCurrIteration() % keepEvery == 0 || fabs(simEndTime - simulation->getTime_sec()) < epsilon) { - updateLog(percentile,simulation->getTime_sec(), simulation->getCurrIteration()); + if (simulation->getCurrIteration() % keepEvery == 0 || fabs(simEndTime - simulation->getTime_sec(this)) < epsilon) { + updateLog(percentile,simulation->getTime_sec(this), simulation->getCurrIteration()); } clock_t currentTime = clock(); double duration = (double)(currentTime - oldTime) / CLOCKS_PER_SEC; if (duration >= 2){ - percentile = (simulation->getTime_sec() - simStartTime)/(simEndTime - simStartTime); + percentile = (simulation->getTime_sec(this) - simStartTime)/(simEndTime - simStartTime); if (percentile - lastSentPercentile >= increment) { - std::cout << "SimTool.start1() sending JOB_PROGRESS to SimulationMessaging: percentile=" << percentile << ", time=" << simulation->getTime_sec() << std::endl; - SimulationMessaging::getInstVar()->setWorkerEvent(new WorkerEvent(JOB_PROGRESS, percentile, simulation->getTime_sec())); + std::cout << "SimTool.start1() sending JOB_PROGRESS to SimulationMessaging: percentile=" << percentile << ", time=" << simulation->getTime_sec(this) << std::endl; + SimulationMessaging::getInstVar()->setWorkerEvent(new WorkerEvent(JOB_PROGRESS, percentile, simulation->getTime_sec(this))); lastSentPercentile = percentile; oldTime = currentTime; } @@ -949,24 +897,24 @@ void SimTool::start1() { } } if (uniform) { - cout << endl << "!!!Spatially Uniform reached at " << simulation->getTime_sec() << endl; + cout << endl << "!!!Spatially Uniform reached at " << simulation->getTime_sec(this) << endl; if (simulation->getCurrIteration() % keepEvery != 0) { - updateLog(percentile,simulation->getTime_sec(), simulation->getCurrIteration()); + updateLog(percentile,simulation->getTime_sec(this), simulation->getCurrIteration()); } break; } } } - if (smoldynSim != NULL) { + if (smoldynSim != nullptr) { vcellhybrid::smoldynEnd(smoldynSim); } if (checkStopRequested()) { return; } - SimulationMessaging::getInstVar()->setWorkerEvent(new WorkerEvent(JOB_PROGRESS, 1.0, simulation->getTime_sec())); - SimulationMessaging::getInstVar()->setWorkerEvent(new WorkerEvent(JOB_COMPLETED, percentile, simulation->getTime_sec())); + SimulationMessaging::getInstVar()->setWorkerEvent(new WorkerEvent(JOB_PROGRESS, 1.0, simulation->getTime_sec(this))); + SimulationMessaging::getInstVar()->setWorkerEvent(new WorkerEvent(JOB_COMPLETED, percentile, simulation->getTime_sec(this))); showSummary(stdout); } @@ -976,7 +924,7 @@ void SimTool::copyParticleCountsToConcentration(){ //concentration in terms of p int numVars = simulation->getNumVariables(); for (int i=0; igetVariable(i)->getVarType() == VAR_VOLUME_PARTICLE){ - VolumeParticleVariable* var = (VolumeParticleVariable*)simulation->getVariable(i); + auto* var = (VolumeParticleVariable*)simulation->getVariable(i); for(int j=0; jgetMesh()->getNumVolumeElements(); j++){ double count = var->getMoleculeCounts()[j]; if(count>0){ @@ -987,7 +935,7 @@ void SimTool::copyParticleCountsToConcentration(){ //concentration in terms of p } } else if(simulation->getVariable(i)->getVarType() == VAR_MEMBRANE_PARTICLE){ - MembraneParticleVariable* var = (MembraneParticleVariable*)simulation->getVariable(i); + auto* var = (MembraneParticleVariable*)simulation->getVariable(i); for(int j=0; jgetMesh()->getNumMembraneElements(); j++){ double count = var->getMoleculeCounts()[j]; if(count>0){ @@ -1003,7 +951,7 @@ void SimTool::copyParticleCountsToConcentration(){ //concentration in terms of p bool SimTool::checkSpatiallyUniform(Variable* var) { double* currSol = var->getCurr(); - CartesianMesh* mesh = (CartesianMesh*)simulation->getMesh(); + CartesianMesh* mesh = simulation->getMesh(); switch (var->getVarType()) { case VAR_VOLUME: for (int r = 0; r < mesh->getNumVolumeRegions(); r ++) { diff --git a/VCell/src/Simulation.cpp b/VCell/src/Simulation.cpp index 9e82e0afb..43aadc2ac 100644 --- a/VCell/src/Simulation.cpp +++ b/VCell/src/Simulation.cpp @@ -3,22 +3,17 @@ * All rights reserved. */ #include -#include #include -#include #include #include #include -#include #include #include #include -#include -#include #include #include -Simulation::Simulation(Mesh *mesh) +Simulation::Simulation(CartesianMesh *mesh) { ASSERTION(mesh); // @@ -60,9 +55,9 @@ void Simulation::advanceTimeOff() { _advanced=false; } -void Simulation::iterate() +void Simulation::iterate(SimTool* sim_tool) { - _scheduler->iterate(); + _scheduler->iterate(sim_tool); } void Simulation::update() @@ -154,51 +149,22 @@ void Simulation::readData(const char *filename) //------------------------------------------------------- // determines scheduler and resets _time_sec and initializes var's //------------------------------------------------------- -void Simulation::resolveReferences() { - VCellModel *model = SimTool::getInstance()->getModel(); - model->resolveReferences(); +void Simulation::resolveReferences(SimTool* sim_tool) { + sim_tool->getModel()->resolveReferences(sim_tool->getSimulation()); if (postProcessingBlock != NULL) { postProcessingBlock->resolveReferences(); } } -void Simulation::initSimulation() -{ - if (_scheduler == 0) { - int odeCount = 0, pdeCount = 0; - for (int i = 0; i < (int)varList.size(); i ++) { - Variable* var = varList[i]; - if (var->isDiffusing()){ - pdeCount ++; - } else { - odeCount ++; - } - } - - printf("pdeCount=%d, odeCount=%d\n", pdeCount, odeCount); - - SimTool* simTool = SimTool::getInstance(); - if (simTool->isSundialsPdeSolver()) { - _scheduler = new SundialsPdeScheduler(this, simTool->getSundialsSolverOptions(), - simTool->getNumDiscontinuityTimes(), simTool->getDiscontinuityTimes(), simTool->isSundialsOneStepOutput()); - } else { - _scheduler = new SerialScheduler(this); - } - } - - _scheduler->initValues(); - currIteration = 0; -} - -double Simulation::getTime_sec() { - if (SimTool::getInstance()->isSundialsPdeSolver()) { +double Simulation::getTime_sec(SimTool* sim_tool) { + if (sim_tool->isSundialsPdeSolver()) { return ((SundialsPdeScheduler*)_scheduler)->getCurrentTime(); } return _advanced ? (currIteration + 1) * _dT_sec : currIteration * _dT_sec; } -void Simulation::setSimStartTime(double st) { - if (SimTool::getInstance()->isSundialsPdeSolver()) { +void Simulation::setSimStartTime(SimTool* sim_tool, double st) { + if (sim_tool->isSundialsPdeSolver()) { ((SundialsPdeScheduler*)_scheduler)->setSimStartTime(st); } } diff --git a/VCell/src/SimulationExpression.cpp b/VCell/src/SimulationExpression.cpp index b75a1882b..ef54f4b82 100644 --- a/VCell/src/SimulationExpression.cpp +++ b/VCell/src/SimulationExpression.cpp @@ -26,6 +26,9 @@ #include #include + +#include "VCELL/SerialScheduler.h" +#include "VCELL/SundialsPdeScheduler.h" using std::endl; using std::stringstream; @@ -88,7 +91,7 @@ class VolumeRegionMembraneValueProxy : public ValueProxy Feature* feature; }; -SimulationExpression::SimulationExpression(Mesh *mesh) : Simulation(mesh) { +SimulationExpression::SimulationExpression(CartesianMesh *mesh) : Simulation(mesh) { symbolTable = NULL; indices = new int[NUM_VAR_INDEX]; @@ -222,19 +225,19 @@ void SimulationExpression::addMembraneRegionVariable(MembraneRegionVariable *var memRegionVarSize ++; } -void SimulationExpression::advanceTimeOn() { +void SimulationExpression::advanceTimeOn(SimTool* sim_tool) { Simulation::advanceTimeOn(); - valueProxyTime->setValue(getTime_sec()); + valueProxyTime->setValue(getTime_sec(sim_tool)); } -void SimulationExpression::advanceTimeOff() { +void SimulationExpression::advanceTimeOff(SimTool* sim_tool) { Simulation::advanceTimeOff(); - valueProxyTime->setValue(getTime_sec()); + valueProxyTime->setValue(getTime_sec(sim_tool)); } -void SimulationExpression::update() { +void SimulationExpression::update(SimTool* sim_tool) { Simulation::update(); - valueProxyTime->setValue(getTime_sec()); + valueProxyTime->setValue(getTime_sec(sim_tool)); } void SimulationExpression::addParameter(string& param) { @@ -242,18 +245,18 @@ void SimulationExpression::addParameter(string& param) { paramValueProxies.push_back(new ScalarValueProxy()); } -void SimulationExpression::createSymbolTable() { +void SimulationExpression::createSymbolTable(SimTool* sim_tool) { if (symbolTable != NULL) { return; } - if (SimTool::getInstance()->isSundialsPdeSolver()) { + if (sim_tool->isSundialsPdeSolver()) { if (volParticleVarSize > 0 || memParticleVarSize > 0) { throw "Fully implicit solver does not support hybrid simulations"; } } CartesianMesh *mesh = (CartesianMesh*)_mesh; - VCellModel* model = SimTool::getInstance()->getModel(); + VCellModel* model = sim_tool->getModel(); // volume size, membrane size, for each feature numRegionSizeVars = 1 + 1 + model->getNumFeatures(); @@ -357,7 +360,7 @@ void SimulationExpression::createSymbolTable() { } } - bool bSundialsPdeSolver = SimTool::getInstance()->isSundialsPdeSolver(); + bool bSundialsPdeSolver = sim_tool->isSundialsPdeSolver(); // t, x, y, z, VOL, VOL_Feature1_membrane, VOL_Feature2_membrane, ... (for computing membrane flux for volume variables), // MEM, VOLREG, VOLREG_Feature1_membrane, VOLREG_Feature2_membrane, ... (for computing membrane flux for volume region variables), @@ -495,7 +498,7 @@ void SimulationExpression::createSymbolTable() { // add random variable if (randomVarList.size() > 0) { char rvfile[512]; - sprintf(rvfile, "%s%s", SimTool::getInstance()->getBaseFileName(), RANDOM_VARIABLE_FILE_EXTENSION); + sprintf(rvfile, "%s%s", sim_tool->getBaseFileName().c_str(), RANDOM_VARIABLE_FILE_EXTENSION); FVDataSet::readRandomVariables(rvfile, this); for (int i = 0; i < (int)randomVarList.size(); i ++) { RandomVariable* rv = randomVarList[i]; @@ -527,7 +530,36 @@ void SimulationExpression::createSymbolTable() { numSymbols = variableIndex; symbolTable = new SimpleSymbolTable(variableNames, variableIndex, oldValueProxies); delete[] variableNames; -} +} + +void SimulationExpression::initSimulation(SimTool* sim_tool) +{ + if (_scheduler == 0) { + int odeCount = 0, pdeCount = 0; + for (int i = 0; i < (int)varList.size(); i ++) { + Variable* var = varList[i]; + if (var->isDiffusing()){ + pdeCount ++; + } else { + odeCount ++; + } + } + + printf("pdeCount=%d, odeCount=%d\n", pdeCount, odeCount); + + if (sim_tool->isSundialsPdeSolver()) { + _scheduler = new SundialsPdeScheduler(this, sim_tool->getSundialsSolverOptions(), + sim_tool->getNumDiscontinuityTimes(), sim_tool->getDiscontinuityTimes(), sim_tool->isSundialsOneStepOutput()); + } else { + _scheduler = new SerialScheduler(this); + } + } + + _scheduler->initValues(sim_tool->getModel()); + currIteration = 0; +} + + string* SimulationExpression::getFieldSymbols() { string* symbols = new string[fieldDataList.size()]; @@ -574,14 +606,14 @@ void SimulationExpression::populateParameterValues(double* darray) { } } -void SimulationExpression::resolveReferences() { +void SimulationExpression::resolveReferences(SimTool* sim_tool) { if (symbolTable == 0) { - createSymbolTable(); + createSymbolTable(sim_tool); } - Simulation::resolveReferences(); + Simulation::resolveReferences(sim_tool); } -void SimulationExpression::setParameterValues(double* paramValues) { +void SimulationExpression::setParameterValues(SimTool* sim_tool, double* paramValues) { if (paramValues == 0) { if (paramList.size() != 0) { throw "SimulationExpression::setParameterValues(), empty values for parameters"; @@ -591,7 +623,7 @@ void SimulationExpression::setParameterValues(double* paramValues) { for (int i = 0; i < (int)paramList.size(); i ++) { paramValueProxies[i]->setValue(paramValues[i]); } - reinitConstantValues(); + reinitConstantValues(sim_tool); } RandomVariable* SimulationExpression::getRandomVariableFromName(char* varName) @@ -621,16 +653,14 @@ bool SimulationExpression::isParameter(string& symbol) { return false; } -void SimulationExpression::reinitConstantValues() { - VCellModel* model = SimTool::getInstance()->getModel(); - - for (int i = 0; i < model->getNumFeatures(); i ++) { - Feature* feature = model->getFeatureFromIndex(i); - feature->reinitConstantValues(); +void SimulationExpression::reinitConstantValues(SimTool* sim_tool) { + for (int i = 0; i < sim_tool->getModel()->getNumFeatures(); i ++) { + Feature* feature = sim_tool->getModel()->getFeatureFromIndex(i); + feature->reinitConstantValues(sim_tool->getSimulation()); } - for (int i = 0; i < model->getNumMembranes(); i ++) { - Membrane* membrane = model->getMembraneFromIndex(i); - membrane->reinitConstantValues(); + for (int i = 0; i < sim_tool->getModel()->getNumMembranes(); i ++) { + Membrane* membrane = sim_tool->getModel()->getMembraneFromIndex(i); + membrane->reinitConstantValues(sim_tool->getSimulation()); } } diff --git a/VCell/src/Solver.cpp b/VCell/src/Solver.cpp index 399275b39..9c402244a 100644 --- a/VCell/src/Solver.cpp +++ b/VCell/src/Solver.cpp @@ -12,7 +12,7 @@ Solver::Solver(Variable *variable) eqnBuilder = NULL; } -void Solver::initEqn(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime) +void Solver::initEqn(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime) { if (!bFirstTime) { return; @@ -20,12 +20,12 @@ void Solver::initEqn(double deltaTime, int volumeIndexStart, int volumeIndexSize ASSERTION(eqnBuilder); - eqnBuilder->initEquation(deltaTime, volumeIndexStart, volumeIndexSize, membraneIndexStart, membraneIndexSize); + eqnBuilder->initEquation(model, deltaTime, volumeIndexStart, volumeIndexSize, membraneIndexStart, membraneIndexSize); } -void Solver::buildEqn(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime) +void Solver::buildEqn(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime) { ASSERTION(eqnBuilder); - eqnBuilder->buildEquation(deltaTime, volumeIndexStart, volumeIndexSize, membraneIndexStart, membraneIndexSize); + eqnBuilder->buildEquation(sim, deltaTime, volumeIndexStart, volumeIndexSize, membraneIndexStart, membraneIndexSize); } diff --git a/src/SolverMain.cpp b/VCell/src/SolverMain.cpp similarity index 60% rename from src/SolverMain.cpp rename to VCell/src/SolverMain.cpp index 1b8308e39..65b24082d 100644 --- a/src/SolverMain.cpp +++ b/VCell/src/SolverMain.cpp @@ -6,7 +6,7 @@ #include #endif -#include "SolverMain.h" +#include "../include/VCELL/SolverMain.h" #include #include @@ -16,25 +16,19 @@ using namespace std; #undef USE_MESSAGING -#include +#include <../include/VCELL/FVSolver.h> #include -#include -#include -#include -#include -#include +#include <../include/VCELL/SimTool.h> +#include <../../VCellMessaging/include/VCELL/SimulationMessaging.h> +#include <../../VCellMessaging/include/VCELL/GitDescribe.h> +#include <../../ExpressionParser/Exception.h> +#include <../../bridgeVCellSmoldyn/vcellhybrid.h> void vcellExit(int returnCode, string& errorMsg) { - if (SimulationMessaging::getInstVar() == 0) { - if (returnCode != 0) { - cerr << errorMsg << endl; - } - } else if (!SimulationMessaging::getInstVar()->isStopRequested()) { - if (returnCode != 0) { - SimulationMessaging::getInstVar()->setWorkerEvent(new WorkerEvent(JOB_FAILURE, errorMsg.c_str())); - } + if (returnCode != 0) { + cerr << errorMsg << endl; } } @@ -43,11 +37,11 @@ std::string version() return "Finite Volume version " + std::string(g_GIT_DESCRIBE) + " with smoldyn version " + std::string(VERSION); } -int solve(const std::string& inputFilename, const std::string& outputDir) { +int solve(const std::string& inputFilename, const std::string& vcgFilename, const std::string& outputDir) { // Check if output directory exists, if not create it std::filesystem::path dirPath(outputDir); - if (!std::filesystem::exists(dirPath)) { - std::filesystem::create_directories(dirPath); + if (!exists(dirPath)) { + create_directories(dirPath); } // Open the input file @@ -56,6 +50,12 @@ int solve(const std::string& inputFilename, const std::string& outputDir) { throw std::runtime_error("Could not open input file: " + inputFilename); } + // Open the vcg file + std::ifstream vcgFile(vcgFilename); + if (!vcgFile.is_open()) { + throw std::runtime_error("Could not open vcg file: " + vcgFilename); + } + vcellhybrid::setHybrid(); //get smoldyn library in correct state int returnCode = 0; string errorMsg = "Exception : "; @@ -63,21 +63,26 @@ int solve(const std::string& inputFilename, const std::string& outputDir) { bool bSimZip = true; int taskID = 0; - SimulationMessaging::create(); + if (SimulationMessaging::getInstVar() == nullptr) + { + SimulationMessaging::create(); + } FVSolver* fvSolver = nullptr; + SimTool* sim_tool = nullptr; try { - fvSolver = new FVSolver(inputFile, taskID, outputDir.c_str(), bSimZip); + fvSolver = new FVSolver(outputDir.c_str()); + sim_tool = fvSolver->createSimTool(inputFile, vcgFile, taskID, bSimZip); inputFile.close(); - if(fvSolver->getNumVariables() == 0){ + if (FVSolver::getNumVariables(sim_tool) == 0){ //sims with no reactions and no diffusing species cause exit logic to 'wait' forever //never sending a job failed or job finished message and never cleaning up threads throw invalid_argument("FiniteVolume error: Must have at least 1 variable or reaction to solve"); } - fvSolver->solve(); + FVSolver::solve(sim_tool); } catch (const char *exStr){ errorMsg += exStr; @@ -100,6 +105,11 @@ int solve(const std::string& inputFilename, const std::string& outputDir) { inputFile.close(); } vcellExit(returnCode, errorMsg); - delete fvSolver; + if (sim_tool != nullptr) { + delete sim_tool; + } + if (fvSolver != nullptr) { + delete fvSolver; + } return returnCode; } diff --git a/VCell/src/SparseLinearSolver.cpp b/VCell/src/SparseLinearSolver.cpp index 7dc0be3c2..fafcd2c42 100644 --- a/VCell/src/SparseLinearSolver.cpp +++ b/VCell/src/SparseLinearSolver.cpp @@ -85,11 +85,11 @@ SparseLinearSolver::~SparseLinearSolver() delete[] pcg_workspace; } -void SparseLinearSolver::solveEqn(double dT_sec, +void SparseLinearSolver::solveEqn(SimTool* sim_tool, double dT_sec, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize, bool bFirstTime) { - int* IParm = PCGSolve(bFirstTime); + int* IParm = PCGSolve(sim_tool, bFirstTime); int returnCode = IParm[50]; int additionalSpace = IParm[53]; delete[] IParm; @@ -106,7 +106,7 @@ void SparseLinearSolver::solveEqn(double dT_sec, enableRetry = false; cout << endl << "!!Note: Insufficient PCG workspace (" << nWork << "), need additional (" << additionalSpace << "), for variable " << var->getName() << ", try again" << endl; initPCGWorkspace(additionalSpace); - solveEqn(dT_sec, volumeIndexStart, volumeIndexSize, membraneIndexStart, membraneIndexSize, true); + solveEqn(sim_tool, dT_sec, volumeIndexStart, volumeIndexSize, membraneIndexStart, membraneIndexSize, true); } else { throwPCGExceptions(returnCode, additionalSpace); // this throws exception } @@ -120,7 +120,7 @@ void SparseLinearSolver::solveEqn(double dT_sec, } // -------------------------------------------------- -int* SparseLinearSolver::PCGSolve(bool bRecomputeIncompleteFactorization) +int* SparseLinearSolver::PCGSolve(SimTool* sim_tool, bool bRecomputeIncompleteFactorization) // -------------------------------------------------- { // prepare data to call pcgpak @@ -137,7 +137,7 @@ int* SparseLinearSolver::PCGSolve(bool bRecomputeIncompleteFactorization) long size = smEqnBuilder->getSize(); string timername = var->getName() + " PCG"; - TimerHandle tHndPCG = SimTool::getInstance()->getTimerHandle(timername); + TimerHandle tHndPCG = sim_tool->getTimerHandle(timername); int* IParm = new int[75]; memset(IParm, 0, 75 * sizeof(int)); @@ -169,14 +169,14 @@ int* SparseLinearSolver::PCGSolve(bool bRecomputeIncompleteFactorization) // Set tolerance double PCG_Tolerance = pcgRelErr; - SimTool::getInstance()->startTimer(tHndPCG); + sim_tool->startTimer(tHndPCG); string varname = var->getName(); double RHSscale = computeRHSscale(size, pRHS, varname); #ifdef SHOW_MATRIX cout << setprecision(10); - cout << "----SparseLinearSolver----Variable " << var->getName() << " at " << SimTool::getInstance()->getSimulation()->getTime_sec() << "---------------" << endl; + cout << "----SparseLinearSolver----Variable " << var->getName() << " at " << sim_tool->getSimulation()->getTime_sec() << "---------------" << endl; A->show(); cout << "--------SparseLinearSolver----RHS-----------" << endl; for (int index = 0; index < size; index++){ @@ -187,7 +187,7 @@ int* SparseLinearSolver::PCGSolve(bool bRecomputeIncompleteFactorization) #endif PCGWRAPPER(&size, &Nrsp, &symmetricflg, ija, sa, pRHS, pNew, &PCG_Tolerance, IParm, RParm, pcg_workspace, pcg_workspace, &RHSscale); - SimTool::getInstance()->stopTimer(tHndPCG); + sim_tool->stopTimer(tHndPCG); #ifdef SHOW_IPARM cout << endl << "-------SparseLinearSolver----numNonZeros=" << A->getNumNonZeros() << "--------------------" << endl; diff --git a/VCell/src/SparseVolumeEqnBuilder.cpp b/VCell/src/SparseVolumeEqnBuilder.cpp index b4bd820df..2edf54dae 100644 --- a/VCell/src/SparseVolumeEqnBuilder.cpp +++ b/VCell/src/SparseVolumeEqnBuilder.cpp @@ -438,7 +438,7 @@ void SparseVolumeEqnBuilder::computeLHS(int index, double* lambdas, double& Aii, // Left Hand Side // //------------------------------------------------------------------ -void SparseVolumeEqnBuilder::initEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) +void SparseVolumeEqnBuilder::initEquation(VCellModel* model, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) { /** * we decided to scale both sides with deltaT/VOLUME @@ -461,7 +461,7 @@ void SparseVolumeEqnBuilder::initEquation(double deltaTime, int volumeIndexStart double lambdas[] = {LAMBDAX, LAMBDAY, LAMBDAZ, LAMBDAAREAX, LAMBDAAREAY, LAMBDAAREAZ}; if (!bPreProcessed) { - preProcess(); + preProcess(model); } ASSERTION(solver->getVar() == var); @@ -542,10 +542,10 @@ void SparseVolumeEqnBuilder::initEquation(double deltaTime, int volumeIndexStart delete[] columnValues; } -double SparseVolumeEqnBuilder::computeRHS(int index, double deltaTime, double* lambdas, double bInit) { +double SparseVolumeEqnBuilder::computeRHS(Simulation* sim, int index, double deltaTime, double* lambdas, double bInit) { string varname = var->getName(); double b = bInit; - Simulation* sim = SimTool::getInstance()->getSimulation(); + VolumeElement *pVolumeElement = mesh->getVolumeElements(); MembraneElement *pMembraneElement = mesh->getMembraneElements(); @@ -725,7 +725,7 @@ double SparseVolumeEqnBuilder::computeRHS(int index, double deltaTime, double* l // Right Hand side // //------------------------------------------------------------------ -void SparseVolumeEqnBuilder::buildEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) +void SparseVolumeEqnBuilder::buildEquation(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) { double lambdaAreaX = deltaTime/DELTAX; double lambdaAreaY = deltaTime/DELTAY; @@ -735,7 +735,7 @@ void SparseVolumeEqnBuilder::buildEquation(double deltaTime, int volumeIndexStar if (bSolveWholeMesh) { memcpy(B, var->getCurr(), var->getSize() * sizeof(double)); for (int index = volumeIndexStart; index < volumeIndexStart + volumeIndexSize; index ++){ - B[index] = computeRHS(index, deltaTime, lambdas, B[index]); + B[index] = computeRHS(sim, index, deltaTime, lambdas, B[index]); } } else { double* currVal = var->getCurr(); @@ -744,7 +744,7 @@ void SparseVolumeEqnBuilder::buildEquation(double deltaTime, int volumeIndexStar // to initialize X, which will be passed to solver as intial guess and final solution. // or set initial guess to zero (need to revisit, we also want to revisit fill-in parameter) X[localIndex] = currVal[globalIndex]; - B[localIndex] = computeRHS(globalIndex, deltaTime, lambdas, X[localIndex]); + B[localIndex] = computeRHS(sim, globalIndex, deltaTime, lambdas, X[localIndex]); } } // to make the matrix symmetric @@ -781,7 +781,7 @@ bool SparseVolumeEqnBuilder::checkPeriodicCoupledPairsInRegions(int indexm, int #define PERIODIC_GEOMETRY_ERROR_MESSAGE "geometry is not compatible with periodic boundary conditions. Couldn't find corresponding volume element." -void SparseVolumeEqnBuilder::preProcess() { +void SparseVolumeEqnBuilder::preProcess(VCellModel* model) { if (bPreProcessed) { return; } @@ -790,7 +790,6 @@ void SparseVolumeEqnBuilder::preProcess() { // check if there is periodic boundary condition in the model bool bPeriodic = false; - VCellModel* model = SimTool::getInstance()->getModel(); for (int i = 0; i < model->getNumFeatures(); i ++) { Feature* feature = model->getFeatureFromIndex(i); if (feature->getXmBoundaryType() == BOUNDARY_PERIODIC diff --git a/VCell/src/Structure.cpp b/VCell/src/Structure.cpp index 99daad5c9..869a9c282 100644 --- a/VCell/src/Structure.cpp +++ b/VCell/src/Structure.cpp @@ -10,10 +10,6 @@ Structure::Structure(string& Aname) { } } -Structure::~Structure(void) -{ -} - int Structure::getNumElements() { if (numElements == 0) { int num = 0; diff --git a/VCell/src/SundialsPdeScheduler.cpp b/VCell/src/SundialsPdeScheduler.cpp index baa52d16f..91bffe9fc 100644 --- a/VCell/src/SundialsPdeScheduler.cpp +++ b/VCell/src/SundialsPdeScheduler.cpp @@ -71,10 +71,10 @@ static int plusMasks[3] = {NEIGHBOR_XP_MASK, NEIGHBOR_YP_MASK, NEIGHBOR_ZP_MASK} #define PRECOMPUTE_DIFFUSION_COEFFICIENT -SundialsPdeScheduler::SundialsPdeScheduler(Simulation *sim, const SundialsSolverOptions& sso, int numDisTimes, double* disTimes, bool bDefaultOuptput) : Scheduler(sim) +SundialsPdeScheduler::SundialsPdeScheduler(SimulationExpression *sim, const SundialsSolverOptions& sso, int numDisTimes, double* disTimes, bool bDefaultOuptput) : Scheduler(sim) { - simulation = (SimulationExpression*)sim; - mesh = (CartesianMesh*)simulation->getMesh(); + simulation = sim; + mesh = simulation->getMesh(); y = 0; sundialsSolverMemory = 0; @@ -141,12 +141,12 @@ SundialsPdeScheduler::~SundialsPdeScheduler() delete[] rhsGradients; } -void SundialsPdeScheduler::iterate() { +void SundialsPdeScheduler::iterate(SimTool* sim_tool) { if (bFirstTime) { setupOrderMaps(); - initSundialsSolver(); + initSundialsSolver(sim_tool->getModel()); } - solve(); + solve(sim_tool); bFirstTime = false; } @@ -419,7 +419,7 @@ int SundialsPdeScheduler::CVodeRHS(double t, double* yinput, double* rhs) { return 0; } -void SundialsPdeScheduler::initSundialsSolver() { +void SundialsPdeScheduler::initSundialsSolver(VCellModel* model) { if (!bFirstTime) { return; } @@ -429,7 +429,7 @@ void SundialsPdeScheduler::initSundialsSolver() { int numVolRegionVar = simulation->getNumVolRegionVariables(); int numMemRegionVar = simulation->getNumMemRegionVariables(); if (sundialsSolverMemory == 0) { - numSymbolsPerVolVar = SimTool::getInstance()->getModel()->getNumFeatures() + 1; + numSymbolsPerVolVar = model->getNumFeatures() + 1; // t, x, y, z, (U, U_Feature1_membrane, U_Feature2_membrane, ...), (M), // (VR, VR_Feature1_membrane, ...), (MR), (RegionSize), (FieldData), (RandomVariable), (Parameter) @@ -643,10 +643,10 @@ void SundialsPdeScheduler::printCVodeStats() printf("last step = %f\n\n", hlast); } -void SundialsPdeScheduler::solve() { +void SundialsPdeScheduler::solve(SimTool* sim_tool) { double endTime = 0; if (bSundialsOneStepOutput) { - endTime = SimTool::getInstance()->getEndTime(); + endTime = sim_tool->getEndTime(); } else { endTime = currentTime + simulation->getDT_sec(); } @@ -671,8 +671,8 @@ void SundialsPdeScheduler::solve() { if (bStop && fabs(stopTime - currentTime) < epsilon) { currDiscontinuityTimeCount ++; - if (SimTool::getInstance()->getEndTime() - currentTime > epsilon // currentTime less than endTime - || fabs(SimTool::getInstance()->getEndTime() - currentTime) > epsilon) { // if this is the last solve, we don't have to reinit + if (sim_tool->getEndTime() - currentTime > epsilon // currentTime less than endTime + || fabs(sim_tool->getEndTime() - currentTime) > epsilon) { // if this is the last solve, we don't have to reinit cout << endl << "SundialsPdeScheduler::solve() : cvode reinit at time " << currentTime << endl; returnCode = CVodeReInit(sundialsSolverMemory, RHS_callback, currentTime, y, ToleranceType, sundialsSolverOptions.relTol, &sundialsSolverOptions.absTol); @@ -698,8 +698,8 @@ void SundialsPdeScheduler::solve() { } } - if (currentTime - SimTool::getInstance()->getEndTime() > epsilon // currentTime greater than endTime - || (SimTool::getInstance()->getEndTime() - currentTime) < epsilon) { + if (currentTime - sim_tool->getEndTime() > epsilon // currentTime greater than endTime + || (sim_tool->getEndTime() - currentTime) < epsilon) { printCVodeStats(); } updateSolutions(); diff --git a/VCell/src/VCellModel.cpp b/VCell/src/VCellModel.cpp index eb4d2e0b4..8e5f74866 100644 --- a/VCell/src/VCellModel.cpp +++ b/VCell/src/VCellModel.cpp @@ -73,9 +73,8 @@ Feature *VCellModel::getFeatureFromName(string& name) throw ss.str(); } -void VCellModel::resolveReferences() -{ - Simulation* sim = SimTool::getInstance()->getSimulation(); +void VCellModel::resolveReferences(SimulationExpression* sim) +{ for (int i = 0; i < (int)featureList.size(); i ++) { Feature* feature = featureList[i]; feature->resolveReferences(sim); diff --git a/VCell/src/VarContext.cpp b/VCell/src/VarContext.cpp index 5f6d01da8..ae12c3904 100644 --- a/VCell/src/VarContext.cpp +++ b/VCell/src/VarContext.cpp @@ -13,41 +13,36 @@ using VCell::Expression; #include #include #include -#include -#include #include #include #include #include -#include +#include VarContext::VarContext(Structure *s, Variable* var) +: expressions(TOTAL_NUM_EXPRESSIONS, nullptr), + constantValues(TOTAL_NUM_EXPRESSIONS, nullptr), + dependencyMask(TOTAL_NUM_EXPRESSIONS, DEPENDENCY_MASK_UNDEFINED), + jumpConditionList(0) { structure = s; - ASSERTION(feature); + ASSERTION(structure != nullptr); species = var; - sim = NULL; - - expressions = new Expression*[TOTAL_NUM_EXPRESSIONS]; - constantValues = new double*[TOTAL_NUM_EXPRESSIONS]; - dependencyMask = new unsigned char[TOTAL_NUM_EXPRESSIONS]; - memset(expressions, 0, TOTAL_NUM_EXPRESSIONS * sizeof(Expression*)); - memset(constantValues, 0, TOTAL_NUM_EXPRESSIONS * sizeof(double*)); - memset(dependencyMask, DEPENDENCY_MASK_UNDEFINED, TOTAL_NUM_EXPRESSIONS * sizeof(unsigned char)); + sim = nullptr; } -void VarContext::resolveReferences(Simulation *Asim) +void VarContext::resolveReferences(SimulationExpression *Asim) { sim = Asim; - if (sim == 0) { - throw "VarContext::resolveReference(), simulation can't be null"; + if (sim == nullptr) { + throw std::runtime_error("VarContext::resolveReference(), simulation can't be null"); } } double VarContext::getInitialValue(long index) { - throw "Application Error: neither initialValue nor getInitialValue() specified for VarContext"; + throw std::runtime_error("Application Error: neither initialValue nor getInitialValue() specified for VarContext"); } VarContext::~VarContext() @@ -55,70 +50,70 @@ VarContext::~VarContext() for (int i = 0; i < TOTAL_NUM_EXPRESSIONS; i ++) { delete expressions[i]; delete constantValues[i]; - } - delete[] expressions; - delete[] constantValues; - delete[] dependencyMask; + } + expressions.clear(); + constantValues.clear(); + dependencyMask.clear(); - for (int i = 0; i < (int)jumpConditionList.size(); i ++) { - delete jumpConditionList[i]; + for (auto & i : jumpConditionList) { + delete i; } jumpConditionList.clear(); } void VarContext::setExpression(Expression* newexp, int expIndex) { - if (expressions[expIndex] != 0) { + if (expressions[expIndex] != nullptr) { stringstream ss; ss << "Expression " << String_Expression_Index[expIndex] << " for variable " << species->getName() << " in Structure " << structure->getName() << " has been set already"; - throw ss.str(); + throw std::runtime_error(ss.str()); } expressions[expIndex] = newexp; } -void VarContext::bindAll(SimulationExpression* simulation) { - SymbolTable* symbolTable = simulation->getSymbolTable(); +void VarContext::bindAll(SimulationExpression* sim) { + SymbolTable* symbolTable = sim->getSymbolTable(); for (int i = 0; i < TOTAL_NUM_EXPRESSIONS; i ++) { - if (expressions[i] == 0) { + if (expressions[i] == nullptr) { if (isNullExpressionOK(i)) { continue; } else { stringstream ss; ss << "VarContext::bindAll(), expression " << String_Expression_Index[i] << " for variable " << species->getName() << " not defined"; - throw ss.str(); + throw std::runtime_error(ss.str()); } } try { - //cout << expressions[i]->infix() << endl; + //cout << expressions[i]->infix() << std::endl; double d = expressions[i]->evaluateConstant(); constantValues[i] = new double[1]; constantValues[i][0] = d; } catch (...) { expressions[i]->bindExpression(symbolTable); } - computeDependencyMask(i); + computeDependencyMask(sim, i); } - for (int i = 0; i < (int)jumpConditionList.size(); i ++) { - jumpConditionList[i]->bindExpression(symbolTable); + for (auto & i : jumpConditionList) { + i->bindExpression(symbolTable); } } -double VarContext::evaluateExpression(MembraneElement* element, long expIndex) +double VarContext::evaluateExpression(MembraneElement* element, long expIndex) const { - if (expressions[expIndex] == 0) { // not defined - throw "VarContext::evaluateExpression(MembaneElement), expression not defined"; + if (expressions[expIndex] == nullptr) { // not defined + throw std::runtime_error("VarContext::evaluateExpression(MembaneElement), expression not defined"); } - if (constantValues[expIndex] != NULL) { + if (constantValues[expIndex] != nullptr) { return constantValues[expIndex][0]; } //cout << expressions[expIndex]->infix() << endl; if (dependencyMask[expIndex] & DEPENDENCY_MASK_XYZ) { WorldCoord wc = sim->getMesh()->getMembraneWorldCoord(element); - ((SimulationExpression*)sim)->setCurrentCoordinate(wc); + sim->setCurrentCoordinate(wc); } - int* indices = ((SimulationExpression*)sim)->getIndices(); + int* indices = sim->getIndices(); indices[VAR_VOLUME_INDEX] = -1; indices[VAR_VOLUME_REGION_INDEX] = -1; indices[VAR_MEMBRANE_INDEX] = element->index; @@ -126,31 +121,33 @@ double VarContext::evaluateExpression(MembraneElement* element, long expIndex) return expressions[expIndex]->evaluateProxy(); } -double VarContext::evaluateConstantExpression(long expIndex) { +double VarContext::evaluateConstantExpression(long expIndex) const +{ // pure constant - if (constantValues[expIndex] != 0) { + if (constantValues[expIndex] != nullptr) { return constantValues[expIndex][0]; } stringstream ss; ss << "VarContext::evaluateConstantExpression(), for variable " << species->getName() << " expression " << String_Expression_Index[expIndex] << " not defined OR not a constant expression"; - throw ss.str(); + throw std::runtime_error(ss.str()); } -double VarContext::evaluateExpression(long volIndex, long expIndex) { - if (expressions[expIndex] == 0) { // not defined +double VarContext::evaluateExpression(long volIndex, long expIndex) const +{ + if (expressions[expIndex] == nullptr) { // not defined stringstream ss; ss << "VarContext::evaluateExpression(VolIndex), for variable " << species->getName() << " expression " << String_Expression_Index[expIndex] << " not defined"; - throw ss.str(); + throw std::runtime_error(ss.str()); } - if (constantValues[expIndex] != NULL) { + if (constantValues[expIndex] != nullptr) { return constantValues[expIndex][0]; } if (dependencyMask[expIndex] & DEPENDENCY_MASK_XYZ) { WorldCoord wc = ((CartesianMesh*)sim->getMesh())->getVolumeWorldCoord(volIndex); - ((SimulationExpression*)sim)->setCurrentCoordinate(wc); + sim->setCurrentCoordinate(wc); } - int* indices = ((SimulationExpression*)sim)->getIndices(); + int* indices = sim->getIndices(); indices[VAR_MEMBRANE_INDEX] = -1; indices[VAR_MEMBRANE_REGION_INDEX] = -1; indices[VAR_VOLUME_INDEX] = volIndex; @@ -158,86 +155,89 @@ double VarContext::evaluateExpression(long volIndex, long expIndex) { return expressions[expIndex]->evaluateProxy(); } -double VarContext::evaluateExpression(long expIndex, double* values) { - if (expressions[expIndex] == 0) { // not defined +double VarContext::evaluateExpression(long expIndex, double* values) const +{ + if (expressions[expIndex] == nullptr) { // not defined stringstream ss; ss << "VarContext::evaluateExpression(), for variable " << species->getName() << " expression " << String_Expression_Index[expIndex] << " not defined"; - throw ss.str(); + throw std::runtime_error(ss.str()); } - if (constantValues[expIndex] != NULL) { + if (constantValues[expIndex] != nullptr) { return constantValues[expIndex][0]; } return expressions[expIndex]->evaluateVector(values); } void VarContext::addJumpCondition(Membrane* membrane, Expression* exp) { - JumpCondition* jc = new JumpCondition(membrane, exp); + auto* jc = new JumpCondition(membrane, exp); jumpConditionList.push_back(jc); } -double VarContext::evaluateJumpCondition(MembraneElement* element) +double VarContext::evaluateJumpCondition(MembraneElement* element) const { - for (int i = 0; i < (int)jumpConditionList.size(); i ++) { - if (jumpConditionList[i]->getMembrane() == element->getMembrane()) { - return jumpConditionList[i]->evaluateExpression(((SimulationExpression*)sim), element); + for (auto & i : jumpConditionList) { + if (i->getMembrane() == element->getMembrane()) { + return i->evaluateExpression(sim, element); } } stringstream ss; ss << "Jump Condition for variable " << species->getName() << " in Feature " << structure->getName() << " not found for Membrane " << element->getMembrane()->getName(); - throw ss.str(); + throw std::runtime_error(ss.str()); } -double VarContext::evaluateJumpCondition(MembraneElement* element, double* values) +double VarContext::evaluateJumpCondition(MembraneElement* element, double* values) const { - for (int i = 0; i < (int)jumpConditionList.size(); i ++) { - if (jumpConditionList[i]->getMembrane() == element->getMembrane()) { - return jumpConditionList[i]->evaluateExpression(values); + for (auto & i : jumpConditionList) { + if (i->getMembrane() == element->getMembrane()) { + return i->evaluateExpression(values); } } stringstream ss; ss << "Jump Condition for variable " << species->getName() << " in Feature " << structure->getName() << " not found for Membrane " << element->getMembrane()->getName(); - throw ss.str(); + throw std::runtime_error(ss.str()); } -bool VarContext::isConstantExpression(long expIndex) { +bool VarContext::isConstantExpression(long expIndex) const +{ if (dependencyMask[expIndex] == DEPENDENCY_MASK_CONSTANT) { return true; } // not defined - if (expressions[expIndex] == NULL) { + if (expressions[expIndex] == nullptr) { stringstream ss; ss << "VarContext::isConstantExpression(), for variable " << species->getName() << " expression " << String_Expression_Index[expIndex] << " not defined"; - throw ss.str(); + throw std::runtime_error(ss.str()); } return false; } -bool VarContext::isXYZOnlyExpression(long expIndex) { +bool VarContext::isXYZOnlyExpression(long expIndex) const +{ if (dependencyMask[expIndex] == DEPENDENCY_MASK_XYZ) { return true; } // not defined - if (expressions[expIndex] == NULL) { + if (expressions[expIndex] == nullptr) { stringstream ss; ss << "VarContext::isConstantExpression(), for variable " << species->getName() << " expression " << String_Expression_Index[expIndex] << " not defined"; - throw ss.str(); + throw std::runtime_error(ss.str()); } return false; } -void VarContext::computeDependencyMask(int expIndex) { +void VarContext::computeDependencyMask(SimulationExpression* sim, int expIndex) { // not defined - if (constantValues[expIndex] != NULL) { // constant 0 + if (constantValues[expIndex] != nullptr) { // constant 0 dependencyMask[expIndex] = DEPENDENCY_MASK_CONSTANT; return; } - if (expressions[expIndex] == NULL) { // expression is not defined + if (expressions[expIndex] == nullptr) { // expression is not defined return; } @@ -246,13 +246,12 @@ void VarContext::computeDependencyMask(int expIndex) { bool bHasTime = false; vector symbols; expressions[expIndex]->getSymbols(symbols); - SimulationExpression* simExp = (SimulationExpression*)SimTool::getInstance()->getSimulation(); - for (int i = 0; i < (int)symbols.size(); i ++) { - if (symbols[i] == "t") { + for (auto & symbol : symbols) { + if (symbol == "t") { bHasTime = true; - } else if (simExp->isVariable(symbols[i])) { + } else if (sim->isVariable(symbol)) { bHasVariable = true; - } else if (!simExp->isParameter(symbols[i])) { + } else if (!sim->isParameter(symbol)) { bHasXYZ = true; } } @@ -272,9 +271,9 @@ void VarContext::computeDependencyMask(int expIndex) { } } -double VarContext::evaluateMembraneRegionExpression(long memRegionIndex, long expIndex) +double VarContext::evaluateMembraneRegionExpression(long memRegionIndex, long expIndex) const { - int* indices = ((SimulationExpression*)sim)->getIndices(); + int* indices = sim->getIndices(); indices[VAR_VOLUME_INDEX] = -1; indices[VAR_VOLUME_REGION_INDEX] = -1; indices[VAR_MEMBRANE_INDEX] = -1; @@ -282,9 +281,9 @@ double VarContext::evaluateMembraneRegionExpression(long memRegionIndex, long ex return expressions[expIndex]->evaluateProxy(); } -double VarContext::evaluateVolumeRegionExpression(long volRegionIndex, long expIndex) +double VarContext::evaluateVolumeRegionExpression(long volRegionIndex, long expIndex) const { - int* indices = ((SimulationExpression*)sim)->getIndices(); + int* indices = sim->getIndices(); indices[VAR_VOLUME_INDEX] = -1; indices[VAR_MEMBRANE_INDEX] = -1; indices[VAR_MEMBRANE_REGION_INDEX] = -1; @@ -292,19 +291,19 @@ double VarContext::evaluateVolumeRegionExpression(long volRegionIndex, long expI return expressions[expIndex]->evaluateProxy(); } -void VarContext::reinitConstantValues() { +void VarContext::reinitConstantValues(SimulationExpression* sim) { for (int i = 0; i < TOTAL_NUM_EXPRESSIONS; i ++) { - if (expressions[i] == 0 || !isConstantExpression(i)) { + if (expressions[i] == nullptr || !isConstantExpression(i)) { continue; } double d = expressions[i]->evaluateProxy(); - if (constantValues[i] == 0) { + if (constantValues[i] == nullptr) { constantValues[i] = new double[1]; } constantValues[i][0] = d; } - for (int i = 0; i < (int)jumpConditionList.size(); i ++) { - jumpConditionList[i]->reinitConstantValues(); + for (auto & i : jumpConditionList) { + i->reinitConstantValues(sim); } } diff --git a/VCell/src/VolumeRegionEqnBuilder.cpp b/VCell/src/VolumeRegionEqnBuilder.cpp index a0cfb4170..bfd2d8f4e 100644 --- a/VCell/src/VolumeRegionEqnBuilder.cpp +++ b/VCell/src/VolumeRegionEqnBuilder.cpp @@ -17,7 +17,7 @@ VolumeRegionEqnBuilder::VolumeRegionEqnBuilder(VolumeRegionVariable *Avar, Carte odeSolver = Asolver; } -void VolumeRegionEqnBuilder::buildEquation(double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) { +void VolumeRegionEqnBuilder::buildEquation(Simulation* sim, double deltaTime, int volumeIndexStart, int volumeIndexSize, int membraneIndexStart, int membraneIndexSize) { int size = ((CartesianMesh*)mesh)->getNumVolumeRegions(); ASSERTION(size==var->getSize()); double *pRate = odeSolver->getRates(); diff --git a/VCell/src/VolumeRegionVarContextExpression.cpp b/VCell/src/VolumeRegionVarContextExpression.cpp index 08614ed71..ff61aa44f 100644 --- a/VCell/src/VolumeRegionVarContextExpression.cpp +++ b/VCell/src/VolumeRegionVarContextExpression.cpp @@ -3,7 +3,6 @@ * All rights reserved. */ #include -#include #include #include #include @@ -13,28 +12,31 @@ VolumeRegionVarContextExpression:: VolumeRegionVarContextExpression(Feature *fea { } -void VolumeRegionVarContextExpression::resolveReferences(Simulation* sim) { +void VolumeRegionVarContextExpression::resolveReferences(SimulationExpression* sim) { VarContext::resolveReferences(sim); - bindAll((SimulationExpression*)sim); + bindAll(sim); } double VolumeRegionVarContextExpression::getInitialValue(long regionIndex) { return evaluateVolumeRegionExpression(regionIndex, INITIAL_VALUE_EXP); } -double VolumeRegionVarContextExpression::getReactionRate(long volumeIndex) { +double VolumeRegionVarContextExpression::getReactionRate(long volumeIndex) const +{ return evaluateExpression(volumeIndex, REACT_RATE_EXP); } -double VolumeRegionVarContextExpression::getUniformRate(VolumeRegion *region){ +double VolumeRegionVarContextExpression::getUniformRate(VolumeRegion *region) const +{ return evaluateVolumeRegionExpression(region->getIndex(), UNIFORM_RATE_EXP); } -double VolumeRegionVarContextExpression::getFlux(MembraneElement *element){ +double VolumeRegionVarContextExpression::getFlux(MembraneElement *element) const +{ return evaluateJumpCondition(element); } -bool VolumeRegionVarContextExpression::isNullExpressionOK(int expIndex) { +bool VolumeRegionVarContextExpression::isNullExpressionOK(int expIndex) const { if (expIndex == INITIAL_VALUE_EXP || expIndex == REACT_RATE_EXP || expIndex == UNIFORM_RATE_EXP) { return false; } diff --git a/VCell/src/VolumeVarContextExpression.cpp b/VCell/src/VolumeVarContextExpression.cpp index f5baf0d45..5167d28c3 100644 --- a/VCell/src/VolumeVarContextExpression.cpp +++ b/VCell/src/VolumeVarContextExpression.cpp @@ -2,12 +2,11 @@ * (C) Copyright University of Connecticut Health Center 2001. * All rights reserved. */ +#include #include #include #include -#include #include -#include #include #include @@ -16,109 +15,128 @@ VolumeVarContextExpression::VolumeVarContextExpression(Feature *feature, VolumeV { } -void VolumeVarContextExpression::resolveReferences(Simulation* sim) { +void VolumeVarContextExpression::resolveReferences(SimulationExpression* sim) { VarContext::resolveReferences(sim); - bindAll((SimulationExpression*)sim); + bindAll(sim); } double VolumeVarContextExpression::getInitialValue(long volIndex) { return evaluateExpression(volIndex, INITIAL_VALUE_EXP); } -double VolumeVarContextExpression::getDiffusionRate(long volIndex) +double VolumeVarContextExpression::getDiffusionRate(long volIndex) const { return evaluateExpression(volIndex, DIFF_RATE_EXP); } -double VolumeVarContextExpression::getReactionRate(long volIndex) +double VolumeVarContextExpression::getReactionRate(long volIndex) const { return evaluateExpression(volIndex, REACT_RATE_EXP); } -double VolumeVarContextExpression::getXmBoundaryValue(long volIndex){ +double VolumeVarContextExpression::getXmBoundaryValue(long volIndex) const +{ return evaluateExpression(volIndex, BOUNDARY_XM_EXP); } -double VolumeVarContextExpression::getXpBoundaryValue(long volIndex) { +double VolumeVarContextExpression::getXpBoundaryValue(long volIndex) const +{ return evaluateExpression(volIndex, BOUNDARY_XP_EXP); } -double VolumeVarContextExpression::getYmBoundaryValue(long volIndex) { +double VolumeVarContextExpression::getYmBoundaryValue(long volIndex) const +{ return evaluateExpression(volIndex, BOUNDARY_YM_EXP); } -double VolumeVarContextExpression::getYpBoundaryValue(long volIndex){ +double VolumeVarContextExpression::getYpBoundaryValue(long volIndex) const +{ return evaluateExpression(volIndex, BOUNDARY_YP_EXP); } -double VolumeVarContextExpression::getZmBoundaryValue(long volIndex){ +double VolumeVarContextExpression::getZmBoundaryValue(long volIndex) const +{ return evaluateExpression(volIndex, BOUNDARY_ZM_EXP); } -double VolumeVarContextExpression::getZpBoundaryValue(long volIndex){ +double VolumeVarContextExpression::getZpBoundaryValue(long volIndex) const +{ return evaluateExpression(volIndex, BOUNDARY_ZP_EXP); } -double VolumeVarContextExpression::getXmBoundaryFlux(long volIndex){ +double VolumeVarContextExpression::getXmBoundaryFlux(long volIndex) const +{ return evaluateExpression(volIndex, BOUNDARY_XM_EXP); } -double VolumeVarContextExpression::getXpBoundaryFlux(long volIndex){ +double VolumeVarContextExpression::getXpBoundaryFlux(long volIndex) const +{ return evaluateExpression(volIndex, BOUNDARY_XP_EXP); } -double VolumeVarContextExpression::getYmBoundaryFlux(long volIndex){ +double VolumeVarContextExpression::getYmBoundaryFlux(long volIndex) const +{ return evaluateExpression(volIndex, BOUNDARY_YM_EXP); } -double VolumeVarContextExpression::getYpBoundaryFlux(long volIndex){ +double VolumeVarContextExpression::getYpBoundaryFlux(long volIndex) const +{ return evaluateExpression(volIndex, BOUNDARY_YP_EXP); } -double VolumeVarContextExpression::getZmBoundaryFlux(long volIndex){ +double VolumeVarContextExpression::getZmBoundaryFlux(long volIndex) const +{ return evaluateExpression(volIndex, BOUNDARY_ZM_EXP); } -double VolumeVarContextExpression::getZpBoundaryFlux(long volIndex){ +double VolumeVarContextExpression::getZpBoundaryFlux(long volIndex) const +{ return evaluateExpression(volIndex, BOUNDARY_ZP_EXP); } -double VolumeVarContextExpression::getXBoundaryPeriodicConstant() { +double VolumeVarContextExpression::getXBoundaryPeriodicConstant() const +{ return evaluateConstantExpression(BOUNDARY_XP_EXP); } -double VolumeVarContextExpression::getYBoundaryPeriodicConstant() { +double VolumeVarContextExpression::getYBoundaryPeriodicConstant() const +{ return evaluateConstantExpression(BOUNDARY_YP_EXP); } -double VolumeVarContextExpression::getZBoundaryPeriodicConstant() { +double VolumeVarContextExpression::getZBoundaryPeriodicConstant() const +{ return evaluateConstantExpression(BOUNDARY_ZP_EXP); } -double VolumeVarContextExpression::getConvectionVelocity_X(long volIndex){ +double VolumeVarContextExpression::getConvectionVelocity_X(long volIndex) const +{ return evaluateExpression(volIndex, VELOCITY_X_EXP); } -double VolumeVarContextExpression::getConvectionVelocity_Y(long volIndex){ +double VolumeVarContextExpression::getConvectionVelocity_Y(long volIndex) const +{ return evaluateExpression(volIndex, VELOCITY_Y_EXP); } -double VolumeVarContextExpression::getConvectionVelocity_Z(long volIndex){ +double VolumeVarContextExpression::getConvectionVelocity_Z(long volIndex) const +{ return evaluateExpression(volIndex, VELOCITY_Z_EXP); } -double VolumeVarContextExpression::getFlux(MembraneElement *element) { +double VolumeVarContextExpression::getFlux(MembraneElement *element) const +{ //return evaluateExpression(element, FLUX_EXP); return evaluateJumpCondition(element); } -bool VolumeVarContextExpression::isNullExpressionOK(int expIndex) { +bool VolumeVarContextExpression::isNullExpressionOK(int expIndex) const { if (expIndex == INITIAL_VALUE_EXP || expIndex == REACT_RATE_EXP) { return false; } Solver* solver = sim->getSolverFromVariable(species); - if (solver != NULL && solver->isPDESolver()) { + if (solver != nullptr && solver->isPDESolver()) { if (expIndex == DIFF_RATE_EXP) { return false; } @@ -132,14 +150,16 @@ bool VolumeVarContextExpression::isNullExpressionOK(int expIndex) { return true; } -bool VolumeVarContextExpression::hasConstantDiffusion() { +bool VolumeVarContextExpression::hasConstantDiffusion() const +{ if (!species->isDiffusing()) { - throw "hasConstantDiffusion() is only for PDE variables"; + throw std::runtime_error("hasConstantDiffusion() is only for PDE variables"); } return isConstantExpression(DIFF_RATE_EXP); } -bool VolumeVarContextExpression::hasConstantCoefficients(int dimension) { +bool VolumeVarContextExpression::hasConstantCoefficients(int dimension) const +{ if (!hasConstantDiffusion() || ((VolumeVariable*)species)->hasGradient()) { return false; } @@ -151,14 +171,16 @@ bool VolumeVarContextExpression::hasConstantCoefficients(int dimension) { && (dimension < 3 || isConstantExpression(VELOCITY_Z_EXP)); } -bool VolumeVarContextExpression::hasXYZOnlyDiffusion() { +bool VolumeVarContextExpression::hasXYZOnlyDiffusion() const +{ if (!species->isDiffusing()) { - throw "hasConstantDiffusion() is only for PDE variables"; + throw std::runtime_error("hasConstantDiffusion() is only for PDE variables"); } return isXYZOnlyExpression(DIFF_RATE_EXP); } -bool VolumeVarContextExpression::hasGradient(int dir) { +bool VolumeVarContextExpression::hasGradient(int dir) const +{ int expIndex = GRADIENT_X_EXP; if (dir == 1) { expIndex = GRADIENT_Y_EXP; diff --git a/VCell/tests/CMakeLists.txt b/VCell/tests/CMakeLists.txt new file mode 100644 index 000000000..bb6fa1d93 --- /dev/null +++ b/VCell/tests/CMakeLists.txt @@ -0,0 +1,33 @@ +cmake_minimum_required (VERSION 3.22) +project(VCellTest CXX) +set (CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +include(FetchContent) +FetchContent_Declare( + googletest + GIT_REPOSITORY https://github.com/google/googletest.git + GIT_TAG v1.14.0 +) + +# For Windows: Prevent overriding the parent project's compiler/linker settings +set(gtest_force_shared_crt ON CACHE BOOL "" FORCE) +FetchContent_MakeAvailable(googletest) + +if (NOT WINDOWS) + set(gtest_force_shared_crt ON CACHE BOOL "" FORCE) +endif() + +set(SRC_FILES + TestVCell.cpp +) + +enable_testing() + +add_executable(vcell_test ${SRC_FILES} ) +target_link_libraries(vcell_test vcell GTest::gtest_main) +#target_link_libraries(TestVCell vcell GTest::gtest_main) +#install(TARGETS TestVCell +# RUNTIME DESTINATION bin) + +include(GoogleTest) +gtest_discover_tests(vcell_test) \ No newline at end of file diff --git a/VCell/tests/TestVCell.cpp b/VCell/tests/TestVCell.cpp new file mode 100644 index 000000000..687af7160 --- /dev/null +++ b/VCell/tests/TestVCell.cpp @@ -0,0 +1,24 @@ +#include + +#include + +TEST(VCellTest, BasicAssertions) { + std::string inputDir = "../../testFiles/input/"; + std::string outputDir = "../../testFiles/output/"; + + std::string fvInputFile = inputDir + "SimID_11538992_0_.fvinput"; + std::string vcgInputFile = inputDir + "SimID_11538992_0_.vcg"; + std::string testOutputDir_1 = outputDir + "test_output_1"; + std::string testOutputDir_2 = outputDir + "test_output_2"; + + auto retcode_1 = solve(fvInputFile, vcgInputFile, testOutputDir_1); + EXPECT_EQ(retcode_1, 0); + + auto retcode_2 = solve(fvInputFile, vcgInputFile, testOutputDir_2); + EXPECT_EQ(retcode_2, 0); + } + +int main(int argc, char **argv) { + testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/VCellMessaging/src/SimulationMessaging.cpp b/VCellMessaging/src/SimulationMessaging.cpp index 55f241580..f93337da9 100644 --- a/VCellMessaging/src/SimulationMessaging.cpp +++ b/VCellMessaging/src/SimulationMessaging.cpp @@ -112,6 +112,7 @@ namespace { SimulationMessaging::~SimulationMessaging() throw() { if (workerEventOutputMode == WORKEREVENT_OUTPUT_MODE_STDOUT) { + m_inst = nullptr; return; } std::for_each(events.begin(), events.end( ),cleanupWorkerEvent); @@ -134,6 +135,7 @@ SimulationMessaging::~SimulationMessaging() throw() #endif curl_global_cleanup(); #endif + m_inst = nullptr; } SimulationMessaging* SimulationMessaging::getInstVar() { @@ -326,7 +328,7 @@ void SimulationMessaging::sendStatus() { #endif } #else - throw "OUPUT_MODE_STANDOUT must be using if not using messaging!"; + throw runtime_error("OUPUT_MODE_STANDOUT must be using if not using messaging!"); #endif } delete workerEvent; @@ -334,7 +336,7 @@ void SimulationMessaging::sendStatus() { void SimulationMessaging::setWorkerEvent(WorkerEvent* arg_workerEvent) { if (m_inst == nullptr) { - throw "SimulationMessaging is not initialized"; + throw runtime_error("SimulationMessaging is not initialized"); } if (workerEventOutputMode == WORKEREVENT_OUTPUT_MODE_STDOUT) { events.push_back(arg_workerEvent); diff --git a/VCellZipUtils/include/VCELL/ZipUtils.h b/VCellZipUtils/include/VCELL/ZipUtils.h index c98474253..a8e61d080 100644 --- a/VCellZipUtils/include/VCELL/ZipUtils.h +++ b/VCellZipUtils/include/VCELL/ZipUtils.h @@ -8,7 +8,9 @@ #ifndef ZIPUTILS_H #define ZIPUTILS_H -extern void extractFileFromZip(const char *zipFilename, const char *zipEntryName); -extern void addFilesToZip(const char *zipFilename, const char *filename1, const char *filename2=NULL); +#include + +extern void extractFileFromZip(const std::filesystem::path& zipFilename, const std::filesystem::path& zipEntryName); +extern void addFilesToZip(const std::filesystem::path& zipFilename, const std::filesystem::path& filename1, const std::filesystem::path& filename2=""); #endif diff --git a/VCellZipUtils/src/ZipUtils.cpp b/VCellZipUtils/src/ZipUtils.cpp index 9ce1b7d91..4476b858d 100644 --- a/VCellZipUtils/src/ZipUtils.cpp +++ b/VCellZipUtils/src/ZipUtils.cpp @@ -4,22 +4,17 @@ */ #include -#include -using std::stringstream; + using std::cout; using std::endl; #include #include -#include #include -#include +#include #include -#include #include -#include -#include -#include +#include #include #include @@ -45,7 +40,7 @@ ziptool -c test.zip add_file file4.txt ../VCellZipUtils/test/file2.txt 0 0 set_f */ static zip_t * -read_from_file(const char *archive, int flags, zip_error_t *error) +read_from_file(const filesystem::path& archive, int flags, zip_error_t *error) { zip_uint64_t offset = 0; @@ -56,9 +51,9 @@ read_from_file(const char *archive, int flags, zip_error_t *error) int err; if (offset == 0 && length == 0) { - if ((zaa = zip_open(archive, flags, &err)) == NULL) { + if ((zaa = zip_open(archive.string().c_str(), flags, &err)) == nullptr) { zip_error_set(error, err, errno); - return NULL; + return nullptr; } } else { @@ -66,10 +61,10 @@ read_from_file(const char *archive, int flags, zip_error_t *error) // zip_error_set(error, ZIP_ER_INVAL, 0); // return NULL; // } - if ((source = zip_source_file_create(archive, offset, (zip_int64_t)length, error)) == NULL - || (zaa = zip_open_from_source(source, flags, error)) == NULL) { + if ((source = zip_source_file_create(archive.string().c_str(), offset, (zip_int64_t)length, error)) == nullptr + || (zaa = zip_open_from_source(source, flags, error)) == nullptr) { zip_source_free(source); - return NULL; + return nullptr; } } @@ -80,7 +75,7 @@ read_from_file(const char *archive, int flags, zip_error_t *error) adds one or two files as uncompressed entries into a zip archive (creating the archive if necessary). entry names are the stripped filenames without full path. */ -void addFilesToZip(const char *ziparchive, const char *filepath1, const char *filepath2) +void addFilesToZip(const filesystem::path& ziparchive, const filesystem::path& filepath1, const filesystem::path& filepath2) { #if ( !defined(WIN32) && !defined(WIN64) && defined(USE_MESSAGING)) // UNIX #define USE_SHELL_ZIP 1 @@ -89,18 +84,13 @@ void addFilesToZip(const char *ziparchive, const char *filepath1, const char *fi #ifdef USE_SHELL_ZIP //printf("---------- using shell zip --------\n"); - char zipcommand [1024]; - if (filepath2 != NULL){ - - sprintf(zipcommand, "zip -j -Z store -g %s %s %s", ziparchive, filepath1, filepath2); - }else{ - sprintf(zipcommand, "zip -j -Z store -g %s %s", ziparchive, filepath1); + string zipcommand = string("zip -j -Z store -g ") + ziparchive.string() + " " + filepath1.string(); + if (!filepath2.string().empty()){ + zipcommand = zipcommand + " " + filepath2.string(); } - int retcode = system(zipcommand); + int retcode = system(zipcommand.c_str()); if (retcode != 0){ - char errmsg [2048]; - sprintf(errmsg, "zip command failed: %s", zipcommand); - throw errmsg; + throw runtime_error("zip command failed: " + zipcommand); } #else @@ -111,56 +101,53 @@ void addFilesToZip(const char *ziparchive, const char *filepath1, const char *fi int existingEntries = 0; zip_t *za = read_from_file(ziparchive, 0, &error); - if (za != NULL){ + if (za != nullptr){ existingEntries = zip_get_num_entries(za, 0); zip_close(za); } - char indexFirstAddedEntry [128]; - sprintf(indexFirstAddedEntry, "%d", existingEntries); // // strip filename from filepath1 // - std::string entryName1 = filepath1; + std::string entryName1 = filepath1.string(); std::size_t dirPos1 = entryName1.find_last_of("/\\"); if (dirPos1 > 0){ entryName1 = entryName1.substr(dirPos1+1, entryName1.length()); } - + string indexFirstAddedEntry = to_string(existingEntries); + const char* argv[50]; int argc = 0; argv[argc++] = "ziptool_main"; argv[argc++] = "-cn"; - argv[argc++] = ziparchive; + argv[argc++] = ziparchive.string().c_str(); argv[argc++] = "add_file"; argv[argc++] = entryName1.c_str(); - argv[argc++] = filepath1; + argv[argc++] = filepath1.string().c_str(); argv[argc++] = "0"; argv[argc++] = "0"; argv[argc++] = "set_file_compression"; - argv[argc++] = indexFirstAddedEntry; + argv[argc++] = indexFirstAddedEntry.c_str(); argv[argc++] = "store"; argv[argc++] = "none"; - if (filepath2 != NULL){ + if (!filepath2.string().empty()){ // // strip filename from filepath2 // - std::string entryName2 = filepath2; + std::string entryName2 = filepath2.string(); std::size_t dirPos2 = entryName2.find_last_of("/\\"); if (dirPos2 > 0){ entryName2 = entryName1.substr(dirPos2+1, entryName2.length()); } - - char indexSecondAddedEntry [128]; - sprintf(indexSecondAddedEntry, "%d", existingEntries + 1); + string indexSecondAddedEntry = to_string(existingEntries + 1); argv[argc++] = "add_file"; argv[argc++] = entryName2.c_str(); - argv[argc++] = filepath2; + argv[argc++] = filepath2.string().c_str(); argv[argc++] = "0"; argv[argc++] = "0"; argv[argc++] = "set_file_compression"; - argv[argc++] = indexSecondAddedEntry; + argv[argc++] = indexSecondAddedEntry.c_str(); argv[argc++] = "store"; argv[argc++] = "none"; } @@ -174,7 +161,7 @@ name_locate(zip_t *za, const char* entryName) { zip_int64_t idx; if ((idx=zip_name_locate(za, entryName, flags)) < 0) { - throw "can't find named entry"; + throw runtime_error(string("can't find named entry: ") + string(entryName)); } return idx; @@ -229,22 +216,21 @@ void extractFileFromZip(const char *zipFilename, const char *zipEntryName){ } */ -void extractFileFromZip(const char *zipFilename, const char *zipEntryName) +void extractFileFromZip(const std::filesystem::path& zipFilename, const std::filesystem::path& zipEntryName) { - char buf[100]; + char buf[1000]; int err; - int fd; - long long sum; - zip_t* za = NULL; + zip_t* za = nullptr; - if ((za = zip_open(zipFilename, 0, &err)) == NULL) { + if ((za = zip_open(zipFilename.string().c_str(), 0, &err)) == nullptr) { zip_error_to_str(buf, sizeof(buf), err, errno); - fprintf(stderr, "can't open zip archive `%s': %s\n", zipFilename, buf); - throw "cannot open zip archive"; + string errmsg = string("can't open zip archive ") + zipFilename.string() + " : " + buf; + cerr << errmsg << endl; + throw runtime_error(errmsg); } - zip_int64_t i = name_locate(za, zipEntryName); + zip_int64_t i = name_locate(za, zipEntryName.string().c_str()); zip_stat_t sb; if (zip_stat_index(za, i, 0, &sb) == 0) { @@ -253,22 +239,25 @@ void extractFileFromZip(const char *zipFilename, const char *zipEntryName) printf("mtime: [%u]\n", (unsigned int)sb.mtime); zip_file_t *zf = zip_fopen_index(za, i, 0); if (!zf) { - fprintf(stderr, "failed to open zip archive %s\n", zipFilename); - throw "failed to open archive"; + string errmsg = string("failed to open zip archive ") + zipFilename.string(); + std::cerr << errmsg << std::endl; + throw runtime_error(errmsg); } - fd = open(sb.name, O_RDWR | O_TRUNC | O_CREAT, 0644); + int fd = open(sb.name, O_RDWR | O_TRUNC | O_CREAT, 0644); if (fd < 0) { - fprintf(stderr, "failed to open zip entry %s : %s\n",zipFilename,zipEntryName); - throw "failed to open zip entry"; + string errmsg = string("failed to open zip entry ") + zipFilename.string() + " : " + zipEntryName.string(); + cerr << errmsg << endl; + throw runtime_error(errmsg); } - sum = 0; + long long sum = 0; while (sum != sb.size) { int len = zip_fread(zf, buf, 100); if (len < 0) { - fprintf(stderr, "failed to read zip entry %s : %s\n", zipFilename, zipEntryName); - throw "failed to read zip entry"; + string errmsg = string("failed to read zip entry ") + zipFilename.string() + " : " + zipEntryName.string(); + cerr << errmsg << endl; + throw runtime_error(errmsg); } write(fd, buf, len); sum += len; @@ -281,8 +270,9 @@ void extractFileFromZip(const char *zipFilename, const char *zipEntryName) if (zip_close(za) == -1) { - fprintf(stderr, "%can't close zip archive `%s'\n", zipFilename); - throw "failed to close zip archive"; + string errMsg = string("can't close zip archive ") + zipFilename.string(); + cerr << errMsg << endl; + throw runtime_error(errMsg); } } diff --git a/bridgeVCellSmoldyn/VCellMesh.cpp b/bridgeVCellSmoldyn/VCellMesh.cpp index a43a22de4..8a70dc747 100644 --- a/bridgeVCellSmoldyn/VCellMesh.cpp +++ b/bridgeVCellSmoldyn/VCellMesh.cpp @@ -1,16 +1,16 @@ -#include "VCellMesh.h" -#include -#include -#include - -using namespace std; - - -VCellMesh::VCellMesh(SimTool* simTool) -{ - this->simTool = simTool; -} - +#include "VCellMesh.h" +#include +#include +#include + +using namespace std; + + +VCellMesh::VCellMesh(SimTool* simTool) +{ + this->simTool = simTool; +} + void VCellMesh::getCenterCoordinates(int volIndex, double* coords) { CartesianMesh* mesh = ((CartesianMesh*)simTool->getSimulation()->getMesh()); diff --git a/bridgeVCellSmoldyn/VCellValueProvider.cpp b/bridgeVCellSmoldyn/VCellValueProvider.cpp index 383002245..2a8575554 100644 --- a/bridgeVCellSmoldyn/VCellValueProvider.cpp +++ b/bridgeVCellSmoldyn/VCellValueProvider.cpp @@ -19,7 +19,7 @@ double VCellValueProvider::getConstantValue() { } double VCellValueProvider::getValue(double t, double x, double y, double z, rxnptr rxn) { - SimulationExpression* vcellSim = (SimulationExpression*)simTool->getSimulation(); + SimulationExpression* vcellSim = simTool->getSimulation(); WorldCoord wc(x, y, z); vcellSim->setCurrentCoordinate(wc); @@ -35,7 +35,7 @@ double VCellValueProvider::getValue(double t, double x, double y, double z, rxnp double VCellValueProvider::getValue(double t, double x, double y, double z, rxnptr rxn, char* panelName) { - SimulationExpression* vcellSim = (SimulationExpression*)simTool->getSimulation(); + SimulationExpression* vcellSim = simTool->getSimulation(); int* indices = vcellSim->getIndices(); WorldCoord wc(x, y, z); vcellSim->setCurrentCoordinate(wc); @@ -60,7 +60,7 @@ double VCellValueProvider::getValue(double t, double x, double y, double z, rxnp double VCellValueProvider::getValue(double t, double x, double y, double z, surfactionptr actiondetails, char* panelName){ - SimulationExpression* vcellSim = (SimulationExpression*)simTool->getSimulation(); + SimulationExpression* vcellSim = simTool->getSimulation(); int* indices = vcellSim->getIndices(); WorldCoord wc(x, y, z); vcellSim->setCurrentCoordinate(wc); diff --git a/bridgeVCellSmoldyn/vcellhybrid.cpp b/bridgeVCellSmoldyn/vcellhybrid.cpp index 2ad8c8410..968cdfb22 100644 --- a/bridgeVCellSmoldyn/vcellhybrid.cpp +++ b/bridgeVCellSmoldyn/vcellhybrid.cpp @@ -48,7 +48,7 @@ simptr vcellhybrid::smoldynInit(SimTool* simTool, string& fileName) { sim->clockstt=time(NULL); er=simdocommands(sim); - SimulationExpression* vcellSim = (SimulationExpression*)simTool->getSimulation(); + SimulationExpression* vcellSim = simTool->getSimulation(); SymbolTable* symbolTable = vcellSim->getSymbolTable(); char erstr[1024]; //initialization for reaction rates (as expression) diff --git a/pyproject.toml b/pyproject.toml index 97930b334..308c7b9de 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,19 +5,17 @@ build-backend = "scikit_build_core.build" [project] name = "pyvcell_fvsolver" -version = "0.0.1" +version = "0.0.4" description="Virtual Cell Finite Volume PDE solver" readme = "README.md" authors = [ { name = "Jim Schaff", email = "schaff@uchc.edu" }, ] -requires-python = ">=3.7" +requires-python = ">=3.9" classifiers = [ "Development Status :: 4 - Beta", "License :: OSI Approved :: MIT License", "Programming Language :: Python :: 3 :: Only", - "Programming Language :: Python :: 3.7", - "Programming Language :: Python :: 3.8", "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", @@ -48,6 +46,7 @@ test-command = "pytest {project}/tests" test-extras = ["test"] test-skip = ["*universal2:arm64"] build-verbosity = 1 +archs = ["native"] [tool.ruff] diff --git a/src/main.cpp b/src/main.cpp index 77073fbcc..31dbac271 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,6 +1,6 @@ #include -#include "SolverMain.h" +#include #define STRINGIFY(x) #x #define MACRO_STRINGIFY(x) STRINGIFY(x) @@ -32,7 +32,7 @@ PYBIND11_MODULE(_core, m) { The inputFilename expects a .fvinput file, the outputDir will be created as needed. )pbdoc", - py::arg("inputFilename"), py::arg("outputDir")); + py::arg("fvInputFilename"), py::arg("vcgInputFilename"), py::arg("outputDir")); #ifdef VERSION_INFO m.attr("__version__") = MACRO_STRINGIFY(VERSION_INFO); diff --git a/tests/test_basic.py b/tests/test_basic.py index 1077a16ea..ba029cb28 100644 --- a/tests/test_basic.py +++ b/tests/test_basic.py @@ -1,7 +1,9 @@ import pyvcell_fvsolver as fv -def test_version(): - assert fv.__version__ == "0.0.1" def test_version(): + assert fv.__version__ == "0.0.4" + + +def test_version_function(): assert fv.version() is not None