diff --git a/doc/README.md b/doc/README.md index a9375c554..a9ee5ff35 100644 --- a/doc/README.md +++ b/doc/README.md @@ -18,5 +18,4 @@ To build the documentation for all releases and the master branch, run: `sphinx-multiversion ./ ./build/`. -This documentation is published under https://cadet.github.io Any changes to the documentation will automatically be pushed to the github-pages repository (https://github.com/cadet/cadet.github.io) using github actions. diff --git a/doc/interface/binding/bi_steric_mass_action.rst b/doc/interface/binding/bi_steric_mass_action.rst index d48890e6c..7202b544f 100644 --- a/doc/interface/binding/bi_steric_mass_action.rst +++ b/doc/interface/binding/bi_steric_mass_action.rst @@ -12,9 +12,9 @@ Bi Steric Mass Action bound states. Otherwise, the adsorption mode is set for each bound state separately. -=================== ========================= ========================================= -**Type:** int **Range:** {0,1} **Length:** 1/NTOTALBND -=================== ========================= ========================================= +=================== ========================= ======================= +**Type:** int **Range:** {0,1} **Length:** 1/NTOTALBND +=================== ========================= ======================= ``BISMA_KA`` Adsorption rate constants in state-major ordering diff --git a/doc/interface/binding/freundlich_ldf.rst b/doc/interface/binding/freundlich_ldf.rst new file mode 100644 index 000000000..65f256855 --- /dev/null +++ b/doc/interface/binding/freundlich_ldf.rst @@ -0,0 +1,52 @@ +.. _freundlich_ldf_config: + +Freundlich LDF +~~~~~~~~~~~~~~~ + +**Group /input/model/unit_XXX/adsorption – ADSORPTION_MODEL = FREUNDLICH_LDF** + + +``IS_KINETIC`` + Selects kinetic or quasi-stationary adsorption mode: 1 = kinetic, 0 = + quasi-stationary. If a single value is given, the mode is set for all + bound states. Otherwise, the adsorption mode is set for each bound + state separately. + +=================== ========================= ================================== +**Type:** int **Range:** {0,1} **Length:** 1/NTOTALBND +=================== ========================= ================================== + +``FLDF_KKIN`` + Driving force coefficient for each component + + +**Unit:** :math:`s^{-1}` + +=================== ========================= ================================== +**Type:** double **Range:** :math:`\ge 0` **Length:** 1/NTOTALBND +=================== ========================= ================================== + + +``FLDF_KF`` + Freundlich coefficient for each component + +**Unit:** :math:`m_{MP}^3~mol^{-1}` + +=================== ========================= ================================== +**Type:** double **Range:** :math:`\ge 0` **Length:** 1/NTOTALBND +=================== ========================= ================================== + +``FLDF_N`` + Freundlich exponent for each component + +**Unit:** :[-] + +=================== ========================= ================================== +**Type:** double **Range:** :math:`> 0` **Length:** 1/NTOTALBND +=================== ========================= ================================== + + + + + + diff --git a/doc/interface/binding/multi_component_bi_langmuir.rst b/doc/interface/binding/multi_component_bi_langmuir.rst index 22c922456..cdeef95aa 100644 --- a/doc/interface/binding/multi_component_bi_langmuir.rst +++ b/doc/interface/binding/multi_component_bi_langmuir.rst @@ -13,11 +13,11 @@ Multi Component Bi-Langmuir state separately. =================== ========================= ========================================= -**Type:** int **Range:** {0,1} **Length:** 1/NTOTALBND +**Type:** int **Range:** {0,1} **Length:** 1/NTOTALBND =================== ========================= ========================================= ``MCBL_KA`` - Adsorption rate constants in state-major ordering + Adsorption rate constants in state-major ordering (see :ref:`ordering_multi_dimensional_data`) **Unit:** :math:`m_{MP}^3~mol^{-1}~s^{-1}` diff --git a/doc/interface/binding/multi_component_bi_langmuir_ldf.rst b/doc/interface/binding/multi_component_bi_langmuir_ldf.rst new file mode 100644 index 000000000..0ca97bb02 --- /dev/null +++ b/doc/interface/binding/multi_component_bi_langmuir_ldf.rst @@ -0,0 +1,44 @@ +.. _multi_component_bi_langmuir_ldf_config: + +Multi Component Bi-Langmuir LDF +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +**Group /input/model/unit_XXX/adsorption – ADSORPTION_MODEL = MULTI_COMPONENT_BILANGMUIR_LDF** + + +``IS_KINETIC`` + Selects kinetic or quasi-stationary adsorption mode: 1 = kinetic, 0 = + quasi-stationary. If a single value is given, the mode is set for all + bound states. Otherwise, the adsorption mode is set for each bound + state separately. + +=================== ========================= ========================================= +**Type:** int **Range:** {0,1} **Length:** 1/NTOTALBND +=================== ========================= ========================================= + +``MCBLLDF_KEQ`` + Equillibrium loading constants in state-major ordering (see :ref:`ordering_multi_dimensional_data`) + +**Unit:** :math:`m_{MP}^3~mol^{-1}` + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP :math:`\cdot` NSTATES +=================== ========================= ========================================= + +``MCBLLDF_KKIN`` + Linear driving force coefficients in state-major ordering + +**Unit:** :math:`s^{-1}` + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP :math:`\cdot` NSTATES +=================== ========================= ========================================= + +``MCBLLDF_QMAX`` + Maximum adsorption capacities in state-major ordering + +**Unit:** :math:`mol~m_{SP}^{-3}` + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\gt 0` **Length:** NCOMP :math:`\cdot` NSTATES +=================== ========================= ========================================= diff --git a/doc/interface/binding/multi_component_langmuir_ldf.rst b/doc/interface/binding/multi_component_langmuir_ldf.rst new file mode 100644 index 000000000..34df73b3f --- /dev/null +++ b/doc/interface/binding/multi_component_langmuir_ldf.rst @@ -0,0 +1,44 @@ +.. _multi_component_langmuir_ldf_config: + +Multi Component Langmuir LDF +============================ + +**Group /input/model/unit_XXX/adsorption – ADSORPTION_MODEL = MULTI_COMPONENT_LANGMUIR_LDF** + + +``IS_KINETIC`` + Selects kinetic or quasi-stationary adsorption mode: 1 = kinetic, 0 = + quasi-stationary. If a single value is given, the mode is set for all + bound states. Otherwise, the adsorption mode is set for each bound + state separately. + +=================== ========================= ========================================= +**Type:** int **Range:** {0,1} **Length:** 1/NTOTALBND +=================== ========================= ========================================= + +``MCLLDF_KEQ`` + Equillibrium loading constants + +**Unit:** :math:`m_{MP}^3~mol^{-1}` + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP +=================== ========================= ========================================= + +``MCLLDF_KKIN`` + Linear driving force coefficients + +**Unit:** :math:`s^{-1}` + +=================== ========================= ================================== +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP +=================== ========================= ================================== + +``MCLLDF_QMAX`` + Maximum adsorption capacities + +**Unit:** :math:`mol~m_{SP}^{-3}` + +=================== ========================= ================================== +**Type:** double **Range:** :math:`\gt 0` **Length:** NCOMP +=================== ========================= ================================== diff --git a/doc/interface/binding/multi_component_langmuir_ldf_liquid_phase.rst b/doc/interface/binding/multi_component_langmuir_ldf_liquid_phase.rst new file mode 100644 index 000000000..69bb8aeb4 --- /dev/null +++ b/doc/interface/binding/multi_component_langmuir_ldf_liquid_phase.rst @@ -0,0 +1,44 @@ +.. _multi_component_langmuir_ldf_liquid_phase_config: + +Multi Component Langmuir LDF Liquid Phase +========================================== + +**Group /input/model/unit_XXX/adsorption – ADSORPTION_MODEL = MULTI_COMPONENT_LANGMUIR_LDF_LIQUID_PHASE** + + +``IS_KINETIC`` + Selects kinetic or quasi-stationary adsorption mode: 1 = kinetic, 0 = + quasi-stationary. If a single value is given, the mode is set for all + bound states. Otherwise, the adsorption mode is set for each bound + state separately. + +=================== ========================= ========================================= +**Type:** int **Range:** {0,1} **Length:** 1/NTOTALBND +=================== ========================= ========================================= + +``MCLLDFC_KEQ`` + Equillibrium loading constants + +**Unit:** :math:`m_{MP}^3~mol^{-1}` + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP +=================== ========================= ========================================= + +``MCLLDFC_KKIN`` + Linear driving force coefficients + +**Unit:** :math:`s^{-1}` + +=================== ========================= ================================== +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP +=================== ========================= ================================== + +``MCLLDFC_QMAX`` + Maximum adsorption capacities + +**Unit:** :math:`mol~m_{SP}^{-3}` + +=================== ========================= ================================== +**Type:** double **Range:** :math:`\gt 0` **Length:** NCOMP +=================== ========================= ================================== diff --git a/doc/interface/introduction.rst b/doc/interface/introduction.rst index 2f25903e2..b0f7ac9c6 100644 --- a/doc/interface/introduction.rst +++ b/doc/interface/introduction.rst @@ -77,6 +77,7 @@ Common notation and identifiers that are used in the subsequent description are * - PARAM_VALUE - Value of a generic unspecified parameter +.. _ordering_multi_dimensional_data: Ordering of multi dimensional data ---------------------------------- @@ -85,14 +86,21 @@ Some model parameters, especially in certain binding models, require multi dimen Since CADET only reads one dimensional arrays, the layout of the data has to be specified (i.e., the way how the data is linearized in memory). The term “*xyz*-major” means that the index corresponding to *xyz* changes the slowest. -For instance, suppose a model with :math:`2` components and :math:`3` bound states has a “component-major” dataset. -Then, the requested matrix is stored in memory such that all bound states are listed for each component (i.e., the component index changes the slowest and the bound state index the fastest): +For instance, suppose a model with :math:`2` components and :math:`3` bound states has a “state-major” dataset. +Then, the requested matrix is stored in memory such that all components are listed for each bound state (i.e., the bound state index changes the slowest and the component index the fastest): :: - comp0bnd0, comp0bnd1, comp0bnd2, comp1bnd0, comp1bnd1, comp1bnd2. + comp0bnd0, comp1bnd0, comp0bnd1, comp1bnd1, comp0bnd2, comp1bnd2 -This linear array can be represented as a :math:`2 \times 3` matrix in “row-major” storage format, or a :math:`3 \times 2` matrix in “column-major” ordering. + +This linear array can also be represented as a :math:`3 \times 2` matrix in “row-major” storage format: + +:: + + comp0bnd0, comp1bnd0 + comp0bnd1, comp1bnd1 + comp0bnd2, comp1bnd2 .. _section_dependent_parameters: diff --git a/doc/literature.bib b/doc/literature.bib index 7ec8aaca5..95c12f1fc 100644 --- a/doc/literature.bib +++ b/doc/literature.bib @@ -374,4 +374,43 @@ @article{Huuk2017 url = {http://doi.wiley.com/10.1002/biot.201600336}, volume = {12}, year = {2017} +} +@article{Alpay1992, +title = {The linear driving force model for fast-cycle adsorption and desorption in a spherical particle}, +journal = {Chemical Engineering Science}, +volume = {47}, +number = {2}, +pages = {499-502}, +year = {1992}, +issn = {0009-2509}, +doi = {https://doi.org/10.1016/0009-2509(92)80041-A}, +url = {https://www.sciencedirect.com/science/article/pii/000925099280041A}, +author = {E. Alpay and D.M. Scott} +} +@incollection{Singh2016, +title = {Chapter 8 - Nanoparticle Ecotoxicology}, +editor = {Ashok K. Singh}, +booktitle = {Engineered Nanoparticles}, +publisher = {Academic Press}, +address = {Boston}, +pages = {343-450}, +year = {2016}, +isbn = {978-0-12-801406-6}, +doi = {https://doi.org/10.1016/B978-0-12-801406-6.00008-X}, +url = {https://www.sciencedirect.com/science/article/pii/B978012801406600008X}, +author = {Ashok K. Singh}, +keywords = {Anthropogenic sphere, Atmosphere, Biosphere, Climate, Environment, Fate of nanoparticles, Human activity, Hydrosphere, Lithosphere, Man made, Nanoparticles, Nanotoxicology, Natural}, +abstract = {Although the physicochemical and beneficial properties of nanoparticles have drawn widespread attention in recent years, their potential for causing environmental damage and ecological toxicity (ecotoxicity) is not fully understood. This chapter provides contemporary information regarding characteristics of the environment (atmosphere, hydrosphere, lithosphere, and biosphere), distribution of bulk and nanoparticles in different spheres, fate of nanoparticles (in air, water, and soil), and adverse effects of different groups of nanoparticles on various types of ecosystems. Although many nanoparticles have been shown to cause potent ecological damage, the results from studies often differ, possibly due to heterogeneity in the physical properties of nanoparticles and the environment, the species and habitat differences, nanoparticles' unique electronic and chemical properties, differences in biomarkers used, and presence of toxins and other contaminants associated with the nanoparticle that may modulate its toxicity. Elucidating these issues may represent a critical future research direction to establish the safe use of nanoparticles.} +} +@article{Benedikt2019, +author = {Aumeier, Benedikt M. and Dang, Anh H. Q. and Ohs, Burkhard and Yüce, Süleyman and Wessling, Matthias}, +title = {Aqueous-Phase Temperature Swing Adsorption for Pesticide Removal}, +journal = {Environmental Science \& Technology}, +volume = {53}, +number = {2}, +pages = {919-927}, +year = {2019}, +doi = {10.1021/acs.est.8b05873}, +URL = { https://doi.org/10.1021/acs.est.8b05873}, +eprint = { https://doi.org/10.1021/acs.est.8b05873} } \ No newline at end of file diff --git a/doc/modelling/binding/freundlich_ldf.rst b/doc/modelling/binding/freundlich_ldf.rst new file mode 100644 index 000000000..4351ce196 --- /dev/null +++ b/doc/modelling/binding/freundlich_ldf.rst @@ -0,0 +1,44 @@ +.. _freundlich_ldf_model: + +Freundlich LDF +~~~~~~~~~~~~~~~ + +The Freundlich isotherm model is an empirical model that considers changes in the equilibrium constant of the binding process due to the heterogeneity of the surface and the variation of the interaction strength :cite:`Benedikt2019,Singh2016`. +This variant of the model is based on the linear driving force approximation (see section :ref:`ldf_model`) and is given as + +.. math:: + \begin{aligned} + q^*_i= k_{F,i}c_{p,i}^{\frac{1}{n_{i}}} && i = 0, \dots, N_{\text{comp}} - 1. + \end{aligned} + +No interaction between the components is considered when the model has multiple components. +One of the limitation of this isotherm is the first order Jacobian :math:`\left(\frac{dq^*}{dc_p}\right)` tends to infinity as :math:`c_{p} \rightarrow 0` for :math:`n>1`. +To address this issue an approximation of isotherm is considered near the origin. +This approximation matches the isotherm in such a way that :math:`q=0` at :math:`c_p=0` and also matches the first derivative of the istotherm at :math:`c_p = \epsilon`, where :math:`\epsilon` is a very small number, for example :math:`1e-14`. +The form of approximation and its derivative is given below for :math:`c_p < \varepsilon` and :math:`n>1`: + +.. math:: + + \begin{aligned} + q^* = \alpha_0+\alpha_1 c+\alpha_2 c_p^2 + \end{aligned} + + \begin{aligned} + \frac{dq^*}{dc_p} = \alpha_1+ 2 \alpha_2 c_p + \end{aligned} + +where :math:`\alpha_0=0` and :math:`\alpha_1` and :math:`\alpha_2` are determined from :math:`\alpha_1 \varepsilon+\alpha_2 \varepsilon^2 = k_f \varepsilon^{1/n}` and :math:`\alpha_1 + 2 \alpha_2 \varepsilon = \frac{1}{n}k_f \varepsilon^{\frac{1-n}{n}}`. + +.. math:: + \begin{aligned} + \alpha_1 = \frac{2 n-1}{n}k_f \varepsilon^{\frac{1-n}{n}} + \end{aligned} +.. math:: + \begin{aligned} + \alpha_2 = \frac{1-n}{n}k_f \varepsilon^{\frac{1-2 n}{n}} + \end{aligned} + +This approximation can be used for any pore phase cocentration :math:`c_p < \epsilon` given :math:`n>1`. +For the case, when :math:`n \le 1` no special treament near the origin is required. +For more information on model parameters required to define in CADET file format, see :ref:`freundlich_ldf_config`. + diff --git a/doc/modelling/binding/index.rst b/doc/modelling/binding/index.rst index 65b444cf1..7a61b4347 100644 --- a/doc/modelling/binding/index.rst +++ b/doc/modelling/binding/index.rst @@ -46,6 +46,34 @@ This can be achieved by a (nonlinear) parameter transform F\left( k_{\text{eq},i}, k_{d,i} \right) = \begin{pmatrix} k_{\text{eq},i} k_{d,i} \\ k_{d,i} \end{pmatrix} \text{ with Jacobian } J_F\left( k_{\text{eq},i}, k_{d,i} \right) = \begin{pmatrix} k_{d,i} & k_{\text{eq},i} \\ 0 & 1 \end{pmatrix}. \end{aligned} + +.. _ldf_model: + +Linear Driving Force (LDF) +--------------------------- +Some authors use the linear driving force (LDF) approximation instead of the native kinetic form of an isotherm. +All three approaches are equivalent in rapid equilibrium (``IS_KINETIC = 0``) but not equivalent when binding kinetics are considered (``IS_KINETIC = 1``). + +1. In the native approach, :math:`\frac{dq}{dt}` is an explicit funtion of :math:`c` and :math:`q`. For example :math:`\frac{dq}{dt}=k_a c (q_m - q)-k_d q` in the Langmuir model. + +2. A linear driving force approximation is based on the equilibrium concentration :math:`q^*` for given :math:`c`. +For example :math:`q^*= \frac{q_m k_{eq} c}{1 + k_{eq} c}` in the Langmuir model. +Here, :math:`\frac{dq}{dt}` is proportional to the actual difference from equilibrium, i.e. :math:`\frac{dq}{dt} = k_{kin}(q^*-q)`. +Note that, the sign of :math:`\frac{dq}{dt}` causes the resulting flux to act towards the equilibrium. +In CADET, this approach is denoted by ``LDF``, for example in ``MULTI_COMPONENT_LANGMUIR_LDF``. + +3. An alterniative linear driving force approximation is based on the equilibrium concentration :math:`c^*` for given :math:`q`. +For example :math:`c^*=\frac{q}{k_{eq} \left(q_{m}-q\right)}` in the Langmuir model. +Here, :math:`\frac{dq}{dt}` is proportional to the actual difference from equilibrium concentrations, i.e. :math:`\frac{dq}{dt} = k_{kin}(c-c^*)`. +Note that, the sign of :math:`\frac{dq}{dt}` causes the resulting flux to act towards the equilibrium. +In CADET, this approach is denoted by ``LDF_LIQUID_PHASE``, for example in ``MULTI_COMPONENT_LANGMUIR_LDF_LIQUID_PHASE``. + +In both LDF examples, the original rate constants :math:`k_a` and :math:`k_d` are replaced by the equilibrium contant :math:`k_{eq}=\frac{k_a}{k_d}`. +The linear driving force approximations are based on a new kinetic constant, :math:`k_{kin}`. + +Note that some isotherms, such as the Freundlich model, don't have a native representation in the above sense. +LDF versions are availabe for some but not all binding models implemented in CADET. + .. _reference_concentrations: Reference concentrations @@ -144,11 +172,26 @@ The models also differ in whether a mobile phase modifier (e.g., salt) is suppor - × - ✓ - × + * - :ref:`freundlich_ldf_model` + - × + - × + - ✓ + - × * - :ref:`multi_component_langmuir_model` - ✓ - × - ✓ - × + * - :ref:`multi_component_langmuir_model_ldf` + - ✓ + - × + - ✓ + - × + * - :ref:`multi_component_langmuir_model_ldf_liquid_phase` + - ✓ + - × + - ✓ + - × * - :ref:`multi_component_anti_langmuir_model` - ✓ - × @@ -189,6 +232,11 @@ The models also differ in whether a mobile phase modifier (e.g., salt) is suppor - × - ✓ - ✓ + * - :ref:`multi_component_bi_langmuir_model_ldf` + - ✓ + - × + - ✓ + - ✓ * - :ref:`multi_component_spreading_model` - ✓ - × diff --git a/doc/modelling/binding/multi_component_bi_langmuir.rst b/doc/modelling/binding/multi_component_bi_langmuir.rst index 540ee7b69..261f25595 100644 --- a/doc/modelling/binding/multi_component_bi_langmuir.rst +++ b/doc/modelling/binding/multi_component_bi_langmuir.rst @@ -3,7 +3,7 @@ Multi Component Bi-Langmuir ~~~~~~~~~~~~~~~~~~~~~~~~~~~ -The multi component Bi-Langmuir model :cite:`Guiochon2006` adds :math:`M - 1` *additional* types of binding sites :math:`q_{i,j}` (:math:`0 \leq j \leq M - 1`) to the Langmuir model (see Section :ref:`multi_component_langmuir_model`) without allowing an exchange between the different sites :math:`q_{i,j}` and :math:`q_{i,k}` (:math:`k \neq j`). +The multi component Bi-Langmuir model :cite:`Guiochon2006` adds :math:`M - 1` additional types of binding sites :math:`q_{i,j}` (:math:`0 \leq j \leq M - 1`) to the Langmuir model (see Section :ref:`multi_component_langmuir_model`) without allowing an exchange between the different sites :math:`q_{i,j}` and :math:`q_{i,k}` (:math:`k \neq j`). Therefore, there are no competitivity effects between the different types of binding sites and they have independent capacities. .. math:: @@ -18,5 +18,4 @@ See the Section :ref:`multi_component_langmuir_model`. Originally, the Bi-Langmuir model is limited to two different binding site types. Here, the model has been extended to arbitrary many binding site types. - For more information on model parameters required to define in CADET file format, see :ref:`multi_component_bi_langmuir_config`. diff --git a/doc/modelling/binding/multi_component_bi_langmuir_ldf.rst b/doc/modelling/binding/multi_component_bi_langmuir_ldf.rst new file mode 100644 index 000000000..fe871973c --- /dev/null +++ b/doc/modelling/binding/multi_component_bi_langmuir_ldf.rst @@ -0,0 +1,15 @@ +.. _multi_component_bi_langmuir_model_ldf: + +Multi Component Bi-Langmuir LDF +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This a linear driving force model variant of the :ref:`multi_component_bi_langmuir_model` model. +It is based on the equilibrium concentration :math:`q^*` for a given liquid phase concentration :math:`c` (see also :ref:`ldf_model`). + +.. math:: + \begin{aligned} + q_{i,j}^*=\frac{q_{m,i,j} k_{eq,i,j} c_i}{1 + \sum_{k=1}^{N_{comp}}{k_{eq,k,j} c_k}} && i = 0, \dots, N_{\text{comp}} - 1, \: j = 0, \dots, M - 1.% (0 \leq i \leq N_{\text{comp}} - 1, \: 0 \leq j \leq M - 1). + \end{aligned} + + +For more information on model parameters required to define in CADET file format, see :ref:`multi_component_bi_langmuir_ldf_config`. diff --git a/doc/modelling/binding/multi_component_langmuir_ldf.rst b/doc/modelling/binding/multi_component_langmuir_ldf.rst new file mode 100644 index 000000000..a356e21c4 --- /dev/null +++ b/doc/modelling/binding/multi_component_langmuir_ldf.rst @@ -0,0 +1,15 @@ +.. _multi_component_langmuir_model_ldf: + +Multi Component Langmuir LDF +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This a linear driving force model variant of the :ref:`multi_component_langmuir_model` model. +It is based on the equilibrium concentration :math:`q^*` for a given liquid phase concentration :math:`c` (see also :ref:`ldf_model`). + +.. math:: + + \begin{aligned} + q_i^*=\frac{q_{m,i} k_{eq,i} c_i}{1 + \sum_{j=1}^{n_{comp}}{k_{eq,j} c_j}} && i = 0, \dots, N_{\text{comp}} - 1. + \end{aligned} + +For more information on model parameters required to define in CADET file format, see :ref:`multi_component_langmuir_ldf_config`. diff --git a/doc/modelling/binding/multi_component_langmuir_ldf_liquid_phase.rst b/doc/modelling/binding/multi_component_langmuir_ldf_liquid_phase.rst new file mode 100644 index 000000000..b8d7ace16 --- /dev/null +++ b/doc/modelling/binding/multi_component_langmuir_ldf_liquid_phase.rst @@ -0,0 +1,16 @@ +.. _multi_component_langmuir_model_ldf_liquid_phase: + +Multi Component Langmuir LDF Liquid Phase +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This a linear driving force model variant of the :ref:`multi_component_langmuir_model` model. +It is based on the equilibrium concentration :math:`c^*` for a given solid phase concentration :math:`q` (see also :ref:`ldf_model`). + +.. math:: + + \begin{aligned} + c_i^*=\frac{q_{i}}{k_{eq,i} q_{m,i} \left(1 - \sum_{j=1}^{N_{\text{comp}}} \frac{q_j}{q_{m,j}}\right) } && i = 0, \dots, N_{\text{comp}} - 1. + \end{aligned} + + +For more information on model parameters required to define in CADET file format, see :ref:`multi_component_langmuir_ldf_liquid_phase_config`. diff --git a/src/libcadet/BindingModelFactory.cpp b/src/libcadet/BindingModelFactory.cpp index ecfe15bf4..5c6eacf63 100644 --- a/src/libcadet/BindingModelFactory.cpp +++ b/src/libcadet/BindingModelFactory.cpp @@ -38,6 +38,10 @@ namespace cadet void registerMultiComponentSpreadingModel(std::unordered_map>& bindings); void registerGeneralizedIonExchangeModel(std::unordered_map>& bindings); void registerColloidalModel(std::unordered_map>& bindings); + void registerFreundlichLDFModel(std::unordered_map>& bindings); + void registerLangmuirLDFModel(std::unordered_map>& bindings); + void registerLangmuirLDFCModel(std::unordered_map>& bindings); + void registerBiLangmuirLDFModel(std::unordered_map>& bindings); } } @@ -61,6 +65,10 @@ namespace cadet model::binding::registerMultiComponentSpreadingModel(_bindingModels); model::binding::registerGeneralizedIonExchangeModel(_bindingModels); model::binding::registerColloidalModel(_bindingModels); + model::binding::registerFreundlichLDFModel(_bindingModels); + model::binding::registerLangmuirLDFModel(_bindingModels); + model::binding::registerLangmuirLDFCModel(_bindingModels); + model::binding::registerBiLangmuirLDFModel(_bindingModels); registerModel(); } diff --git a/src/libcadet/CMakeLists.txt b/src/libcadet/CMakeLists.txt index 9615bafaf..74419153f 100644 --- a/src/libcadet/CMakeLists.txt +++ b/src/libcadet/CMakeLists.txt @@ -94,6 +94,10 @@ set(LIBCADET_BINDINGMODEL_SOURCES ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/MultiComponentSpreadingBinding.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/ColloidalBinding.cpp + ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/FreundlichLDFBinding.cpp + ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/LangmuirLDFBinding.cpp + ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/LangmuirLDFCBinding.cpp + ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/BiLangmuirLDFBinding.cpp ) # LIBCADET_REACTIONMODEL_SOURCES holds all source files of reaction models @@ -219,7 +223,7 @@ if (LAPACK_FOUND) set(LIB_LAPACK_DEFINE "CADET_LAPACK_TRAILING_UNDERSCORE") # Add the build target for the static nonlinalg library - add_library(libcadet_nonlinalg_static STATIC ${LIBCADET_NONLINALG_SOURCES} ${LIBCADET_NONLINALG_SPARSE_SOURCES}) + add_library(libcadet_nonlinalg_static STATIC ${LIBCADET_NONLINALG_SOURCES} ${LIBCADET_NONLINALG_SPARSE_SOURCES} "model/binding/LangmuirLDFCBinding.cpp") set_target_properties(libcadet_nonlinalg_static PROPERTIES OUTPUT_NAME cadet_nonlinalg_static) target_compile_definitions(libcadet_nonlinalg_static PRIVATE libcadet_nonlinalg_static_EXPORTS ${LIB_LAPACK_DEFINE}) target_link_libraries(libcadet_nonlinalg_static PUBLIC CADET::CompileOptions PRIVATE SUNDIALS::sundials_idas ${SUNDIALS_NVEC_TARGET} ${LAPACK_LIBRARIES}) @@ -235,7 +239,7 @@ if (LAPACK_FOUND) # Add the build target for CADET object library - add_library(libcadet_object OBJECT ${LIBCADET_SOURCES}) + add_library(libcadet_object OBJECT ${LIBCADET_SOURCES} "model/binding/LangmuirLDFCBinding.cpp") target_compile_definitions(libcadet_object PRIVATE libcadet_EXPORTS ${LIB_LAPACK_DEFINE}) target_link_libraries(libcadet_object PUBLIC CADET::CompileOptions CADET::LibOptions PRIVATE CADET::AD libcadet_nonlinalg_static SUNDIALS::sundials_idas ${SUNDIALS_NVEC_TARGET} ${TBB_TARGET}) @@ -362,4 +366,4 @@ install(DIRECTORY ${CMAKE_SOURCE_DIR}/include/cadet DESTINATION include) unset(LIBCADET_TARGETS) # Info message -message(STATUS "Added LIBCADET module") +message(STATUS "Added LIBCADET module") \ No newline at end of file diff --git a/src/libcadet/model/binding/BiLangmuirBinding.cpp b/src/libcadet/model/binding/BiLangmuirBinding.cpp index 8e949e3f1..61cfe4bdb 100644 --- a/src/libcadet/model/binding/BiLangmuirBinding.cpp +++ b/src/libcadet/model/binding/BiLangmuirBinding.cpp @@ -112,7 +112,7 @@ class BiLangmuirBindingBase : public ParamHandlerBindingModelBase +#include +#include + +/* +{ + "name": "BiLangmuirLDFParamHandler", + "externalName": "ExtBiLangmuirLDFParamHandler", + "parameters": + [ + { "type": "ComponentBoundStateMajorDependentParameter", "varName": "keq", "confName": "MCBLLDF_KEQ"}, + { "type": "ComponentBoundStateMajorDependentParameter", "varName": "kkin", "confName": "MCBLLDF_KKIN"}, + { "type": "ComponentBoundStateMajorDependentParameter", "varName": "qMax", "confName": "MCBLLDF_QMAX"} + ] +} +*/ + +/* Parameter description + ------------------------ + keq = Equillibrium constant in binding site-major ordering + kkkin = Linear driving force coefficient in binding site-major ordering + qMax = Capacity in binding site-major ordering +*/ + +namespace cadet +{ + + namespace model + { + + inline const char* BiLangmuirLDFParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_BILANGMUIR_LDF"; } + + inline bool BiLangmuirLDFParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) + { + if ((_keq.size() != _kkin.size()) || (_keq.size() != _qMax.size()) || (_keq.size() < nComp)) + throw InvalidParameterException("MCBLLDF_KEQ, MCBLLDF_KKIN, and MCBLLDF_QMAX have to have the same size"); + + return true; + } + + inline const char* ExtBiLangmuirLDFParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_BILANGMUIR_LDF"; } + + inline bool ExtBiLangmuirLDFParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) + { + if ((_keq.size() != _kkin.size()) || (_keq.size() != _qMax.size()) || (_keq.size() < nComp)) + throw InvalidParameterException("EXT_MCBLLDF_KEQ, EXT_MCBLLDF_KKIN, and EXT_MCBLLDF_QMAX have to have the same size"); + + return true; + } + + + /** + * @brief Defines the multi component Bi-Langmuir binding model + * @details Implements the Bi-Langmuir adsorption model: \f[ \begin{align} + * \frac{\mathrm{d}q_i^A}{\mathrm{d}t} &= k_{a,i}^A c_{p,i} q_{\text{max},i}^A \left( 1 - \sum_j \frac{q_j^A}{q_{\text{max},j}^A} \right) - k_{d,i}^A q_i^A \\ + * \frac{\mathrm{d}q_i^B}{\mathrm{d}t} &= k_{a,i}^B c_{p,i} q_{\text{max},i}^B \left( 1 - \sum_j \frac{q_j^B}{q_{\text{max},j}^B} \right) - k_{d,i}^B q_i^B \\ + * \vodts & \vdots + * \end{align} \f] + * Here, several different types of binding sites @f$ q^A @f$, @f$ q^B @f$, etc. are considered. A molecule can either bind to + * site A or B (or C, etc.). A direct exchange between the different binding sites does not occur. + * While components without bound state (i.e., non-binding components) are supported, all other components must have + * the same number of bound states (i.e., binding sites). + * + * Internal state vector order is component-major. The state vector is composed of all components and within each component + * all bound states are listed. + * + * See @cite Guiochon2006. + * @tparam ParamHandler_t Type that can add support for external function dependence + */ + template + class BiLangmuirLDFBindingBase : public ParamHandlerBindingModelBase + { + public: + + BiLangmuirLDFBindingBase() { } + virtual ~BiLangmuirLDFBindingBase() 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); + + unsigned int numSlices = 0; + for (unsigned int i = 0; i < nComp; ++i) + { + if (nBound[i] == 0) + continue; + + if (numSlices == 0) + numSlices = nBound[i]; + + if (nBound[i] != numSlices) + throw InvalidParameterException("Bi-Langmuir LDF binding model requires all components to have the same number of bound states or zero"); + } + + _numBindingComp = numBindingComponents(_nBoundStates, _nComp); + + // Allocate space for parameters + _paramHandler.reserve(nComp * numSlices, numSlices, nComp, nBound); + + return res; + } + + virtual bool hasSalt() const CADET_NOEXCEPT { return false; } + virtual bool supportsMultistate() const CADET_NOEXCEPT { return true; } + virtual bool supportsNonBinding() const CADET_NOEXCEPT { return true; } + 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 flux: -k_{a,i}^j * c_{p,i} * (1 - \sum q_i^j / q_{max,i}^j) + k_{d,i}^j * q_i^j + + const unsigned int nSites = p->kA.slices(); + + // Ordering of the states is (q_{comp,state}) + // q_{0,0}, q{0,1}, q_{0,2}, q_{1,0}, q_{1,1}, q_{1,2}, ... + // A state corresponds to a type of binding site. It is assumed that all components have either 0 + // or the same number of states. Thus, a component is either non-binding or has nSites bound states. + // + // The same ordering is used for the fluxes. That is, q_{0,0}, q_{1,0} and q_{0,1}, q_{1,1} and ... each + // form one Langmuir binding model system. + + // Loop over all binding site types + for (unsigned int site = 0; site < nSites; ++site, ++y, ++dResDt) + { + // y, and dResDt point to q_{0,site} + + // Get parameter slice for current binding site type + active const* const localKa = p->kA[site]; + // active const* const localKd = p->kD[site]; + active const* const localQmax = p->qMax[site]; + active const* const localKaT = dpDt->kA[site]; + active const* const localKdT = dpDt->kD[site]; + active const* const localQmaxT = dpDt->qMax[site]; + + double qSum = 1.0; + double qSumT = 0.0; + unsigned int bndIdx = 0; + + // bndIdx is used as a counter inside one binding site type + // Getting from one component to another requires a step size of nSites (stride) + + 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 * nSites] / static_cast(localQmax[i]); + qSum -= summand; + qSumT += summand / static_cast(localQmax[i]) * static_cast(localQmaxT[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 * nSites] = static_cast(localKdT[i]) * y[bndIdx * nSites] + - yCp[i] * (static_cast(localKaT[i]) * static_cast(localQmax[i]) * qSum + + static_cast(localKa[i]) * static_cast(localQmaxT[i]) * qSum + + static_cast(localKa[i]) * static_cast(localQmax[i]) * qSumT); + + // Next bound component + ++bndIdx; + } + } + }*/ + + CADET_BINDINGMODELBASE_BOILERPLATE + + protected: + using ParamHandlerBindingModelBase::_paramHandler; + using ParamHandlerBindingModelBase::_reactionQuasistationarity; + using ParamHandlerBindingModelBase::_nComp; + using ParamHandlerBindingModelBase::_nBoundStates; + + unsigned int _numBindingComp; //!< Number of binding components + + 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 flux: \frac{dq_{i,j}}{dt} = k_{kin,i,j}(q_{i,j}^*-q_{i,j}) with + //q_{i,j}^*=\frac{q_{m,i,j} k_{eq,i,j} c_i}{1 + \sum_{k=1}^{n_{comp}}{k_{eq,k,j} c_k}} + + const unsigned int nSites = p->keq.slices(); + + // Ordering of the states is (q_{comp,state}) + // q_{0,0}, q{0,1}, q_{0,2}, q_{1,0}, q_{1,1}, q_{1,2}, ... + // A state corresponds to a type of binding site. It is assumed that all components have either 0 + // or the same number of states. Thus, a component is either non-binding or has nSites bound states. + // + // The same ordering is used for the fluxes. That is, q_{0,0}, q_{1,0} and q_{0,1}, q_{1,1} and ... each + // form one Langmuir binding model system. + + // Loop over all binding site types + for (unsigned int site = 0; site < nSites; ++site, ++y, ++res) + { + // y, and res point to q_{0,site} + + // Get parameter slice for current binding site type + active const* const localKeq = p->keq[site]; + active const* const localKkin = p->kkin[site]; + active const* const localQmax = p->qMax[site]; + + ResidualType cpSum = 1.0; + unsigned int bndIdx = 0; + + // bndIdx is used as a counter inside one binding site type + // Getting from one component to another requires a step size of nSites (stride) + + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + cpSum += yCp[i] * static_cast(localKeq[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 * nSites] = static_cast(localKkin[i]) * (y[bndIdx * nSites] - static_cast(localKeq[i]) * yCp[i] * static_cast(localQmax[i]) / cpSum); + + // 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 flux: \frac{dq_{i,j}}{dt} = k_{kin,i,j}(q_{i,j}^*-q_{i,j}) with + //q_{i,j}^*=\frac{q_{m,i,j} k_{eq,i,j} c_i}{1 + \sum_{k=1}^{n_{comp}}{k_{eq,k,j} c_k}} + + // Ordering of the states is (q_{comp,state}, example uses 2 components, 3 binding sites) + // q_{0,0}, q{0,1}, q_{0,2}, q_{1,0}, q_{1,1}, q_{1,2}, ... + // Ordering of the fluxes is the same, that is, we need nSites steps to jump from one flux of a Langmuir + // binding model system to the next. + + const int nSites = static_cast(p->keq.slices()); + + // Loop over all binding site types + for (int site = 0; site < nSites; ++site, ++y) + { + // Get parameter slice for current binding site type + active const* const localKeq = p->keq[site]; + active const* const localKkin = p->kkin[site]; + active const* const localQmax = p->qMax[site]; + + double cpSum = 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; + + cpSum += yCp[i] * static_cast(localKeq[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 keq = static_cast(localKeq[i]); + const double kkin = static_cast(localKkin[i]); + + // dres_i / dc_{p,i} + jac[i - site - offsetCp - nSites * bndIdx] = -kkin * keq * static_cast(localQmax[i]) / cpSum; + // Getting to c_{p,i}: -nSites * bndIdx takes us to q_{0,site}, another -site to q_{0,0}. From there, we + // take a -offsetCp to reach c_{p,0} and a +i to arrive at c_{p,i}. + // This means jac[i - site - offsetCp - nSites * bndIdx] corresponds to c_{p,i}. + + // Fill dres_i / dc_j + const double commonFactor = kkin * keq * static_cast(localQmax[i]) * yCp[i] / (cpSum * cpSum); + int bndIdx2 = 0; + for (int j = 0; j < _nComp; ++j) + { + // Skip components without bound states (bound state index bndIdx2 is not advanced) + if (_nBoundStates[j] == 0) + continue; + + // dres_i / dc_j + jac[j - site - offsetCp - nSites * bndIdx] += static_cast(localKeq[j]) * commonFactor; + // Getting to c_{p,j}: -nSites * bndIdx takes us to q_{0,site}, another -site to q_{0,0}. From there, we + // take a -offsetCp to reach c_{p,0} and a +j to arrive at c_{p,j}. + // This means jac[j - site - offsetCp - nSites * bndIdx] corresponds to c_{p,j}. + + ++bndIdx2; + } + + // Add to dres_i / dq_{i,site} + jac[0] += kkin; + + // Advance to next flux and Jacobian row + ++bndIdx; + jac += nSites; + // Note that there is a spacing of nSites between the fluxes inside one binding site type + } + // We are at the end of the flux block q_{_numBindingComp,site} and need to jump back to the beginning + // by using -_numBindingComp * nSites steps. From there, we take one step forward to arrive at q_{0,site+1}. + jac -= _numBindingComp * nSites - 1; + } + } + }; + + typedef BiLangmuirLDFBindingBase BiLangmuirLDFBinding; + typedef BiLangmuirLDFBindingBase ExternalBiLangmuirLDFBinding; + + namespace binding + { + void registerBiLangmuirLDFModel(std::unordered_map>& bindings) + { + bindings[BiLangmuirLDFBinding::identifier()] = []() { return new BiLangmuirLDFBinding(); }; + bindings[ExternalBiLangmuirLDFBinding::identifier()] = []() { return new ExternalBiLangmuirLDFBinding(); }; + } + } // namespace binding + + } // namespace model + +} // namespace cadet diff --git a/src/libcadet/model/binding/FreundlichLDFBinding.cpp b/src/libcadet/model/binding/FreundlichLDFBinding.cpp new file mode 100644 index 000000000..d28bbba90 --- /dev/null +++ b/src/libcadet/model/binding/FreundlichLDFBinding.cpp @@ -0,0 +1,198 @@ +// ============================================================================= +// CADET +// +// Copyright © 2008-2022: 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 +#include +/* +{ + "name": "FreundlichLDFParamHandler", + "externalName": "ExtFreundlichLDFParamHandler", + "parameters": + [ + { "type": "ScalarComponentDependentParameter", "varName": "kkin", "confName": "FLDF_KKIN"}, + { "type": "ScalarComponentDependentParameter", "varName": "kF", "confName": "FLDF_KF"}, + { "type": "ScalarComponentDependentParameter", "varName": "n", "confName": "FLDF_N"} + ] +} +*/ + +/* Parameter description + ------------------------ + k_kin : Kinetic coefficient + k_F : Freundlich coefficient + n : Freundlich exponent +*/ + +namespace cadet +{ + + namespace model + { + + inline const char* FreundlichLDFParamHandler::identifier() CADET_NOEXCEPT { return "FREUNDLICH_LDF"; } + + inline bool FreundlichLDFParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) + { + if ((_kkin.size() != _kF.size()) || (_kkin.size() != _n.size()) || (_kkin.size() < nComp)) + throw InvalidParameterException("FLDF_KKIN, FLDF_KF and FLDF_N, have to have the same size"); + + return true; + } + + inline const char* ExtFreundlichLDFParamHandler::identifier() CADET_NOEXCEPT { return "EXT_FREUNDLICH_LDF"; } + + inline bool ExtFreundlichLDFParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) + { + if ((_kkin.size() != _kF.size()) || (_kkin.size() != _n.size()) || (_kkin.size() < nComp)) + throw InvalidParameterException("FLDF_KKIN, FLDF_KF and FLDF_N, have to have the same size"); + + + return true; + } + + + /** + * @brief Defines the multi component Freundlich LDF (Linear Driving Force) isotherm + * @details Implements the Freundlich LDF adsorption model: \f[ \begin{align} + * \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{f,i} c_{p,i}^{\frac{1}{n_{i}}} \end{align} \f] + * Components without bound state (i.e., non-binding components) are supported. + * @tparam ParamHandler_t Type that can add support for external function dependence + */ + template + class FreundlichLDFBindingBase : public ParamHandlerBindingModelBase + { + public: + + FreundlichLDFBindingBase() { } + virtual ~FreundlichLDFBindingBase() 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; + + const double _threshold = 1e-14; + + template + int fluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* y, + CpStateType const* yCp, ResidualType* res, LinearBufferAllocator workSpace) const + { + using std::abs; + using std::pow; + + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + // Protein Flux: dq/dt = k_kin * (q* - q) with q* = k_F * c^(1/n) + 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 ParamType n_param = static_cast(p->n[i]); + const ParamType kF = static_cast(p->kF[i]); + const ParamType kkin = static_cast(p->kkin[i]); + + // Residual + if ((n_param > 1) && (abs(yCp[i]) < _threshold)) + { + const ParamType alpha_1 = ((2.0 * n_param - 1.0) / n_param) * kF * pow(_threshold, (1.0 - n_param) / n_param); + const ParamType alpha_2 = ((1.0 - n_param) / n_param) * kF * pow(_threshold, (1.0 - 2.0 * n_param) / n_param); + res[bndIdx] = kkin * (y[bndIdx] - yCp[i] * (alpha_1 + alpha_2 * yCp[i])); + } + else + { + res[bndIdx] = kkin * (y[bndIdx] - kF * pow(abs(yCp[i]), 1.0 / n_param)); + } + + // 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 + { + using std::abs; + using std::pow; + + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + 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 n_param = static_cast(p->n[i]); + const double kF = static_cast(p->kF[i]); + const double kkin = static_cast(p->kkin[i]); + + // dres / dq_i + jac[0] = kkin; + + // dres / dc_{p,i} + if ((n_param > 1) && (abs(yCp[i]) < _threshold)) + { + double const alpha_1 = ((2.0 * n_param - 1.0) / n_param) * kF * pow(_threshold, (1.0 - n_param) / n_param); + double const alpha_2 = ((1.0 - n_param) / n_param) * kF * pow(_threshold, (1.0 - 2.0 * n_param) / n_param); + jac[i - bndIdx - offsetCp] = -kkin * (alpha_1 + 2.0 * alpha_2 * yCp[i]); + } + else + { + jac[i - bndIdx - offsetCp] = -(1.0 / n_param) * kkin * kF * pow(abs(yCp[i]), (1.0 - n_param) / n_param); + } + + ++bndIdx; + ++jac; + } + } + + }; + + typedef FreundlichLDFBindingBase FreundlichLDFBinding; + typedef FreundlichLDFBindingBase ExternalFreundlichLDFBinding; + + namespace binding + { + void registerFreundlichLDFModel(std::unordered_map>& bindings) + { + bindings[FreundlichLDFBinding::identifier()] = []() { return new FreundlichLDFBinding(); }; + bindings[ExternalFreundlichLDFBinding::identifier()] = []() { return new ExternalFreundlichLDFBinding(); }; + } + } // namespace binding + + } // namespace model + +} // namespace cadet diff --git a/src/libcadet/model/binding/LangmuirLDFBinding.cpp b/src/libcadet/model/binding/LangmuirLDFBinding.cpp new file mode 100644 index 000000000..3a48a3f23 --- /dev/null +++ b/src/libcadet/model/binding/LangmuirLDFBinding.cpp @@ -0,0 +1,268 @@ +// ============================================================================= +// CADET +// +// Copyright © 2008-2022: 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": "LangmuirLDFParamHandler", + "externalName": "ExtLangmuirLDFParamHandler", + "parameters": + [ + { "type": "ScalarComponentDependentParameter", "varName": "keq", "confName": "MCLLDF_KEQ"}, + { "type": "ScalarComponentDependentParameter", "varName": "kkin", "confName": "MCLLDF_KKIN"}, + { "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "MCLLDF_QMAX"} + ] +} +*/ + +/* Parameter description + ------------------------ + keq = Equillibrium constant + kkin = Linear driving force + qMax = Capacity +*/ + +namespace cadet +{ + +namespace model +{ + +inline const char* LangmuirLDFParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_LANGMUIR_LDF"; } + +inline bool LangmuirLDFParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +{ + if ((_keq.size() != _kkin.size()) || (_keq.size() != _qMax.size()) || (_keq.size() < nComp)) + throw InvalidParameterException("MCLLDF_KEQ, MCLLDF_KKIN, and MCLLDF_QMAX have to have the same size"); + + return true; +} + +inline const char* ExtLangmuirLDFParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_LANGMUIR_LDF"; } + +inline bool ExtLangmuirLDFParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +{ + if ((_keq.size() != _kkin.size()) || (_keq.size() != _qMax.size()) || (_keq.size() < nComp)) + throw InvalidParameterException("EXT_MCLLDF_KEQ, EXT_MCLLDF_KKIN, and EXT_MCLLDF_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 LangmuirLDFBindingBase : public ParamHandlerBindingModelBase +{ +public: + + LangmuirLDFBindingBase() { } + virtual ~LangmuirLDFBindingBase() 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_{kin,i}q_i - k_{kin,i} \frac{q_{m,i}k_{eq,i}c_i}{1+\sum_{j=1}^{n_{comp}} k_{eq,j}c_j} + double cpSum = 1.0; + double cpSumT = 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]); + cpSum -= summand; + cpSumT += 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]) * cpSum + + static_cast(p->kA[i]) * static_cast(dpDt->qMax[i]) * cpSum + + static_cast(p->kA[i]) * static_cast(p->qMax[i]) * cpSumT); + + // 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_{kin,i}q_i - k_{kin,i} \frac{q_{m,i}k_{eq,i}c_i}{1+\sum_{j=1}^{n_{comp}} k_{eq,j}c_j} + ResidualType cpSum = 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; + + cpSum += yCp[i] * static_cast(p->keq[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->kkin[i]) * (y[bndIdx] - static_cast(p->keq[i]) * yCp[i] * static_cast(p->qMax[i]) / cpSum); + + // 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_{kin,i}q_i - k_{kin,i} \frac{q_{m,i}k_{eq,i}c_i}{1+\sum_{j=1}^{n_{comp}} k_{eq,j}c_j} + double cpSum = 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; + + cpSum += yCp[i] * static_cast(p->keq[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 keq = static_cast(p->keq[i]); + const double kkin = static_cast(p->kkin[i]); + //(keq *keq *static_cast(p->qMax[i]) * yCp[i] / (cpSum*cpSum)) + // dres_i / dc_{p,i} + jac[i - bndIdx - offsetCp] = -kkin * keq * static_cast(p->qMax[i]) / cpSum; + // 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 / dc_j + const double commonFactor = keq * kkin * static_cast(p->qMax[i]) * yCp[i] / (cpSum * cpSum); + 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 / dc_j + jac[j - bndIdx - offsetCp] += static_cast(p->keq[j]) * commonFactor; + // Getting to c_{p,j}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0} and a +j to c_{p,j}. + // This means jac[j - bndIdx - offsetCp] corresponds to c_{p,j}. + + + ++bndIdx2; + } + + // Add to dres_i / dq_i + jac[0] += kkin; + + // Advance to next flux and Jacobian row + ++bndIdx; + ++jac; + } + } +}; + +typedef LangmuirLDFBindingBase LangmuirLDFBinding; +typedef LangmuirLDFBindingBase ExternalLangmuirLDFBinding; + +namespace binding +{ + void registerLangmuirLDFModel(std::unordered_map>& bindings) + { + bindings[LangmuirLDFBinding::identifier()] = []() { return new LangmuirLDFBinding(); }; + bindings[ExternalLangmuirLDFBinding::identifier()] = []() { return new ExternalLangmuirLDFBinding(); }; + } +} // namespace binding + +} // namespace model + +} // namespace cadet diff --git a/src/libcadet/model/binding/LangmuirLDFCBinding.cpp b/src/libcadet/model/binding/LangmuirLDFCBinding.cpp new file mode 100644 index 000000000..2cf296cdf --- /dev/null +++ b/src/libcadet/model/binding/LangmuirLDFCBinding.cpp @@ -0,0 +1,267 @@ +// ============================================================================= +// CADET +// +// Copyright © 2008-2022: 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": "LangmuirLDFCParamHandler", + "externalName": "ExtLangmuirLDFCParamHandler", + "parameters": + [ + { "type": "ScalarComponentDependentParameter", "varName": "keq", "confName": "MCLLDFC_KEQ"}, + { "type": "ScalarComponentDependentParameter", "varName": "kkin", "confName": "MCLLDFC_KKIN"}, + { "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "MCLLDFC_QMAX"} + ] +} +*/ + +/* Parameter description + ------------------------ + keq = Equillibrium constant + kkin = Linear driving force + qMax = Capacity +*/ + +namespace cadet +{ + +namespace model +{ + +inline const char* LangmuirLDFCParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_LANGMUIR_LDF_LIQUID_PHASE"; } + +inline bool LangmuirLDFCParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +{ + if ((_keq.size() != _kkin.size()) || (_keq.size() != _qMax.size()) || (_keq.size() < nComp)) + throw InvalidParameterException("MCLLDFC_KEQ, MCLLDFC_KKIN, and MCLLDFC_QMAX have to have the same size"); + + return true; +} + +inline const char* ExtLangmuirLDFCParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_LANGMUIR_LDF_LIQUID_PHASE"; } + +inline bool ExtLangmuirLDFCParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +{ + if ((_keq.size() != _kkin.size()) || (_keq.size() != _qMax.size()) || (_keq.size() < nComp)) + throw InvalidParameterException("EXT_MCLLDFC_KEQ, EXT_MCLLDFC_KKIN, and EXT_MCLLDFC_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 LangmuirLDFCBindingBase : public ParamHandlerBindingModelBase +{ +public: + + LangmuirLDFCBindingBase() { } + virtual ~LangmuirLDFCBindingBase() 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_{kin,i}q_i - k_{kin,i} \frac{q_{m,i}k_{eq,i}c_i}{1+\sum_{j=1}^{n_{comp}} k_{eq,j}c_j} + 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->kkin[i]) * (yCp[i] - y[bndIdx] / (static_cast(p->keq[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_{kin,i}q_i - k_{kin,i} \frac{q_{m,i}k_{eq,i}c_i}{1+\sum_{j=1}^{n_{comp}} k_{eq,j}c_j} + 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 keq = static_cast(p->keq[i]); + const double kkin = static_cast(p->kkin[i]); + const double qMax = static_cast(p->qMax[i]); + + // dres_i / dc_{p,i} + jac[i - bndIdx - offsetCp] = -kkin; + // 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] += y[bndIdx] * kkin / (keq * qMax * static_cast(p->qMax[j]) * qSum * qSum); + // 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] += kkin / (keq * qMax * qSum); + + // Advance to next flux and Jacobian row + ++bndIdx; + ++jac; + } + } +}; + +typedef LangmuirLDFCBindingBase LangmuirLDFCBinding; +typedef LangmuirLDFCBindingBase ExternalLangmuirLDFCBinding; + +namespace binding +{ + void registerLangmuirLDFCModel(std::unordered_map>& bindings) + { + bindings[LangmuirLDFCBinding::identifier()] = []() { return new LangmuirLDFCBinding(); }; + bindings[ExternalLangmuirLDFCBinding::identifier()] = []() { return new ExternalLangmuirLDFCBinding(); }; + } +} // namespace binding + +} // namespace model + +} // namespace cadet diff --git a/test/BindingModels.cpp b/test/BindingModels.cpp index 0145786bd..cebeae684 100644 --- a/test/BindingModels.cpp +++ b/test/BindingModels.cpp @@ -42,6 +42,42 @@ CADET_BINDINGTEST("LINEAR", "EXT_LINEAR", (1,1), (1,0,1), (1.0, 2.0, 0.0, 0.0), )json", \ 1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_UNUSED, CADET_COMPARE_BINDING_VS_NONBINDING) + CADET_BINDINGTEST("FREUNDLICH_LDF", "EXT_FREUNDLICH_LDF", (1, 1), (1, 0, 1), (1.0, 2.0, 0.0, 0.0), (1.0, 3.0, 2.0, 0.0, 0.0), \ + R"json( "FLDF_KKIN": [1.0, 2.0], + "FLDF_KF": [0.1, 0.2], + "FLDF_N": [0.5, 1.2] + )json", \ + R"json( "FLDF_KKIN": [1.0, 0.5, 2.0], + "FLDF_KF": [0.1, 0.3, 0.2], + "FLDF_N": [0.5, 2.2, 1.2] + )json", \ + R"json( "EXT_FLDF_KKIN": [0.0, 0.0], + "EXT_FLDF_KKIN_T": [1.0, 2.0], + "EXT_FLDF_KKIN_TT": [0.0, 0.0], + "EXT_FLDF_KKIN_TTT": [0.0, 0.0], + "EXT_FLDF_KF": [0.0, 0.0], + "EXT_FLDF_KF_T": [0.1, 0.2], + "EXT_FLDF_KF_TT": [0.0, 0.0], + "EXT_FLDF_KF_TTT": [0.0, 0.0], + "EXT_FLDF_N": [0.0, 0.0], + "EXT_FLDF_N_T": [0.5, 1.2], + "EXT_FLDF_N_TT": [0.0, 0.0], + "EXT_FLDF_N_TTT": [0.0, 0.0] + )json", \ + R"json( "EXT_FLDF_KKIN": [0.0, 0.0, 0.0], + "EXT_FLDF_KKIN_T": [1.0, 0.5, 2.0], + "EXT_FLDF_KKIN_TT": [0.0, 0.0, 0.0], + "EXT_FLDF_KKIN_TTT": [0.0, 0.0, 0.0], + "EXT_FLDF_KF": [0.0, 0.0, 0.0], + "EXT_FLDF_KF_T": [0.1, 0.3, 0.2], + "EXT_FLDF_KF_TT": [0.0, 0.0, 0.0], + "EXT_FLDF_KF_TTT": [0.0, 0.0, 0.0], + "EXT_FLDF_N": [0.0, 0.0, 0.0], + "EXT_FLDF_N_T": [0.5, 2.0, 1.2], + "EXT_FLDF_N_TT": [0.0, 0.0, 0.0], + "EXT_FLDF_N_TTT": [0.0, 0.0, 0.0] + )json", \ + 1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_UNUSED, CADET_COMPARE_BINDING_VS_NONBINDING) CADET_BINDINGTEST("MULTI_COMPONENT_LANGMUIR", "EXT_MULTI_COMPONENT_LANGMUIR", (1,1), (1,0,1), (1.0, 2.0, 0.0, 0.0), (1.0, 3.0, 2.0, 0.0, 0.0), \ R"json( "MCL_KA": [1.14, 2.0], @@ -80,6 +116,79 @@ CADET_BINDINGTEST("MULTI_COMPONENT_LANGMUIR", "EXT_MULTI_COMPONENT_LANGMUIR", (1 )json", \ 1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_UNUSED, CADET_COMPARE_BINDING_VS_NONBINDING) + CADET_BINDINGTEST("MULTI_COMPONENT_LANGMUIR_LDF", "EXT_MULTI_COMPONENT_LANGMUIR_LDF", (1, 1), (1, 0, 1), (1.0, 2.0, 0.0, 0.0), (1.0, 3.0, 2.0, 0.0, 0.0), \ + R"json( "MCLLDF_KEQ": [1.14, 2.0], + "MCLLDF_KKIN": [0.004, 0.008], + "MCLLDF_QMAX": [4.88, 3.5] + )json", \ + R"json( "MCLLDF_KEQ": [1.14, 1.0, 2.0], + "MCLLDF_KKIN": [0.004, 2.0, 0.008], + "MCLLDF_QMAX": [4.88, 3.0, 3.5] + )json", \ + R"json( "EXT_MCLLDF_KEQ": [0.0, 0.0], + "EXT_MCLLDF_KEQ_T": [1.14, 2.0], + "EXT_MCLLDF_KEQ_TT": [0.0, 0.0], + "EXT_MCLLDF_KEQ_TTT": [0.0, 0.0], + "EXT_MCLLDF_KKIN": [0.0, 0.0], + "EXT_MCLLDF_KKIN_T": [0.004, 0.008], + "EXT_MCLLDF_KKIN_TT": [0.0, 0.0], + "EXT_MCLLDF_KKIN_TTT": [0.0, 0.0], + "EXT_MCLLDF_QMAX": [0.0, 0.0], + "EXT_MCLLDF_QMAX_T": [4.88, 3.5], + "EXT_MCLLDF_QMAX_TT": [0.0, 0.0], + "EXT_MCLLDF_QMAX_TTT": [0.0, 0.0] + )json", \ + R"json( "EXT_MCLLDF_KEQ": [0.0, 0.0, 0.0], + "EXT_MCLLDF_KEQ_T": [1.14, 1.0, 2.0], + "EXT_MCLLDF_KEQ_TT": [0.0, 0.0, 0.0], + "EXT_MCLLDF_KEQ_TTT": [0.0, 0.0, 0.0], + "EXT_MCLLDF_KKIN": [0.0, 0.0, 0.0], + "EXT_MCLLDF_KKIN_T": [0.004, 2.0, 0.008], + "EXT_MCLLDF_KKIN_TT": [0.0, 0.0, 0.0], + "EXT_MCLLDF_KKIN_TTT": [0.0, 0.0, 0.0], + "EXT_MCLLDF_QMAX": [0.0, 0.0, 0.0], + "EXT_MCLLDF_QMAX_T": [4.88, 3.0, 3.5], + "EXT_MCLLDF_QMAX_TT": [0.0, 0.0, 0.0], + "EXT_MCLLDF_QMAX_TTT": [0.0, 0.0, 0.0] + )json", \ + 1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_UNUSED, CADET_COMPARE_BINDING_VS_NONBINDING) + + CADET_BINDINGTEST("MULTI_COMPONENT_LANGMUIR_LDF_LIQUID_PHASE", "EXT_MULTI_COMPONENT_LANGMUIR_LDF_LIQUID_PHASE", (1, 1), (1, 0, 1), (1.0, 2.0, 0.0, 0.0), (1.0, 3.0, 2.0, 0.0, 0.0), \ + R"json( "MCLLDFC_KEQ": [1.14, 2.0], + "MCLLDFC_KKIN": [0.004, 0.008], + "MCLLDFC_QMAX": [4.88, 3.5] + )json", \ + R"json( "MCLLDFC_KEQ": [1.14, 1.0, 2.0], + "MCLLDFC_KKIN": [0.004, 2.0, 0.008], + "MCLLDFC_QMAX": [4.88, 3.0, 3.5] + )json", \ + R"json( "EXT_MCLLDFC_KEQ": [0.0, 0.0], + "EXT_MCLLDFC_KEQ_T": [1.14, 2.0], + "EXT_MCLLDFC_KEQ_TT": [0.0, 0.0], + "EXT_MCLLDFC_KEQ_TTT": [0.0, 0.0], + "EXT_MCLLDFC_KKIN": [0.0, 0.0], + "EXT_MCLLDFC_KKIN_T": [0.004, 0.008], + "EXT_MCLLDFC_KKIN_TT": [0.0, 0.0], + "EXT_MCLLDFC_KKIN_TTT": [0.0, 0.0], + "EXT_MCLLDFC_QMAX": [0.0, 0.0], + "EXT_MCLLDFC_QMAX_T": [4.88, 3.5], + "EXT_MCLLDFC_QMAX_TT": [0.0, 0.0], + "EXT_MCLLDFC_QMAX_TTT": [0.0, 0.0] + )json", \ + R"json( "EXT_MCLLDFC_KEQ": [0.0, 0.0, 0.0], + "EXT_MCLLDFC_KEQ_T": [1.14, 1.0, 2.0], + "EXT_MCLLDFC_KEQ_TT": [0.0, 0.0, 0.0], + "EXT_MCLLDFC_KEQ_TTT": [0.0, 0.0, 0.0], + "EXT_MCLLDFC_KKIN": [0.0, 0.0, 0.0], + "EXT_MCLLDFC_KKIN_T": [0.004, 2.0, 0.008], + "EXT_MCLLDFC_KKIN_TT": [0.0, 0.0, 0.0], + "EXT_MCLLDFC_KKIN_TTT": [0.0, 0.0, 0.0], + "EXT_MCLLDFC_QMAX": [0.0, 0.0, 0.0], + "EXT_MCLLDFC_QMAX_T": [4.88, 3.0, 3.5], + "EXT_MCLLDFC_QMAX_TT": [0.0, 0.0, 0.0], + "EXT_MCLLDFC_QMAX_TTT": [0.0, 0.0, 0.0] + )json", \ + 1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_UNUSED, CADET_COMPARE_BINDING_VS_NONBINDING) CADET_BINDINGTEST("MULTI_COMPONENT_ANTILANGMUIR", "EXT_MULTI_COMPONENT_ANTILANGMUIR", (1,1), (1,0,1), (1.0, 2.0, 0.0, 0.0), (1.0, 3.0, 2.0, 0.0, 0.0), \ R"json( "MCAL_KA": [1.14, 2.0], @@ -586,7 +695,25 @@ CADET_BINDINGTEST_ALLBINDING("MULTI_COMPONENT_BILANGMUIR", "EXT_MULTI_COMPONENT_ "EXT_MCBL_QMAX_TTT": [0.0, 0.0, 0.0, 0.0] )json", \ 1e-10, 1e-10) - +CADET_BINDINGTEST_ALLBINDING("MULTI_COMPONENT_BILANGMUIR_LDF", "EXT_MULTI_COMPONENT_BILANGMUIR_LDF", (2, 2), (1.0, 2.0, 0.0, 0.0, 0.0, 0.0), \ + R"json( "MCBLLDF_KEQ": [1.14, 2.0, 2.28, 4.0], + "MCBLLDF_KKIN": [0.004, 0.008, 0.002, 0.003], + "MCBLLDF_QMAX": [4.88, 3.5, 3.88, 2.5] + )json", \ + R"json( "EXT_MCBLLDF_KEQ": [0.0, 0.0, 0.0, 0.0], + "EXT_MCBLLDF_KEQ_T": [1.14, 2.0, 2.28, 4.0], + "EXT_MCBLLDF_KEQ_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_MCBLLDF_KEQ_TTT": [0.0, 0.0, 0.0, 0.0], + "EXT_MCBLLDF_KKIN": [0.0, 0.0, 0.0, 0.0], + "EXT_MCBLLDF_KKIN_T": [0.004, 0.008, 0.002, 0.003], + "EXT_MCBLLDF_KKIN_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_MCBLLDF_KKIN_TTT": [0.0, 0.0, 0.0, 0.0], + "EXT_MCBLLDF_QMAX": [0.0, 0.0, 0.0, 0.0], + "EXT_MCBLLDF_QMAX_T": [4.88, 3.5, 3.88, 2.5], + "EXT_MCBLLDF_QMAX_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_MCBLLDF_QMAX_TTT": [0.0, 0.0, 0.0, 0.0] + )json", \ + 1e-10, 1e-10) CADET_BINDINGTEST("STERIC_MASS_ACTION", "EXT_STERIC_MASS_ACTION", (1,1,1), (1,1,0,1), (1.2, 2.0, 1.5, 80.0, 3.5, 2.7), (1.2, 2.0, 3.0, 1.5, 80.0, 3.5, 2.7), \ R"json( "SMA_KA": [0.0, 3.55, 1.59],