From 948538ac08c402358baa086eefd7925bd6f6a800 Mon Sep 17 00:00:00 2001 From: Daniel-H2 <167224060+Daniel-H2@users.noreply.github.com> Date: Wed, 17 Apr 2024 16:05:18 -0400 Subject: [PATCH 1/3] Update README.rst --- README.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.rst b/README.rst index e016ba418..05cfbfd7f 100644 --- a/README.rst +++ b/README.rst @@ -1,4 +1,4 @@ -CADET +CADET (my version) ====== .. image:: https://img.shields.io/github/release/modsim/cadet.svg From b0d4b6a8934c657b7d3d069104407a5e8746673e Mon Sep 17 00:00:00 2001 From: Daniel-H2 Date: Wed, 17 Apr 2024 16:26:17 -0400 Subject: [PATCH 2/3] Custom binding model Test Following the tutorial for custom binding model and trying to add the custom binding model Test. --- .github/workflows/ci.yml | 186 ----------------- .github/workflows/deploy_documentation.yml | 16 -- src/libcadet/BindingModelFactory.cpp | 2 + src/libcadet/CMakeLists.txt | 1 + src/libcadet/model/binding/TestBinding.cpp | 224 +++++++++++++++++++++ 5 files changed, 227 insertions(+), 202 deletions(-) delete mode 100644 .github/workflows/ci.yml delete mode 100644 .github/workflows/deploy_documentation.yml create mode 100644 src/libcadet/model/binding/TestBinding.cpp diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml deleted file mode 100644 index d08355264..000000000 --- a/.github/workflows/ci.yml +++ /dev/null @@ -1,186 +0,0 @@ -name: CI - -on: - push: - branches: - - master - - dev - pull_request: - - -jobs: - Win64: - runs-on: windows-latest - strategy: - fail-fast: true - env: - BASE_DIR: ${{ github.workspace }} - SRC_DIR: ${{ github.workspace }}/src - BUILD_DIR: ${{ github.workspace }}/build - INSTALL_PREFIX: ${{ github.workspace }}/cadet - steps: - - uses: actions/checkout@v4 - with: - submodules: true - path: src - - uses: ilammy/msvc-dev-cmd@v1 - - name: Get MSVC version - id: get-msvc - run: | - $MSVC_VERSION = cl 2>&1 | Select-String -Pattern 'Version (\d\d\.\d\d)' | % { $_.matches.groups[1].Value } - echo "version=$MSVC_VERSION" >> $Env:GITHUB_OUTPUT - - uses: actions/cache@v4 - id: cache - with: - path: | - ${{ github.workspace }}/hdf5 - ${{ github.workspace }}/suitesparse/install - key: ${{ runner.os }}-deps-${{ steps.get-msvc.outputs.version }} - - name: Download MKL and TBB - run: | - cd "${env:BASE_DIR}" - $base_dir = $($env:BASE_DIR.Replace('\', '/')) - nuget install inteltbb.devel.win -Version 2021.8.0.25874 - nuget install intelmkl.static.win-x64 -Version 2023.0.0.25930 - Invoke-WebRequest -Uri "https://gitlab.com/libeigen/eigen/-/archive/master/eigen-master.zip" -OutFile eigen.zip - 7z x eigen.zip - - name: Build UMFPACK and HDF5 - if: steps.cache.outputs.cache-hit != 'true' - run: | - $base_dir = $($env:BASE_DIR.Replace('\', '/')) - Invoke-WebRequest -Uri "https://www.hdfgroup.org/package/cmake-hdf5-1-14-0-zip/?wpdmdl=17553" -OutFile hdf5.zip - 7z x hdf5.zip - cd CMake-hdf5-1.14.0 - cmake -DCMAKE_BUILD_TYPE:STRING=Release -DBUILD_SHARED_LIBS:BOOL=ON -DHDF5_BUILD_FORTRAN:BOOL=OFF -DHDF5_ENABLE_F2003:BOOL=OFF -DHDF5_BUILD_JAVA:BOOL=OFF -DCMAKE_INSTALL_PREFIX:PATH="${base_dir}/hdf5" -DCTEST_CONFIGURATION_TYPE:STRING=Release -DBUILD_TESTING=ON -DHDF5_BUILD_TOOLS=OFF -DHDF5_BUILD_EXAMPLES=OFF -DHDF5_BUILD_HL_LIB=OFF -DHDF5_BUILD_CPP_LIB=OFF -DHDF5_ALLOW_EXTERNAL_SUPPORT:STRING=TGZ -DTGZPATH:PATH="${base_dir}/CMake-hdf5-1.14.0" -DHDF5_PACKAGE_EXTLIBS:BOOL=ON -DSITE:STRING=WIN10VS202264.XXXX -DBUILDNAME:STRING=Windows-WIN10-vs2022-STATIC -G "Ninja" hdf5-1.14.0 - ninja install - cd "${env:BASE_DIR}" - Invoke-WebRequest -Uri "https://github.com/DrTimothyAldenDavis/SuiteSparse/archive/v5.7.2.zip" -OutFile suitesparse.zip - 7z x suitesparse.zip -osuitesparse\code -y - Invoke-WebRequest -Uri "https://github.com/jlblancoc/suitesparse-metis-for-windows/archive/e8d953dffb8a99aa8b65ff3ff03e12a3ed72f90c.zip" -OutFile ssb.zip - 7z x ssb.zip -osuitesparse\code -y - cd ${env:BASE_DIR}\suitesparse\code - xcopy .\SuiteSparse-5.7.2 .\suitesparse-metis-for-windows-e8d953dffb8a99aa8b65ff3ff03e12a3ed72f90c\SuiteSparse /s /e /y /q - cd "${env:BASE_DIR}\suitesparse\code\suitesparse-metis-for-windows-e8d953dffb8a99aa8b65ff3ff03e12a3ed72f90c" - Remove-Item -ErrorAction Ignore -Recurse -Force lapack_windows - Remove-Item -ErrorAction Ignore -Recurse -Force lapack_windows - cd "${env:BASE_DIR}\suitesparse\code" - mkdir build - cd build - $ENV:MKLROOT="${env:BASE_DIR}/intelmkl.static.win-x64.2023.0.0.25930/lib/native/win-x64".Replace('\', '/') - cmake -DCMAKE_INSTALL_PREFIX="${base_dir}\suitesparse\install" -DBLA_VENDOR=Intel10_64lp_seq -DBLA_STATIC=ON -G "Ninja" -DCMAKE_C_FLAGS="/GL" -DCMAKE_STATIC_LINKER_FLAGS="/LTCG" -DCMAKE_BUILD_TYPE=Release -DBUILD_METIS=OFF ..\suitesparse-metis-for-windows-e8d953dffb8a99aa8b65ff3ff03e12a3ed72f90c\ - ninja install - - name: Build and Install - run: | - cd "${env:BASE_DIR}" - cmake -E make_directory "${env:BUILD_DIR}" - cmake -E make_directory "${env:INSTALL_PREFIX}" - cd "${env:BUILD_DIR}" - $ENV:MKLROOT="${env:BASE_DIR}/intelmkl.static.win-x64.2023.0.0.25930/lib/native/win-x64".Replace('\', '/') - $ENV:TBB_ROOT="${env:BASE_DIR}/inteltbb.devel.win.2021.8.0.25874/lib/native".Replace('\', '/') - $ENV:UMFPACK_ROOT="${env:BASE_DIR}/suitesparse/install".Replace('\', '/') - $install_prefix = $($env:INSTALL_PREFIX.Replace('\', '/')) - $src_dir = $($env:SRC_DIR.Replace('\', '/')) - $base_dir = $($env:BASE_DIR.Replace('\', '/')) - cmake -G "Ninja" -DCMAKE_INSTALL_PREFIX="${install_prefix}" -DCMAKE_BUILD_TYPE=Release -DHDF5_ROOT="${base_dir}/hdf5" -DENABLE_STATIC_LINK_DEPS=ON -DENABLE_STATIC_LINK_LAPACK=ON -DBLA_VENDOR=Intel10_64lp_seq "${src_dir}" - ninja - ninja install - - name: Include Deps and Package - run: | - cd "${env:BASE_DIR}" - Copy-Item "${env:BASE_DIR}\inteltbb.redist.win.2021.8.0.25874\runtimes\win-x64\native\tbb12.dll" "${env:INSTALL_PREFIX}\bin\" - Copy-Item "${env:BASE_DIR}\inteltbb.redist.win.2021.8.0.25874\runtimes\win-x64\native\tbb12.dll" "${env:BUILD_DIR}\src\cadet-cli\" - Copy-Item "${env:BASE_DIR}\inteltbb.redist.win.2021.8.0.25874\runtimes\win-x64\native\tbb12.dll" "${env:BUILD_DIR}\test\" - 7z a "${env:BASE_DIR}\cadet-${{ runner.os }}.zip" cadet\ - - name: Check if it runs - run: | - cd "${env:INSTALL_PREFIX}\bin" - .\cadet-cli.exe --version - .\createLWE.exe - .\cadet-cli.exe LWE.h5 - - name: Upload Artifacts - uses: actions/upload-artifact@v4 - with: - name: cadet-${{ runner.os }} - retention-days: 5 - path: ${{ github.workspace }}/cadet-${{ runner.os }}.zip - Ubuntu-latest: - runs-on: ubuntu-latest - strategy: - fail-fast: true - defaults: - run: - shell: bash -l {0} - env: - SRC_DIR: ${{ github.workspace }}/src - BUILD_DIR: ${{ github.workspace }}/build - INSTALL_PREFIX: ${{ github.workspace }}/install - BUILD_TYPE: Release - steps: - - uses: actions/checkout@v4 - with: - submodules: true - path: src - - name: Install Dependencies - run: | - sudo apt-get update - sudo apt -y install \ - build-essential \ - libhdf5-dev \ - liblapack-dev \ - libblas-dev \ - libtbb-dev \ - libsuperlu-dev \ - libeigen3-dev; - - name: Build and Install - run: | - cmake -E make_directory "${BUILD_DIR}" - cmake -E make_directory "${INSTALL_PREFIX}" - cd "${BUILD_DIR}" - cmake -DCMAKE_INSTALL_PREFIX="${INSTALL_PREFIX}" -DCMAKE_BUILD_TYPE="${BUILD_TYPE}" "${SRC_DIR}" - make install -j$(nproc) - - name: Check if it runs - run: | - export LD_LIBRARY_PATH=${INSTALL_PREFIX}/lib:$LD_LIBRARY_PATH - ${INSTALL_PREFIX}/bin/cadet-cli --version || true - ${INSTALL_PREFIX}/bin/createLWE - ${INSTALL_PREFIX}/bin/cadet-cli LWE.h5 || true - - name: Run tests - run: | - ${BUILD_DIR}/test/testRunner [CI] - MacOS: - runs-on: macos-latest - strategy: - fail-fast: true - defaults: - run: - shell: bash -l {0} - env: - SRC_DIR: ${{ github.workspace }}/src - BUILD_DIR: ${{ github.workspace }}/build - INSTALL_PREFIX: ${{ github.workspace }}/install - BUILD_TYPE: Release - steps: - - uses: actions/checkout@v4 - with: - submodules: true - path: src - - name: Install Dependencies - run: | - brew update > /dev/null || true - brew install cmake --without-docs - brew install hdf5 - brew install tbb - brew install eigen - - name: Build and Install - run: | - cmake -E make_directory "${BUILD_DIR}" - cmake -E make_directory "${INSTALL_PREFIX}" - cd "${BUILD_DIR}" - cmake -DCMAKE_INSTALL_PREFIX="${INSTALL_PREFIX}" -DCMAKE_BUILD_TYPE="${BUILD_TYPE}" "${SRC_DIR}" - make install -j$(sysctl -n hw.logicalcpu) - - name: Check if it runs - run: | - ${INSTALL_PREFIX}/bin/cadet-cli --version || true - ${INSTALL_PREFIX}/bin/createLWE - ${INSTALL_PREFIX}/bin/cadet-cli LWE.h5 || true - diff --git a/.github/workflows/deploy_documentation.yml b/.github/workflows/deploy_documentation.yml deleted file mode 100644 index 662e7f1e2..000000000 --- a/.github/workflows/deploy_documentation.yml +++ /dev/null @@ -1,16 +0,0 @@ -name: "Trigger for dispatching CADET-Website" -on: - push: - branches: - - 'master' - -jobs: - dispatch_docs: - runs-on: ubuntu-latest - steps: - - name: Repository Dispatch - uses: peter-evans/repository-dispatch@v2 - with: - token: ${{ secrets.DISPATCH_DOCS }} - repository: modsim/CADET-Website - event-type: build_docs diff --git a/src/libcadet/BindingModelFactory.cpp b/src/libcadet/BindingModelFactory.cpp index b937462ea..3105f2091 100644 --- a/src/libcadet/BindingModelFactory.cpp +++ b/src/libcadet/BindingModelFactory.cpp @@ -44,6 +44,7 @@ namespace cadet void registerBiLangmuirLDFModel(std::unordered_map>& bindings); void registerHICWaterOnHydrophobicSurfacesModel(std::unordered_map>& bindings); void registerHICConstantWaterActivityModel(std::unordered_map>& bindings); + void registerTestModel(std::unordered_map>& bindings); } } @@ -73,6 +74,7 @@ namespace cadet model::binding::registerBiLangmuirLDFModel(_bindingModels); model::binding::registerHICWaterOnHydrophobicSurfacesModel(_bindingModels); model::binding::registerHICConstantWaterActivityModel(_bindingModels); + model::binding::registerTestModel(_bindingModels); registerModel(); } diff --git a/src/libcadet/CMakeLists.txt b/src/libcadet/CMakeLists.txt index c1c550a0b..4e93d42f0 100644 --- a/src/libcadet/CMakeLists.txt +++ b/src/libcadet/CMakeLists.txt @@ -100,6 +100,7 @@ set(LIBCADET_BINDINGMODEL_SOURCES ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/BiLangmuirLDFBinding.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/HICWaterOnHydrophobicSurfacesBinding.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/HICConstantWaterActivityBinding.cpp + ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/TestBinding.cpp ) # LIBCADET_REACTIONMODEL_SOURCES holds all source files of reaction models diff --git a/src/libcadet/model/binding/TestBinding.cpp b/src/libcadet/model/binding/TestBinding.cpp new file mode 100644 index 000000000..3e5bea7fc --- /dev/null +++ b/src/libcadet/model/binding/TestBinding.cpp @@ -0,0 +1,224 @@ +// ============================================================================= +// CADET +// +// Copyright © 2008-2021: The CADET Authors +// Please see the AUTHORS and CONTRIBUTORS file. +// +// All rights reserved. This program and the accompanying materials +// are made available under the terms of the GNU Public License v3.0 (or, at +// your option, any later version) which accompanies this distribution, and +// is available at http://www.gnu.org/licenses/gpl.html +// ============================================================================= + +#include "model/binding/BindingModelBase.hpp" +#include "model/ExternalFunctionSupport.hpp" +#include "model/ModelUtils.hpp" +#include "cadet/Exceptions.hpp" +#include "model/Parameters.hpp" +#include "LocalVector.hpp" +#include "SimulationTypes.hpp" + +#include +#include +#include +#include + +/* In the following lines of code the interface to the binding model is defined. + In the interface the binding model specific parameters are defined. Please modify configuration name and + parameter names according to the naming convection explained in the tutorial.*/ + +/* +{ + "name": TestParamHandler", + "externalName": "ExtTestParamHandler", + "parameters": + [ + { "type": "ScalarComponentDependentParameter", "varName": "kA", "confName": "MCL_KA"}, + { "type": "ScalarComponentDependentParameter", "varName": "kD", "confName": "MCL_KD"}, + { "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "MCL_QMAX"} + ] +} +*/ + +/* Parameter description + ------------------------ + kA = Adsorption rate + kD = Desorption rate + qMax = Capacity +*/ + +namespace cadet +{ + + namespace model + { + + inline const char* TestParamHandler::identifier() CADET_NOEXCEPT { return "TEST_MULTI_COMPONENT_LANGMUIR"; } + + inline bool TestParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) + { + if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("MCL_KA, MCL_KD, and MCL_QMAX have to have the same size"); + + return true; + } + + inline const char* ExtTestParamHandler::identifier() CADET_NOEXCEPT { return "EXT_TEST_MULTI_COMPONENT_LANGMUIR"; } + + inline bool ExtTestParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) + { + if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("EXT_MCL_KA, EXT_MCL_KD, and EXT_MCL_QMAX have to have the same size"); + + return true; + } + + + /** + * @brief Defines the multi component Langmuir binding model + * @details Implements the Langmuir adsorption model: \f[ \begin{align} + * \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} q_{\text{max},i} \left( 1 - \sum_j \frac{q_j}{q_{\text{max},j}} \right) - k_{d,i} q_i + * \end{align} \f] + * Multiple bound states are not supported. + * Components without bound state (i.e., non-binding components) are supported. + * + * See @cite Langmuir1916. + * @tparam ParamHandler_t Type that can add support for external function dependence + */ + template + class TestBindingBase : public ParamHandlerBindingModelBase + { + public: + + TestBindingBase() { } + virtual ~TestBindingBase() CADET_NOEXCEPT { } + + static const char* identifier() { return ParamHandler_t::identifier(); } + + virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; } + + CADET_BINDINGMODELBASE_BOILERPLATE + + protected: + using ParamHandlerBindingModelBase::_paramHandler; + using ParamHandlerBindingModelBase::_reactionQuasistationarity; + using ParamHandlerBindingModelBase::_nComp; + using ParamHandlerBindingModelBase::_nBoundStates; + + // In the follwing the class method the binding model mass transfer behavior is implemented. + // Apart from the mentioned places in the tutorial, do not modify anything in this function. + template + int fluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* y, + CpStateType const* yCp, ResidualType* res, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + // Protein fluxes: -k_{a,i} * c_{p,i} * q_{max,i} * (1 - \sum_j q_j / q_{max,j}) + k_{d,i} * q_i + ResidualType qSum = 1.0; + unsigned int bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + qSum -= y[bndIdx] / static_cast(p->qMax[i]); + + // Next bound component + ++bndIdx; + } + + bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + // Residual + res[bndIdx] = static_cast(p->kD[i]) * y[bndIdx] - static_cast(p->kA[i]) * yCp[i] * static_cast(p->qMax[i]) * qSum; + + // Next bound component + ++bndIdx; + } + + return 0; + } + // In the follwing the class method the anlytical Jacobian of the binding model should be implemented. + // Apart from the mentioned places in the tutorial, do not modify anything in this function. + template + void jacobianImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double const* yCp, int offsetCp, RowIterator jac, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + // Protein fluxes: -k_{a,i} * c_{p,i} * q_{max,i} * (1 - \sum_j q_j / q_{max,j}) + k_{d,i} * q_i + double qSum = 1.0; + int bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + qSum -= y[bndIdx] / static_cast(p->qMax[i]); + + // Next bound component + ++bndIdx; + } + + bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + const double ka = static_cast(p->kA[i]); + const double kd = static_cast(p->kD[i]); + + // dres_i / dc_{p,i} + jac[i - bndIdx - offsetCp] = -static_cast(p->kA[i]) * static_cast(p->qMax[i]) * qSum; + // Getting to c_{p,i}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0} and a +i to c_{p,i}. + // This means jac[i - bndIdx - offsetCp] corresponds to c_{p,i}. + + // Fill dres_i / dq_j + int bndIdx2 = 0; + for (int j = 0; j < _nComp; ++j) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[j] == 0) + continue; + + // dres_i / dq_j + jac[bndIdx2 - bndIdx] = static_cast(p->kA[i]) * yCp[i] * static_cast(p->qMax[i]) / static_cast(p->qMax[j]); + // Getting to q_j: -bndIdx takes us to q_0, another +bndIdx2 to q_j. This means jac[bndIdx2 - bndIdx] corresponds to q_j. + + ++bndIdx2; + } + + // Add to dres_i / dq_i + jac[0] += kd; + + // Advance to next flux and Jacobian row + ++bndIdx; + ++jac; + } + } + }; + + // Do not forget to make changes in the following lines of code according to the guidelines given in the tutorial. + typedef TestBindingBase TestBinding; + typedef TestBindingBase ExternalTestBinding; + + namespace binding + { + void registerTestModel(std::unordered_map>& bindings) + { + bindings[TestBinding::identifier()] = []() { return new TestBinding(); }; + bindings[ExternalTestBinding::identifier()] = []() { return new ExternalTestBinding(); }; + } + } // namespace binding + + } // namespace model + +} // namespace cadet From 1c33612d221c6b274a4eeeaf6c1f7cfe3fd78836 Mon Sep 17 00:00:00 2001 From: Daniel-H2 Date: Tue, 14 May 2024 22:36:30 -0400 Subject: [PATCH 3/3] ANN custom binding model --- src/libcadet/BindingModelFactory.cpp | 4 + src/libcadet/CMakeLists.txt | 2 + .../binding/LangmuirHybridDesorbBinding.cpp | 354 ++++++++++++++++++ src/libcadet/model/binding/TestBinding.cpp | 101 ++++- .../binding/TestStericMassActionBinding.cpp | 334 +++++++++++++++++ 5 files changed, 788 insertions(+), 7 deletions(-) create mode 100644 src/libcadet/model/binding/LangmuirHybridDesorbBinding.cpp create mode 100644 src/libcadet/model/binding/TestStericMassActionBinding.cpp diff --git a/src/libcadet/BindingModelFactory.cpp b/src/libcadet/BindingModelFactory.cpp index 3105f2091..0f6566764 100644 --- a/src/libcadet/BindingModelFactory.cpp +++ b/src/libcadet/BindingModelFactory.cpp @@ -45,6 +45,8 @@ namespace cadet void registerHICWaterOnHydrophobicSurfacesModel(std::unordered_map>& bindings); void registerHICConstantWaterActivityModel(std::unordered_map>& bindings); void registerTestModel(std::unordered_map>& bindings); + void registerTestStericMassActionModel(std::unordered_map>& bindings); + void registerLangmuirHybridDesorbModel(std::unordered_map>& bindings); } } @@ -75,6 +77,8 @@ namespace cadet model::binding::registerHICWaterOnHydrophobicSurfacesModel(_bindingModels); model::binding::registerHICConstantWaterActivityModel(_bindingModels); model::binding::registerTestModel(_bindingModels); + model::binding::registerTestStericMassActionModel(_bindingModels); + model::binding::registerLangmuirHybridDesorbModel(_bindingModels); registerModel(); } diff --git a/src/libcadet/CMakeLists.txt b/src/libcadet/CMakeLists.txt index 4e93d42f0..b5cae81b7 100644 --- a/src/libcadet/CMakeLists.txt +++ b/src/libcadet/CMakeLists.txt @@ -101,6 +101,8 @@ set(LIBCADET_BINDINGMODEL_SOURCES ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/HICWaterOnHydrophobicSurfacesBinding.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/HICConstantWaterActivityBinding.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/TestBinding.cpp + ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/TestStericMassActionBinding.cpp + ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/LangmuirHybridDesorbBinding.cpp ) # LIBCADET_REACTIONMODEL_SOURCES holds all source files of reaction models diff --git a/src/libcadet/model/binding/LangmuirHybridDesorbBinding.cpp b/src/libcadet/model/binding/LangmuirHybridDesorbBinding.cpp new file mode 100644 index 000000000..ac7eb9d84 --- /dev/null +++ b/src/libcadet/model/binding/LangmuirHybridDesorbBinding.cpp @@ -0,0 +1,354 @@ +// ============================================================================= +// CADET +// +// Copyright © 2008-2024: The CADET Authors +// Please see the AUTHORS and CONTRIBUTORS file. +// +// All rights reserved. This program and the accompanying materials +// are made available under the terms of the GNU Public License v3.0 (or, at +// your option, any later version) which accompanies this distribution, and +// is available at http://www.gnu.org/licenses/gpl.html +// ============================================================================= + +#include "model/binding/BindingModelBase.hpp" +#include "model/ExternalFunctionSupport.hpp" +#include "model/ModelUtils.hpp" +#include "cadet/Exceptions.hpp" +#include "model/Parameters.hpp" +#include "LocalVector.hpp" +#include "SimulationTypes.hpp" + +#include +#include +#include +#include + +/* +{ + "name": "LangmuirHDParamHandler", + "externalName": "ExtLangmuirHDParamHandler", + "parameters": + [ + { "type": "ScalarComponentDependentParameter", "varName": "kA", "confName": "MCL_KA"}, + { "type": "ScalarComponentDependentParameter", "varName": "kD", "confName": "MCL_KD"}, + { "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "MCL_QMAX"}, + { "type": "ScalarComponentDependentParameter", "varName": "p1", "confName": "ANN_P1"}, + { "type": "ScalarComponentDependentParameter", "varName": "p2", "confName": "ANN_P2"}, + { "type": "ScalarComponentDependentParameter", "varName": "p3", "confName": "ANN_P3"}, + { "type": "ScalarComponentDependentParameter", "varName": "p4", "confName": "ANN_P4"}, + { "type": "ScalarComponentDependentParameter", "varName": "p5", "confName": "ANN_P5"}, + { "type": "ScalarComponentDependentParameter", "varName": "p6", "confName": "ANN_P6"}, + { "type": "ScalarComponentDependentParameter", "varName": "p7", "confName": "ANN_P7"}, + { "type": "ScalarComponentDependentParameter", "varName": "p8", "confName": "ANN_P8"}, + { "type": "ScalarComponentDependentParameter", "varName": "p9", "confName": "ANN_P9"}, + { "type": "ScalarComponentDependentParameter", "varName": "p10", "confName": "ANN_P10"}, + { "type": "ScalarComponentDependentParameter", "varName": "p11", "confName": "ANN_P11"}, + { "type": "ScalarComponentDependentParameter", "varName": "P_norm", "confName": "ANN_P_NORM"}, + { "type": "ScalarComponentDependentParameter", "varName": "q_norm", "confName": "ANN_Q_NORM"}, + { "type": "ScalarComponentDependentParameter", "varName": "S_norm", "confName": "ANN_S_NORM"} + ] +} +*/ + +/* Parameter description + ------------------------ + kA = Adsorption rate + kD = Desorption rate + qMax = Capacity + p1 = ANN param 1 + p2 = ANN param 2 + p3 = ANN param 3 + p4 = ANN param 4 + p5 = ANN param 5 + p6 = ANN param 6 + p7 = ANN param 7 + p8 = ANN param 8 + p9 = ANN param 9 + p10 = ANN param 10 + p11 = ANN param 11 + P_norm = Normative Vaule for Protien pore conc + q_norm = Normative Vaule for Protien adsorbed conc + S_norm = Normative Vaule for Salt pore conc +*/ + +namespace cadet +{ + +namespace model +{ + +inline const char* LangmuirHDParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_LANGMUIR_HD"; } + +inline bool LangmuirHDParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +{ + if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _p1.size()) || (_kA.size() != _p2.size()) + || (_kA.size() != _p3.size()) || (_kA.size() != _p4.size()) || (_kA.size() != _p5.size()) || (_kA.size() != _p6.size()) + || (_kA.size() != _p7.size()) || (_kA.size() != _p8.size()) || (_kA.size() != _p9.size()) || (_kA.size() != _p10.size()) + || (_kA.size() != _p11.size()) || (_kA.size() != _P_norm.size()) || (_kA.size() != _q_norm.size()) || (_kA.size() != _S_norm.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("MCL_KA, MCL_KD, MCL_QMAX, and a ANN parameter have to have the same size"); + + return true; +} + +inline const char* ExtLangmuirHDParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_LANGMUIR_HD"; } + +inline bool ExtLangmuirHDParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +{ + if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _p1.size()) || (_kA.size() != _p2.size()) + || (_kA.size() != _p3.size()) || (_kA.size() != _p4.size()) || (_kA.size() != _p5.size()) || (_kA.size() != _p6.size()) + || (_kA.size() != _p7.size()) || (_kA.size() != _p8.size()) || (_kA.size() != _p9.size()) || (_kA.size() != _p10.size()) + || (_kA.size() != _p11.size()) || (_kA.size() != _P_norm.size()) || (_kA.size() != _q_norm.size()) || (_kA.size() != _S_norm.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("MCL_KA, MCL_KD, MCL_QMAX, and a ANN parameter have to have the same size"); + + return true; +} + + +/** + * @brief Defines the multi component LangmuirHD binding model + * @details Implements the LangmuirHD adsorption model: \f[ \begin{align} + * \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} q_{\text{max},i} \left( 1 - \sum_j \frac{q_j}{q_{\text{max},j}} \right) - k_{d,i} q_i + * \end{align} \f] + * Multiple bound states are not supported. + * Components without bound state (i.e., non-binding components) are supported. + * + * See @cite LangmuirHD1916. + * @tparam ParamHandler_t Type that can add support for external function dependence + */ +template +class LangmuirHybridDesorbBindingBase : public ParamHandlerBindingModelBase +{ +public: + + LangmuirHybridDesorbBindingBase() { } + virtual ~LangmuirHybridDesorbBindingBase() CADET_NOEXCEPT { } + + static const char* identifier() { return ParamHandler_t::identifier(); } + + virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; } + + virtual void timeDerivativeQuasiStationaryFluxes(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yCp, double const* y, double* dResDt, LinearBufferAllocator workSpace) const + { + if (!this->hasQuasiStationaryReactions()) + return; + + if (!ParamHandler_t::dependsOnTime()) + return; + + // Update external function and compute time derivative of parameters + typename ParamHandler_t::ParamsHandle p; + typename ParamHandler_t::ParamsHandle dpDt; + std::tie(p, dpDt) = _paramHandler.updateTimeDerivative(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + // Protein fluxes: -k_{a,i} * c_{p,i} * q_{max,i} * (1 - \sum_j q_j / q_{max,j}) + k_{d,i} * q_i + double qSum = 1.0; + double qSumT = 0.0; + unsigned int bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + const double summand = y[bndIdx] / static_cast(p->qMax[i]); + qSum -= summand; + qSumT += summand / static_cast(p->qMax[i]) * static_cast(dpDt->qMax[i]); + + // Next bound component + ++bndIdx; + } + + bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + // Residual + dResDt[bndIdx] = static_cast(dpDt->kD[i]) * y[bndIdx] + - yCp[i] * (static_cast(dpDt->kA[i]) * static_cast(p->qMax[i]) * qSum + + static_cast(p->kA[i]) * static_cast(dpDt->qMax[i]) * qSum + + static_cast(p->kA[i]) * static_cast(p->qMax[i]) * qSumT); + + // Next bound component + ++bndIdx; + } + } + + CADET_BINDINGMODELBASE_BOILERPLATE + +protected: + using ParamHandlerBindingModelBase::_paramHandler; + using ParamHandlerBindingModelBase::_reactionQuasistationarity; + using ParamHandlerBindingModelBase::_nComp; + using ParamHandlerBindingModelBase::_nBoundStates; + + template + int fluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* y, + CpStateType const* yCp, ResidualType* res, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + // Protein fluxes: -k_{a,i} * c_{p,i} * q_{max,i} * (1 - \sum_j q_j / q_{max,j}) + k_{d,i} * q_i + ResidualType qSum = 1.0; + unsigned int bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + qSum -= y[bndIdx] / static_cast(p->qMax[i]); + + // Next bound component + ++bndIdx; + } + + bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + // Residual + // ANN + const double p1 = static_cast(p->p1[i]); + const double p2 = static_cast(p->p2[i]); + const double p3 = static_cast(p->p3[i]); + const double p4 = static_cast(p->p4[i]); + const double p5 = static_cast(p->p5[i]); + const double p6 = static_cast(p->p6[i]); + const double p7 = static_cast(p->p7[i]); + const double p8 = static_cast(p->p8[i]); + const double p9 = static_cast(p->p9[i]); + const double p10 = static_cast(p->p10[i]); + const double p11 = static_cast(p->p11[i]); + const double P_norm = static_cast(p->P_norm[i]); + const double q_norm = static_cast(p->q_norm[i]); + const double S_norm = static_cast(p->S_norm[i]); + + // const double CF1 = exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)); + // const double CF2 = exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8)); + + // double ANN = (p11 + (p9 / (1 + CF1)) + (p10 / (1 + CF2))); + + res[bndIdx] = static_cast(p->kD[i]) * y[bndIdx] - static_cast(p->kA[i]) * yCp[i] * static_cast(p->qMax[i]) * qSum; + + // Next bound component + ++bndIdx; + } + + return 0; + } + + template + void jacobianImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double const* yCp, int offsetCp, RowIterator jac, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + // Protein fluxes: -k_{a,i} * c_{p,i} * q_{max,i} * (1 - \sum_j q_j / q_{max,j}) + k_{d,i} * q_i + double qSum = 1.0; + int bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + qSum -= y[bndIdx] / static_cast(p->qMax[i]); + + // Next bound component + ++bndIdx; + } + + bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + // Langmuir Params + const double ka = static_cast(p->kA[i]); + const double kd = static_cast(p->kD[i]); + // ANN Params + const double p1 = static_cast(p->p1[i]); + const double p2 = static_cast(p->p2[i]); + const double p3 = static_cast(p->p3[i]); + const double p4 = static_cast(p->p4[i]); + const double p5 = static_cast(p->p5[i]); + const double p6 = static_cast(p->p6[i]); + const double p7 = static_cast(p->p7[i]); + const double p8 = static_cast(p->p8[i]); + const double p9 = static_cast(p->p9[i]); + const double p10 = static_cast(p->p10[i]); + const double p11 = static_cast(p->p11[i]); + const double P_norm = static_cast(p->P_norm[i]); + const double q_norm = static_cast(p->q_norm[i]); + const double S_norm = static_cast(p->S_norm[i]); + // CF = Common Factor + const double CF1 = exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)); + const double CF2 = exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8)); + const double ANN = (p11 + (p9 / (1 + CF1)) + (p10 / (1 + CF2))); + const double ANN_CF = (ANN / abs(ANN)); + + + // dres_i / dc_{p,i} + const double dc_pi_ANN = (((p10 * p2 * CF2) / (P_norm * pow(CF2 + 1, 2))) + ((p1 * p9 * CF1) / (P_norm * pow(CF1 + 1, 2)))) * ANN_CF; + jac[i - bndIdx - offsetCp] = kd * y[bndIdx] * dc_pi_ANN - ka * static_cast(p->qMax[i]) * qSum; + //jac[i - bndIdx - offsetCp] = ka * static_cast(p->qMax[i]) * qSum; + // Getting to c_{p,i}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0} and a +i to c_{p,i}. + // This means jac[i - bndIdx - offsetCp] corresponds to c_{p,i}. + + + // dres_i / dc_{p,0} + const double dc_p0_ANN = (((p10 * p6 * CF2) / (S_norm * pow(CF2 + 1, 2))) + ((p5 * p9 * CF1) / (S_norm * pow(CF1 + 1, 2)))) * ANN_CF; + jac[-bndIdx - offsetCp] = kd * y[bndIdx] * dc_p0_ANN; + // Getting to c_{p,0}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0}. + // This means jac[bndIdx - offsetCp] corresponds to c_{p,0}. + + + + + // Fill dres_i / dq_j + int bndIdx2 = 0; + for (int j = 0; j < _nComp; ++j) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[j] == 0) + continue; + + // dres_i / dq_j + jac[bndIdx2 - bndIdx] = ka * yCp[i] * static_cast(p->qMax[i]) / static_cast(p->qMax[j]); + // Getting to q_j: -bndIdx takes us to q_0, another +bndIdx2 to q_j. This means jac[bndIdx2 - bndIdx] corresponds to q_j. + + ++bndIdx2; + } + + // Add to dres_i / dq_i + const double dc_qi_ANN = (((p10 * p4 * CF2) / (q_norm * pow(CF2 + 1, 2))) + ((p3 * p9 * CF1) / (q_norm * pow(CF1 + 1, 2)))) * ANN_CF; + jac[0] += kd*abs(ANN) + kd* y[bndIdx] * dc_qi_ANN; + // jac[0] += kd; + // Advance to next flux and Jacobian row + ++bndIdx; + ++jac; + } + } +}; + +typedef LangmuirHybridDesorbBindingBase LangmuirHybridDesorbBinding; +typedef LangmuirHybridDesorbBindingBase ExternalLangmuirHybridDesorbBinding; + +namespace binding +{ + void registerLangmuirHybridDesorbModel(std::unordered_map>& bindings) + { + bindings[LangmuirHybridDesorbBinding::identifier()] = []() { return new LangmuirHybridDesorbBinding(); }; + bindings[ExternalLangmuirHybridDesorbBinding::identifier()] = []() { return new ExternalLangmuirHybridDesorbBinding(); }; + } +} // namespace binding + +} // namespace model + +} // namespace cadet diff --git a/src/libcadet/model/binding/TestBinding.cpp b/src/libcadet/model/binding/TestBinding.cpp index 3e5bea7fc..19bca56ab 100644 --- a/src/libcadet/model/binding/TestBinding.cpp +++ b/src/libcadet/model/binding/TestBinding.cpp @@ -29,13 +29,27 @@ /* { - "name": TestParamHandler", + "name": "TestParamHandler", "externalName": "ExtTestParamHandler", "parameters": [ { "type": "ScalarComponentDependentParameter", "varName": "kA", "confName": "MCL_KA"}, { "type": "ScalarComponentDependentParameter", "varName": "kD", "confName": "MCL_KD"}, - { "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "MCL_QMAX"} + { "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "MCL_QMAX"}, + { "type": "ScalarComponentDependentParameter", "varName": "p1", "confName": "ANN_P1"}, + { "type": "ScalarComponentDependentParameter", "varName": "p2", "confName": "ANN_P2"}, + { "type": "ScalarComponentDependentParameter", "varName": "p3", "confName": "ANN_P3"}, + { "type": "ScalarComponentDependentParameter", "varName": "p4", "confName": "ANN_P4"}, + { "type": "ScalarComponentDependentParameter", "varName": "p5", "confName": "ANN_P5"}, + { "type": "ScalarComponentDependentParameter", "varName": "p6", "confName": "ANN_P6"}, + { "type": "ScalarComponentDependentParameter", "varName": "p7", "confName": "ANN_P7"}, + { "type": "ScalarComponentDependentParameter", "varName": "p8", "confName": "ANN_P8"}, + { "type": "ScalarComponentDependentParameter", "varName": "p9", "confName": "ANN_P9"}, + { "type": "ScalarComponentDependentParameter", "varName": "p10", "confName": "ANN_P10"}, + { "type": "ScalarComponentDependentParameter", "varName": "p11", "confName": "ANN_P11"}, + { "type": "ScalarComponentDependentParameter", "varName": "P_norm", "confName": "ANN_P_NORM"}, + { "type": "ScalarComponentDependentParameter", "varName": "q_norm", "confName": "ANN_Q_NORM"}, + { "type": "ScalarComponentDependentParameter", "varName": "S_norm", "confName": "ANN_S_NORM"} ] } */ @@ -45,6 +59,20 @@ kA = Adsorption rate kD = Desorption rate qMax = Capacity + p1 = ANN param 1 + p2 = ANN param 2 + p3 = ANN param 3 + p4 = ANN param 4 + p5 = ANN param 5 + p6 = ANN param 6 + p7 = ANN param 7 + p8 = ANN param 8 + p9 = ANN param 9 + p10 = ANN param 10 + p11 = ANN param 11 + P_norm = Normative Vaule for Protien pore conc + q_norm = Normative Vaule for Protien adsorbed conc + S_norm = Normative Vaule for Salt pore conc */ namespace cadet @@ -57,7 +85,10 @@ namespace cadet inline bool TestParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) { - if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) + if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _p1.size()) || (_kA.size() != _p2.size()) + || (_kA.size() != _p3.size()) || (_kA.size() != _p4.size()) || (_kA.size() != _p5.size()) || (_kA.size() != _p6.size()) + || (_kA.size() != _p7.size()) || (_kA.size() != _p8.size()) || (_kA.size() != _p9.size()) || (_kA.size() != _p10.size()) + || (_kA.size() != _p11.size()) || (_kA.size() != _P_norm.size()) || (_kA.size() != _q_norm.size()) || (_kA.size() != _S_norm.size()) || (_kA.size() < nComp)) throw InvalidParameterException("MCL_KA, MCL_KD, and MCL_QMAX have to have the same size"); return true; @@ -67,7 +98,10 @@ namespace cadet inline bool ExtTestParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) { - if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) + if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _p1.size()) || (_kA.size() != _p2.size()) + || (_kA.size() != _p3.size()) || (_kA.size() != _p4.size()) || (_kA.size() != _p5.size()) || (_kA.size() != _p6.size()) + || (_kA.size() != _p7.size()) || (_kA.size() != _p8.size()) || (_kA.size() != _p9.size()) || (_kA.size() != _p10.size()) + || (_kA.size() != _p11.size()) || (_kA.size() != _P_norm.size()) || (_kA.size() != _q_norm.size()) || (_kA.size() != _S_norm.size()) || (_kA.size() < nComp)) throw InvalidParameterException("EXT_MCL_KA, EXT_MCL_KD, and EXT_MCL_QMAX have to have the same size"); return true; @@ -136,7 +170,28 @@ namespace cadet continue; // Residual - res[bndIdx] = static_cast(p->kD[i]) * y[bndIdx] - static_cast(p->kA[i]) * yCp[i] * static_cast(p->qMax[i]) * qSum; + const double p1 = static_cast(p->p1[i]); + const double p2 = static_cast(p->p2[i]); + const double p3 = static_cast(p->p3[i]); + const double p4 = static_cast(p->p4[i]); + const double p5 = static_cast(p->p5[i]); + const double p6 = static_cast(p->p6[i]); + const double p7 = static_cast(p->p7[i]); + const double p8 = static_cast(p->p8[i]); + const double p9 = static_cast(p->p9[i]); + const double p10 = static_cast(p->p10[i]); + const double p11 = static_cast(p->p11[i]); + const double P_norm = static_cast(p->P_norm[i]); + const double q_norm = static_cast(p->q_norm[i]); + const double S_norm = static_cast(p->S_norm[i]); + + // const double CF1 = exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)); + // const double CF2 = exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8)); + + // double ANN = (p11 + (p9 / (1 + CF1)) + (p10 / (1 + CF2))); + // ANN = (p11 + (p9 / (1 + exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)))) + (p10 / (1 + exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8))))); + + res[bndIdx] = static_cast(p->kD[i]) * y[bndIdx] * abs((p11 + (p9 / (1 + exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)))) + (p10 / (1 + exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8)))))) - static_cast(p->kA[i]) * yCp[i] * static_cast(p->qMax[i]) * qSum; // Next bound component ++bndIdx; @@ -175,12 +230,42 @@ namespace cadet const double ka = static_cast(p->kA[i]); const double kd = static_cast(p->kD[i]); + // ANN Params + const double p1 = static_cast(p->p1[i]); + const double p2 = static_cast(p->p2[i]); + const double p3 = static_cast(p->p3[i]); + const double p4 = static_cast(p->p4[i]); + const double p5 = static_cast(p->p5[i]); + const double p6 = static_cast(p->p6[i]); + const double p7 = static_cast(p->p7[i]); + const double p8 = static_cast(p->p8[i]); + const double p9 = static_cast(p->p9[i]); + const double p10 = static_cast(p->p10[i]); + const double p11 = static_cast(p->p11[i]); + const double P_norm = static_cast(p->P_norm[i]); + const double q_norm = static_cast(p->q_norm[i]); + const double S_norm = static_cast(p->S_norm[i]); + // CF = Common Factor + // const double CF1 = exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)); + // const double CF2 = exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8)); + // const double ANN = (p11 + (p9 / (1 + CF1)) + (p10 / (1 + CF2))); + // const double ANN_CF = (ANN / abs(ANN)); // dres_i / dc_{p,i} - jac[i - bndIdx - offsetCp] = -static_cast(p->kA[i]) * static_cast(p->qMax[i]) * qSum; + // const double dc_pi_ANN = (((p10 * p2 * CF2) / (P_norm * pow(CF2 + 1, 2))) + ((p1 * p9 * CF1) / (P_norm * pow(CF1 + 1, 2)))) * ANN_CF; + + jac[i - bndIdx - offsetCp] = kd * y[bndIdx] * (((p10 * p2 * exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8))) / (P_norm * pow(exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8)) + 1, 2))) + ((p1 * p9 * exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7))) / (P_norm * pow(exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)) + 1, 2)))) * ((p11 + (p9 / (1 + exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)))) + (p10 / (1 + exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8))))) / abs((p11 + (p9 / (1 + exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)))) + (p10 / (1 + exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8))))))) - ka * static_cast(p->qMax[i]) * qSum; // Getting to c_{p,i}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0} and a +i to c_{p,i}. // This means jac[i - bndIdx - offsetCp] corresponds to c_{p,i}. + // dres_i / dc_{p,0} + // const double dc_p0_ANN = (((p10 * p6 * CF2) / (S_norm * pow(CF2 + 1, 2))) + ((p5 * p9 * CF1) / (S_norm * pow(CF1 + 1, 2)))) * ANN_CF; + + jac[-bndIdx - offsetCp] = kd * y[bndIdx] * (((p10 * p6 * exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8))) / (S_norm * pow(exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8)) + 1, 2))) + ((p5 * p9 * exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7))) / (S_norm * pow(exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)) + 1, 2)))) * ((p11 + (p9 / (1 + exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)))) + (p10 / (1 + exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8))))) / abs((p11 + (p9 / (1 + exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)))) + (p10 / (1 + exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8))))))); + // Getting to c_{p,0}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0}. + // This means jac[bndIdx - offsetCp] corresponds to c_{p,0}. + + // Fill dres_i / dq_j int bndIdx2 = 0; for (int j = 0; j < _nComp; ++j) @@ -197,7 +282,9 @@ namespace cadet } // Add to dres_i / dq_i - jac[0] += kd; + // const double dc_qi_ANN = (((p10 * p4 * CF2) / (q_norm * pow(CF2 + 1, 2))) + ((p3 * p9 * CF1) / (q_norm * pow(CF1 + 1, 2)))) * ANN_CF; + + jac[0] += kd * abs((p11 + (p9 / (1 + exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)))) + (p10 / (1 + exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8)))))) + kd * y[bndIdx] * (((p10 * p4 * exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8))) / (q_norm * pow(exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8)) + 1, 2))) + ((p3 * p9 * exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7))) / (q_norm * pow(exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)) + 1, 2)))) * ((p11 + (p9 / (1 + exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)))) + (p10 / (1 + exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8))))) / abs((p11 + (p9 / (1 + exp(-((p1 * yCp[i] / P_norm) + (p3 * y[bndIdx] / q_norm) + (p5 * yCp[0] / S_norm) + p7)))) + (p10 / (1 + exp(-((p2 * yCp[i] / P_norm) + (p4 * y[bndIdx] / q_norm) + (p6 * yCp[0] / S_norm) + p8))))))); // Advance to next flux and Jacobian row ++bndIdx; diff --git a/src/libcadet/model/binding/TestStericMassActionBinding.cpp b/src/libcadet/model/binding/TestStericMassActionBinding.cpp new file mode 100644 index 000000000..8e9ee0f39 --- /dev/null +++ b/src/libcadet/model/binding/TestStericMassActionBinding.cpp @@ -0,0 +1,334 @@ +// ============================================================================= +// CADET +// +// Copyright © 2008-2024: The CADET Authors +// Please see the AUTHORS and CONTRIBUTORS file. +// +// All rights reserved. This program and the accompanying materials +// are made available under the terms of the GNU Public License v3.0 (or, at +// your option, any later version) which accompanies this distribution, and +// is available at http://www.gnu.org/licenses/gpl.html +// ============================================================================= + +#include "model/binding/BindingModelBase.hpp" +#include "model/ExternalFunctionSupport.hpp" +#include "model/binding/BindingModelMacros.hpp" +#include "model/binding/RefConcentrationSupport.hpp" +#include "model/ModelUtils.hpp" +#include "cadet/Exceptions.hpp" +#include "model/Parameters.hpp" +#include "LocalVector.hpp" +#include "SimulationTypes.hpp" + +#include +#include +#include +#include + +/* +{ + "name": "TestSMAParamHandler", + "externalName": "ExtTestSMAParamHandler", + "parameters": + [ + { "type": "ScalarParameter", "varName": "lambda", "confName": "SMA_LAMBDA"}, + { "type": "ScalarComponentDependentParameter", "varName": "kA", "confName": "SMA_KA"}, + { "type": "ScalarComponentDependentParameter", "varName": "kD", "confName": "SMA_KD"}, + { "type": "ScalarComponentDependentParameter", "varName": "nu", "confName": "SMA_NU"}, + { "type": "ScalarComponentDependentParameter", "varName": "sigma", "confName": "SMA_SIGMA"} + ], + "constantParameters": + [ + { "type": "ReferenceConcentrationParameter", "varName": ["refC0", "refQ"], "objName": "refConcentration", "confPrefix": "SMA_"} + ] +} +*/ + +/* Parameter description + ------------------------ + lambda = Ionic capacity + kA = Adsorption rate + kD = Desorption rate + nu = Characteristic charge + sigma = Steric factor + refC0, refQ = Reference concentrations +*/ + +namespace cadet +{ + +namespace model +{ + +inline const char* TestSMAParamHandler::identifier() CADET_NOEXCEPT { return "TEST_STERIC_MASS_ACTION"; } + +inline bool TestSMAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +{ + if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _sigma.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("SMA_KA, SMA_KD, SMA_NU, and SMA_SIGMA have to have the same size"); + + // Assume monovalent salt ions by default + if (_nu.get()[0] <= 0.0) + _nu.get()[0] = 1.0; + + return true; +} + +inline const char* ExtTestSMAParamHandler::identifier() CADET_NOEXCEPT { return "EXT_TEST_STERIC_MASS_ACTION"; } + +inline bool ExtTestSMAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +{ + if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _sigma.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("EXT_SMA_KA, EXT_SMA_KD, EXT_SMA_NU, and EXT_SMA_SIGMA have to have the same size"); + + // Assume monovalent salt ions by default + if (_nu.base()[0] <= 0.0) + _nu.base()[0] = 1.0; + + return true; +} + + +/** + * @brief Defines the steric mass action binding model + * @details Implements the steric mass action adsorption model: \f[ \begin{align} + * q_0 &= \Lambda - \sum_{j} \nu_j q_j \\ + * \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} \left( \Lambda - \sum_j\left( \nu_j + \sigma_j \right) q_j \right)^{\nu_i} - k_{d,i} q_i c_{p,0}^{\nu_i} + * \end{align} \f] + * Component @c 0 is assumed to be salt. Multiple bound states are not supported. + * Components without bound state (i.e., non-binding components) are supported. + * + * See @cite Brooks1992. + * @tparam ParamHandler_t Type that can add support for external function dependence + */ +template +class TestStericMassActionBindingBase : public ParamHandlerBindingModelBase +{ +public: + + TestStericMassActionBindingBase() { } + virtual ~TestStericMassActionBindingBase() CADET_NOEXCEPT { } + + static const char* identifier() { return ParamHandler_t::identifier(); } + + virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBound, unsigned int const* boundOffset) + { + const bool res = BindingModelBase::configureModelDiscretization(paramProvider, nComp, nBound, boundOffset); + + // Guarantee that salt has exactly one bound state + if (nBound[0] != 1) + throw InvalidParameterException("Steric Mass Action binding model requires exactly one bound state for salt component"); + + // First flux is salt, which is always quasi-stationary + _reactionQuasistationarity[0] = true; + + return res; + } + + virtual bool hasSalt() const CADET_NOEXCEPT { return true; } + virtual bool supportsMultistate() const CADET_NOEXCEPT { return false; } + virtual bool supportsNonBinding() const CADET_NOEXCEPT { return true; } + virtual bool hasQuasiStationaryReactions() const CADET_NOEXCEPT { return true; } + virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; } + + virtual bool preConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + // Compute salt component from given bound states q_j + + // Salt equation: nu_0 * q_0 - Lambda + Sum[nu_j * q_j, j] == 0 + // <=> q_0 == (Lambda - Sum[nu_j * q_j, j]) / nu_0 + y[0] = static_cast(p->lambda); + + unsigned int bndIdx = 1; + for (int j = 1; j < _nComp; ++j) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[j] == 0) + continue; + + y[0] -= static_cast(p->nu[j]) * y[bndIdx]; + + // Next bound component + ++bndIdx; + } + + y[0] /= static_cast(p->nu[0]); + + return true; + } + + virtual void postConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const + { + preConsistentInitialState(t, secIdx, colPos, y, yCp, workSpace); + } + + + CADET_BINDINGMODELBASE_BOILERPLATE + +protected: + using ParamHandlerBindingModelBase::_paramHandler; + using ParamHandlerBindingModelBase::_reactionQuasistationarity; + using ParamHandlerBindingModelBase::_nComp; + using ParamHandlerBindingModelBase::_nBoundStates; + + template + int fluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* y, + CpStateType const* yCp, ResidualType* res, LinearBufferAllocator workSpace) const + { + using CpStateParamType = typename DoubleActivePromoter::type; + using StateParamType = typename DoubleActivePromoter::type; + + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + // Salt flux: nu_0 * q_0 - Lambda + Sum[nu_j * q_j, j] == 0 + // <=> nu_0 * q_0 == Lambda - Sum[nu_j * q_j, j] + // Also compute \bar{q}_0 = nu_0 * q_0 - Sum[sigma_j * q_j, j] + res[0] = static_cast(p->nu[0]) * y[0] - static_cast(p->lambda); + StateParamType q0_bar = static_cast(p->nu[0]) * y[0]; + + unsigned int bndIdx = 1; + for (int j = 1; j < _nComp; ++j) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[j] == 0) + continue; + + res[0] += static_cast(p->nu[j]) * y[bndIdx]; + q0_bar -= static_cast(p->sigma[j]) * y[bndIdx]; + + // Next bound component + ++bndIdx; + } + + const ParamType refC0 = static_cast(p->refC0); + const ParamType refQ = static_cast(p->refQ); + const CpStateParamType yCp0_divRef = yCp[0] / refC0; + const StateParamType q0_bar_divRef = q0_bar / refQ; + + // Protein fluxes: -k_{a,i} * c_{p,i} * \bar{q}_0^{nu_i / nu_0} + k_{d,i} * q_i * c_{p,0}^{nu_i / nu_0} + bndIdx = 1; + for (int i = 1; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + const ParamType nu_over_nu0 = static_cast(p->nu[i]) / static_cast(p->nu[0]); + const CpStateParamType c0_pow_nu = pow(yCp0_divRef, nu_over_nu0); + const StateParamType q0_bar_pow_nu = pow(q0_bar_divRef, nu_over_nu0); + + // Residual + res[bndIdx] = static_cast(p->kD[i]) * y[bndIdx] * c0_pow_nu - static_cast(p->kA[i]) * yCp[i] * q0_bar_pow_nu; + + // Next bound component + ++bndIdx; + } + + return 0; + } + + template + void jacobianImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double const* yCp, int offsetCp, RowIterator jac, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + double q0_bar = static_cast(p->nu[0]) * y[0]; + + // Salt flux: nu_0 * q_0 - Lambda + Sum[nu_j * q_j, j] == 0 + jac[0] = static_cast(p->nu[0]); + int bndIdx = 1; + for (int j = 1; j < _nComp; ++j) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[j] == 0) + continue; + + jac[bndIdx] = static_cast(p->nu[j]); + + // Calculate \bar{q}_0 = nu_0 * q_0 - Sum[sigma_j * q_j, j] + q0_bar -= static_cast(p->sigma[j]) * y[bndIdx]; + + // Next bound component + ++bndIdx; + } + + // Advance to protein fluxes + ++jac; + + const double refC0 = static_cast(p->refC0); + const double refQ = static_cast(p->refQ); + const double yCp0_divRef = yCp[0] / refC0; + const double q0_bar_divRef = q0_bar / refQ; + + // Protein fluxes: -k_{a,i} * c_{p,i} * \bar{q}_0^{nu_i} + k_{d,i} * q_i * c_{p,0}^{nu_i} + // We have already computed \bar{q}_0 in the loop above + bndIdx = 1; + for (int i = 1; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + // Getting to c_{p,0}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0}. This means jac[-bndIdx - offsetCp] corresponds to c_{p,0}. + // Getting to c_{p,i}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0} and a +i to c_{p,i}. + // This means jac[i - bndIdx - offsetCp] corresponds to c_{p,i}. + + const double ka = static_cast(p->kA[i]); + const double kd = static_cast(p->kD[i]); + const double nu = static_cast(p->nu[i]) / static_cast(p->nu[0]); + + const double c0_pow_nu = pow(yCp0_divRef, nu); + const double q0_bar_pow_nu = pow(q0_bar_divRef, nu); + const double c0_pow_nu_m1_divRef = pow(yCp0_divRef, nu - 1.0) / refC0; + const double q0_bar_pow_nu_m1_divRef = nu * pow(q0_bar_divRef, nu - 1.0) / refQ; + + // dres_i / dc_{p,0} + jac[-bndIdx - offsetCp] = kd * y[bndIdx] * nu * c0_pow_nu_m1_divRef; + // dres_i / dc_{p,i} + jac[i - bndIdx - offsetCp] = -ka * q0_bar_pow_nu; + // dres_i / dq_0 + jac[-bndIdx] = -ka * yCp[i] * q0_bar_pow_nu_m1_divRef * static_cast(p->nu[0]); + + // Fill dres_i / dq_j + int bndIdx2 = 1; + for (int j = 1; j < _nComp; ++j) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[j] == 0) + continue; + + // dres_i / dq_j + jac[bndIdx2 - bndIdx] = -ka * yCp[i] * q0_bar_pow_nu_m1_divRef * (-static_cast(p->sigma[j])); + // Getting to q_j: -bndIdx takes us to q_0, another +bndIdx2 to q_j. This means jac[bndIdx2 - bndIdx] corresponds to q_j. + + ++bndIdx2; + } + + // Add to dres_i / dq_i + jac[0] += kd * c0_pow_nu; + + // Advance to next flux and Jacobian row + ++bndIdx; + ++jac; + } + } +}; + + +typedef TestStericMassActionBindingBase TestStericMassActionBinding; +typedef TestStericMassActionBindingBase ExternalTestStericMassActionBinding; + +namespace binding +{ + void registerTestStericMassActionModel(std::unordered_map>& bindings) + { + bindings[TestStericMassActionBinding::identifier()] = []() { return new TestStericMassActionBinding(); }; + bindings[ExternalTestStericMassActionBinding::identifier()] = []() { return new ExternalTestStericMassActionBinding(); }; + } +} // namespace binding + +} // namespace model + +} // namespace cadet