diff --git a/doc/interface/reaction/activated_sludge_model.rst b/doc/interface/reaction/activated_sludge_model.rst new file mode 100644 index 000000000..183e4ed47 --- /dev/null +++ b/doc/interface/reaction/activated_sludge_model.rst @@ -0,0 +1,276 @@ +.. _activated_sludge_model_config: + +Activated Sludge Model (ASM3h) +~~~~~~~~~~~~~~~~~~~~~~~~~ + +**Group /input/model/unit_XXX/reaction - REACTION_MODEL = ACTIVATED_SLUDGE_MODEL3** + +For information on model equations, refer to :ref:`activated_sludge_model`. + +**Component configuration** + +``ASM3_COMP_IDX`` + + Optional component indexes. Set this in case the relevant components start at a certain offset or are provided + in a different order than listed below: + + ============= ========= + **Component** **Index** + ============= ========= + SO 0 + SS 1 + SNH 2 + SNO 3 + SN2 4 + SALK 5 + SI 6 + XI 7 + XS 8 + XH 9 + XSTO 10 + XA 11 + XMI 12 + ============= ========= + + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{N}` **Length:** 13 + ================ ============================= ======================================================== + +**Environment/process parameters** + +``ASM3_T`` + + Temperature :math:`T`. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_V`` + + Reactor volume :math:`V`. Set this to the column volume since the model can't compute it itself. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\gt 0` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_IO2`` + + Aeration oxygen input :math:`iO_2`. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\ge 0` **Length:** 1 + ================ ============================= ======================================================== + +**Maximum rates** + +``ASM3_KH20`` + + Hydrolysis rate constant :math:`k_H` at 20 °C. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_KSTO20`` + + Maximum storage rate :math:`K_{STO}` at 20 °C. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\gt 0` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_MU_H20`` + + Heterotrophic max. growth rate :math:`\mu_{H}` at 20 °C. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_BH20`` + + Rate constant for lysis and decay :math:`b_H` at 20 °C. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_MU_AUT20`` + + Autotrophic max. growth rate :math:`\mu_{AUT}` at 20 °C. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_BAUT20`` + + Rate constant :math:`b_{AUT}` for decay of autotrophs at 20 °C. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +**Anoxic reduction factors** + +``ASM3_ETA_HNO3`` + + Anoxic reduction factor :math:`\eta_{HNO_3}` – heterotrophic growth. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_ETAH_END`` + + Anoxic reduction factor :math:`\eta_{H_{end}}` – endogenous respiration XH. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_ETAN_END`` + + Anoxic reduction factor :math:`\eta_{N_{end}}` – endogenous respiration XA. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +**Saturation/inhibition coefficients** + +``ASM3_KX`` + + Saturation/inhibition coefficient :math:`KX` for particulate COD. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_KHO2`` + + Saturation/inhibition coefficient :math:`KHO_2` for oxygen, heterotrophic growth. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_KHSS`` + + Saturation/inhibition coefficient :math:`KHSS` for readily biodegradable substrates. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_KHNO3`` + + Saturation/inhibition coefficient :math:`KHNO_3` for nitrate. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_KHNH4`` + + Saturation/inhibition coefficient :math:`KHNH_4` for ammonium (nutrient). + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_KHALK`` + + Saturation coefficient :math:`KH_{ALK}` for alkalinity (HCO3-). + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_KHSTO`` + + Saturation coefficient :math:`KH_{STO}` for storage products. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_KNO2`` + + Saturation coefficient :math:`K_{NO_2}` for oxygen, autotrophic growth. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_KNNH4`` + + Saturation coefficient :math:`K_{NNH_4}` for ammonium (substrate), autotrophic growth. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +``ASM3_KNALK`` + + Saturation coefficient :math:`K_{NALK}` for alkalinity (HCO3-), autotrophic growth. + + ================ ============================= ======================================================== + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ======================================================== + +**Example** + +.. python:: + + # Example of setting up an ASM3 reaction model in a unit operation with bulk reaction + + # Setup ASM3 reaction for unit 000 with example values + model.root.input.model.unit_000.reaction_model = 'ACTIVATED_SLUDGE_MODEL3' + model.root.input.model.unit_000.reaction_bulk.asm3_insi = 0.01 + model.root.input.model.unit_000.reaction_bulk.asm3_inss = 0.03 + model.root.input.model.unit_000.reaction_bulk.asm3_inxi = 0.04 + model.root.input.model.unit_000.reaction_bulk.asm3_inxs = 0.03 + model.root.input.model.unit_000.reaction_bulk.asm3_inbm = 0.07 + model.root.input.model.unit_000.reaction_bulk.asm3_ivss_xi = 0.751879699 + model.root.input.model.unit_000.reaction_bulk.asm3_ivss_xs = 0.555555556 + model.root.input.model.unit_000.reaction_bulk.asm3_ivss_sto = 0.6 + model.root.input.model.unit_000.reaction_bulk.asm3_ivss_bm = 0.704225352 + model.root.input.model.unit_000.reaction_bulk.asm3_itss_vss_bm = 1.086956522 + + + model.root.input.model.unit_000.reaction_bulk.asm3_fiss_bm_prod = 1 + model.root.input.model.unit_000.reaction_bulk.asm3_fsi = 0 + model.root.input.model.unit_000.reaction_bulk.asm3_yh_aer = 0.8 + model.root.input.model.unit_000.reaction_bulk.asm3_yh_anox = 0.65 + + model.root.input.model.unit_000.reaction_bulk.asm3_ysto_aer = 0.8375 + model.root.input.model.unit_000.reaction_bulk.asm3_ysto_anox = 0.7 + model.root.input.model.unit_000.reaction_bulk.asm3_fxi = 0.2 + model.root.input.model.unit_000.reaction_bulk.asm3_ya = 0.24 + model.root.input.model.unit_000.reaction_bulk.asm3_kh20 = 9 + model.root.input.model.unit_000.reaction_bulk.asm3_kx = 1 + model.root.input.model.unit_000.reaction_bulk.asm3_ksto20 = 12 + model.root.input.model.unit_000.reaction_bulk.asm3_mu_h20 = 3 + model.root.input.model.unit_000.reaction_bulk.asm3_bh20 = 0.33 + model.root.input.model.unit_000.reaction_bulk.asm3_eta_hno3 = 0.5 + model.root.input.model.unit_000.reaction_bulk.asm3_khO2 = 0.2 + model.root.input.model.unit_000.reaction_bulk.asm3_khss = 10 + model.root.input.model.unit_000.reaction_bulk.asm3_khno3 = 0.5 + model.root.input.model.unit_000.reaction_bulk.asm3_khnh4 = 0.01 + model.root.input.model.unit_000.reaction_bulk.asm3_khalk = 0.1 + model.root.input.model.unit_000.reaction_bulk.asm3_khsto = 0.1 + model.root.input.model.unit_000.reaction_bulk.asm3_mu_aut20 = 1.12 + model.root.input.model.unit_000.reaction_bulk.asm3_baut20 = 0.18 + model.root.input.model.unit_000.reaction_bulk.asm3_etah_end = 0.5 + model.root.input.model.unit_000.reaction_bulk.asm3_etan_end = 0.5 + model.root.input.model.unit_000.reaction_bulk.asm3_kno2 = 0.5 + model.root.input.model.unit_000.reaction_bulk.asm3_knnh4 = 0.7 + model.root.input.model.unit_000.reaction_bulk.asm3_knalk = 0.5 + model.root.input.model.unit_000.reaction_bulk.asm3_t = 12 + + + model.root.input.model.unit_000.reaction_bulk.asm3_v = 1000.0 + model.root.input.model.unit_000.reaction_bulk.asm3_io2 = 0.0 diff --git a/doc/interface/reaction/index.rst b/doc/interface/reaction/index.rst index a1b6e14a7..e8e955750 100644 --- a/doc/interface/reaction/index.rst +++ b/doc/interface/reaction/index.rst @@ -8,6 +8,7 @@ Reaction models mass_action_law michaelis_menten_kinetics + activated_sludge_model Externally dependent reaction models ------------------------------------ diff --git a/doc/modelling/reaction/activated_sludge_model.rst b/doc/modelling/reaction/activated_sludge_model.rst new file mode 100644 index 000000000..594c6b83a --- /dev/null +++ b/doc/modelling/reaction/activated_sludge_model.rst @@ -0,0 +1,336 @@ +.. _activated_sludge_model: + +Activated Sludge Model (ASM3h) +=============================== + +Activated Sludge Models (ASM) are a group of reaction models published by the International Water Association (IWA) +and mostly used for modelling biological processes in wastewater treatment plants. + +It consists of 13 components representing different species in wastewater treatment plants. + +The used components with their abbreviations are listed in the following table: + +.. list-table:: + :header-rows: 1 + :widths: 20 80 + + * - Abbreviation + - Component + * - :math:`SO` + - Dissolved oxygen + * - :math:`SS` + - Readily biodegradable substrate + * - :math:`SNH` + - Ammonium + * - :math:`SNO` + - Nitrite and nitrate + * - :math:`SN2` + - Dinitrogen, released by denitrification + * - :math:`SALK` + - Alkalinity, bicarbonate + * - :math:`SI` + - Soluble inert organic + * - :math:`XI` + - Inert Particulate organic + * - :math:`XS` + - Slowly biodegradable substrate + * - :math:`XH` + - Heterotrophic biomass + * - :math:`XSTO` + - Organics stored by heterotrophs + * - :math:`XA` + - Autotrophic, nitrifying biomass + * - :math:`XMI` + - Mineral particulate matter from biomass + +The net flux for each of the components in the ASM is calculated with a stoichiometric +matrix :math:`S \in \mathbb{R}^{13 \times 13}` and a reaction rate vector +:math:`\varphi_j(c)`, where :math:`j` is the reaction number. +The list of reactions included in the ASM3h model is as follows: + +Reaction Rates +~~~~~~~~~~~~~~ + +.. list-table:: + :header-rows: 1 + :widths: 15 35 50 + + * - Reaction Number + - Reaction Description + - Rate Equation + * - :math:`r_1` + - Hydrolysis of organic structures + - :math:`k_{h,20} \cdot f_T^{0.04} \cdot \frac{XS/XH}{XS/XH + K_X} \cdot XH` + * - :math:`r_2` + - Aerobic storage of SS + - :math:`k_{STO} \cdot \frac{SO}{SO + K_{h,O2}} \cdot \frac{SS}{SS + K_{h,SS}} \cdot XH` + * - :math:`r_3` + - Anoxic storage of SS + - :math:`k_{STO} \cdot \eta_{h,NO3} \cdot \frac{K_{h,O2}}{SO + K_{h,O2}} \cdot \frac{SNO}{SNO + K_{h,NO3}} \cdot \frac{SS}{SS + K_{h,SS}} \cdot XH` + * - :math:`r_4` + - Aerobic growth of XH + - :math:`\mu_H \cdot \frac{SO}{SO + K_{h,O2}} \cdot \frac{XSTO/XH}{XSTO/XH + K_{h,STO}} \cdot \frac{SNH}{SNH + K_{h,NH4}} \cdot \frac{SALK}{SALK + K_{h,ALK}} \cdot XH` + * - :math:`r_5` + - Anoxic growth of XH (denitrification) + - :math:`\mu_H \cdot \eta_{h,NO3} \cdot \frac{K_{h,O2}}{SO + K_{h,O2}} \cdot \frac{SNO}{SNO + K_{h,NO3}} \cdot \frac{XSTO/XH}{XSTO/XH + K_{h,STO}} \cdot \frac{SNH}{SNH + K_{h,NH4}} \cdot \frac{SALK}{SALK + K_{h,ALK}} \cdot XH` + * - :math:`r_6` + - Aerobic endogenous respiration of XH + - :math:`b_H \cdot \frac{SO}{SO + K_{h,O2}} \cdot XH` + * - :math:`r_7` + - Anoxic endogenous respiration of XH + - :math:`b_H \cdot \eta_{h,end} \cdot \frac{K_{h,O2}}{SO + K_{h,O2}} \cdot \frac{SNO}{SNO + K_{h,NO3}} \cdot XH` + * - :math:`r_8` + - Aerobic respiration of internal cell storage products + - :math:`b_H \cdot \frac{SO}{SO + K_{h,O2}} \cdot \frac{XSTO}{XSTO + K_{h,STO}} \cdot XH` + * - :math:`r_9` + - Anoxic respiration of internal cell storage products + - :math:`b_H \cdot \eta_{h,end} \cdot \frac{K_{h,O2}}{SO + K_{h,O2}} \cdot \frac{SNO}{SNO + K_{h,NO3}} \cdot \frac{XSTO}{XSTO + K_{h,STO}} \cdot XH` + * - :math:`r_{10}` + - Aerobic growth of XA + - :math:`\mu_{AUT} \cdot \frac{SO}{SO + K_{N,O2}} \cdot \frac{SNH}{SNH + K_{N,NH4}} \cdot \frac{SALK}{SALK + K_{N,ALK}} \cdot XA` + * - :math:`r_{11}` + - Aerobic endogenous respiration of XA + - :math:`b_{AUT} \cdot \frac{SO}{SO + K_{N,O2}} \cdot XA` + * - :math:`r_{12}` + - Anoxic endogenous respiration of XA + - :math:`b_{AUT} \cdot \eta_{N,end} \cdot \frac{K_{N,O2}}{SO + K_{N,O2}} \cdot \frac{SNO}{SNO + K_{h,NO3}} \cdot XA` + * - :math:`r_{13}` + - Aeration + - :math:`\frac{iO_2}{1000}` + + +And the stoichiometric equations are represented in the stoichiometric matrix :math:`S`, which is given by: + +.. list-table:: + :header-rows: 1 + :widths: 10 7 7 7 7 7 7 7 7 7 7 7 7 7 + + * - Component\\Reaction + - :math:`r_1` + - :math:`r_2` + - :math:`r_3` + - :math:`r_4` + - :math:`r_5` + - :math:`r_6` + - :math:`r_7` + - :math:`r_8` + - :math:`r_9` + - :math:`r_{10}` + - :math:`r_{11}` + - :math:`r_{12}` + - :math:`r_{13}` + * - :math:`SO` + - 0 + - :math:`Y_{STO,aer} - 1` + - :math:`Y_{STO,anox} - 1` + - :math:`1 - 1/Y_{H,aer}` + - :math:`1 - 1/Y_{H,anox}` + - :math:` f_{XI} - 1` + - :math:`f_{XI} - 1` + - -1 + - -1 + - :math:`-(64/14) \cdot 1/Y_A + 1` + - :math:`f_{XI} - 1` + - :math:`f_{XI} - 1` + - 1 + * - :math:`SS` + - :math:`1 - f_{SI}` + - -1 + - -1 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + * - :math:`SNH` + - :math:`i_{N,XS} - i_{N,SI} \cdot f_{SI} - (1 - f_{SI}) \cdot i_{N,SS}` + - :math:`i_{N,SS}` + - :math:`i_{N,SS}` + - :math:`-i_{N,BM}` + - :math:`-i_{N,BM}` + - :math:`i_{N,BM} - f_{XI} \cdot i_{N,XI}` + - :math:`i_{N,BM} -f_{XI} \cdot i_{N,XI}` + - 0 + - 0 + - :math:`-1/Y_A - i_{N,BM}` + - :math:`f_{XI} \cdot i_{N,XI} + i_{N,BM}` + - :math:`f_{XI} \cdot i_{N,XI} + i_{N,BM}` + - 0 + * - :math:`SNO` + - 0 + - 0 + - :math:`(Y_{STO,anox} - 1) / (40/14)` + - 0 + - :math:`(1 - 1/Y_{H,anox}) / (40/14)` + - 0 + - :math:`(f_{XI} - 1) / (40/14)` + - 0 + - :math:`-14/40` + - :math:`1/Y_A` + - 0 + - :math:`(f_{XI} - 1) / (40/14)` + - 0 + * - :math:`SN2` + - 0 + - 0 + - :math:`(1 - Y_{STO,anox}) / (40/14)` + - 0 + - :math:`(1/Y_{H,anox} - 1) / (40/14)` + - 0 + - :math:`(-f_{XI} - 1) / (40/14)` + - 0 + - :math:`--14/40` + - 0 + - 0 + - :math:`(1 + f_{XI} ) / (40/14)` + - 0 + * - :math:`SALK` + - :math:`(i_{N,XS} - i_{N,SI} \cdot f_{SI} - (1 - f_{SI}) \cdot i_{N,SS})/14` + - :math:`i_{N,SS} / 14` + - :math:`i_{N,SS} / 14 - (Y_{STO,anox} - 1) / 40` + - :math:`-i_{N,BM} / 14` + - :math:`-i_{N,BM} / 14 - (1 - 1/Y_{H,anox}) / 40` + - :math:`(-f_{XI} \cdot i_{N,XI} + i_{N,BM}) / 14` + - :math:`(-f_{XI} \cdot i_{N,XI} + i_{N,BM}) / 14 (-f_{XI} + 1) / 40` + - 0 + - :math:`1/40` + - :math:`-(1/Y_A) 1/7- i_{N,BM}/14` + - :math:`f_{XI} \cdot i_{N,XI} + i_{N,BM} / 14` + - :math:`(-f_{XI} \cdot i_{N,XI} + i_{N,BM})/14 - (f_{XI} - 1) / 40` + - 0 + * - :math:`SI` + - :math:`f_{SI}` + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + * - :math:`XI` + - 0 + - 0 + - 0 + - 0 + - 0 + - :math:`f_{XI}` + - :math:`f_{XI}` + - 0 + - 0 + - 0 + - :math:`f_{XI}` + - :math:`f_{XI}` + - 0 + * - :math:`XS` + - -1 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + * - :math:`XH` + - 0 + - 0 + - 0 + - 1 + - 1 + - -1 + - -1 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + * - :math:`XSTO` + - 0 + - :math:`Y_{STO,aer}` + - :math:`Y_{STO,anox}` + - :math:`-1/Y_{H,aer}` + - :math:`-1/Y_{H,anox}` + - 0 + - 0 + - -1 + - -1 + - 0 + - 0 + - 0 + - 0 + * - :math:`XA` + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 0 + - 1 + - -1 + - -1 + - 0 + * - :math:`XMI` + - 0 + - 0 + - 0 + - 0 + - 0 + - :math:`f_{XMI,BI}` + - :math:`f_{XMI,BI}` + - 0 + - 0 + - 0 + - :math:`f_{XMI,BI}` + - :math:`f_{XMI,BI}` + - 0 + + +.. figure:: images/netzwerk.jpg + :width: 80% + :align: center + + Network representation of the AMS3h model. + +For more information on model parameters required to define in CADET file format, see :ref:`activated_sludge_model_config`. + +Combining the ASM3 model within CADET +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +The ASM3h model is a konkrete implementation of of the ASM framework which is aimled to be used as ASM3h in SIMBA. +But with different configuations it is also flexible to to be used in the general CADET framework. + +Active Aeration +--------------- +In configuation of the ASM3 model, there is the option to set the volume parameter of the aeration reaction. +This volume refers to the volume of the aeration tank. +But we recomend to model the aeration by setting up the Inlet unit Operation and connect it to the unit operation where the ASM3 model is used. +This way the aeration can be handled more flexible, see :ref:`inlet_operation`. +To deactivate the aeration reaction in the ASM3 model, set the volume parameter to zero. + +Fractionation +------------- + +TODO + + + +References +~~~~~~~~~~ +- TODO \ No newline at end of file diff --git a/doc/modelling/reaction/images/netzwerk.jpg b/doc/modelling/reaction/images/netzwerk.jpg new file mode 100644 index 000000000..fae20be83 Binary files /dev/null and b/doc/modelling/reaction/images/netzwerk.jpg differ diff --git a/doc/modelling/reaction/index.rst b/doc/modelling/reaction/index.rst index 58b8991b8..5c539fbb5 100644 --- a/doc/modelling/reaction/index.rst +++ b/doc/modelling/reaction/index.rst @@ -17,6 +17,10 @@ Historically, a chromatography system is modeled as a reaction system without co - :ref:`thomas_model` - :ref:`rate_constant_distribution_theory` +More specialized models are also implemented: + + - :ref:`activated_sludge_model` + .. _dependence-on-external-function_react: Dependence on external function diff --git a/src/libcadet/CMakeLists.txt b/src/libcadet/CMakeLists.txt index 0d2a1746b..d10286f5e 100644 --- a/src/libcadet/CMakeLists.txt +++ b/src/libcadet/CMakeLists.txt @@ -116,7 +116,7 @@ set(LIBCADET_REACTIONMODEL_SOURCES ${CMAKE_SOURCE_DIR}/src/libcadet/model/reaction/DummyReaction.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/reaction/MassActionLawReaction.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/reaction/MichaelisMentenReaction.cpp - + ${CMAKE_SOURCE_DIR}/src/libcadet/model/reaction/ActivatedSludgeModelThree.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/reaction/CrystallizationReaction.cpp ) diff --git a/src/libcadet/ReactionModelFactory.cpp b/src/libcadet/ReactionModelFactory.cpp index 598662c47..3976b4cb8 100644 --- a/src/libcadet/ReactionModelFactory.cpp +++ b/src/libcadet/ReactionModelFactory.cpp @@ -22,6 +22,7 @@ namespace cadet void registerMassActionLawReaction(std::unordered_map>& reactions); void registerDummyReaction(std::unordered_map>& reactions); void registerMichaelisMentenReaction(std::unordered_map>& reactions); + void registerActivatedSludgeModelThreeReaction(std::unordered_map>& reactions); void registerCrystallizationReaction(std::unordered_map>& reactions); } @@ -33,6 +34,7 @@ namespace cadet model::reaction::registerDummyReaction(_dynamicModels); model::reaction::registerMassActionLawReaction(_dynamicModels); model::reaction::registerMichaelisMentenReaction(_dynamicModels); + model::reaction::registerActivatedSludgeModelThreeReaction(_dynamicModels); model::reaction::registerCrystallizationReaction(_dynamicModels); } diff --git a/src/libcadet/model/reaction/ActivatedSludgeModelThree.cpp b/src/libcadet/model/reaction/ActivatedSludgeModelThree.cpp new file mode 100644 index 000000000..6309310e5 --- /dev/null +++ b/src/libcadet/model/reaction/ActivatedSludgeModelThree.cpp @@ -0,0 +1,925 @@ +// ============================================================================= +// CADET - The Chromatography Analysis and Design Toolkit +// +// Copyright © 2008-present: The CADET-Core Authors +// Please see the AUTHORS.md 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/reaction/ReactionModelBase.hpp" +#include "model/ExternalFunctionSupport.hpp" +#include "model/ModelUtils.hpp" +#include "cadet/Exceptions.hpp" +#include "model/Parameters.hpp" +#include "LocalVector.hpp" +#include "SimulationTypes.hpp" +#include "linalg/ActiveDenseMatrix.hpp" +#include "MathUtil.hpp" +#include "Memory.hpp" + +#include +#include +#include +#include + +/* +{ + "name": "ActivatedSludgeModelThreeParamHandler", + "externalName": "ExtActivatedSludgeModelThreeParamHandler", + "parameters": + [ + { "type": "ScalarParameter", "varName": "kh20", "confName": "ASM3_KH20"}, + { "type": "ScalarParameter", "varName": "T", "confName": "ASM3_T"}, + { "type": "ScalarParameter", "varName": "io2", "confName": "ASM3_IO2"}, + { "type": "ScalarParameter", "varName": "V", "confName": "ASM3_V"}, + { "type": "ScalarParameter", "varName": "k_sto20", "confName": "ASM3_KSTO20"}, + { "type": "ScalarParameter", "varName": "kx", "confName": "ASM3_KX"}, + { "type": "ScalarParameter", "varName": "kho2", "confName": "ASM3_KHO2"}, + { "type": "ScalarParameter", "varName": "khss", "confName": "ASM3_KHSS"}, + { "type": "ScalarParameter", "varName": "khn03", "confName": "ASM3_KHNO3"}, + { "type": "ScalarParameter", "varName": "etahno3", "confName": "ASM3_ETA_HNO3"}, + { "type": "ScalarParameter", "varName": "khnh4", "confName": "ASM3_KHNH4"}, + { "type": "ScalarParameter", "varName": "khalk", "confName": "ASM3_KHALK"}, + { "type": "ScalarParameter", "varName": "khsto", "confName": "ASM3_KHSTO"}, + { "type": "ScalarParameter", "varName": "muh2o", "confName": "ASM3_MU_H20"}, + { "type": "ScalarParameter", "varName": "etahend", "confName": "ASM3_ETAH_END"}, + { "type": "ScalarParameter", "varName": "bh20", "confName": "ASM3_BH20"}, + { "type": "ScalarParameter", "varName": "muAUT20", "confName": "ASM3_MU_AUT20"}, + { "type": "ScalarParameter", "varName": "kno2", "confName": "ASM3_KNO2"}, + { "type": "ScalarParameter", "varName": "knnh4", "confName": "ASM3_KNNH4"}, + { "type": "ScalarParameter", "varName": "knalk", "confName": "ASM3_KNALK"}, + { "type": "ScalarParameter", "varName": "baut20", "confName": "ASM3_BAUT20"}, + { "type": "ScalarParameter", "varName": "etanend", "confName": "ASM3_ETAN_END"} + ] +} +*/ + +/* Parameter description + ------------------------ + + +*/ + + +namespace cadet +{ + +namespace model +{ + +inline const char* ActivatedSludgeModelThreeParamHandler::identifier() CADET_NOEXCEPT { return "ACTIVATED_SLUDGE_MODEL3"; } + +inline bool ActivatedSludgeModelThreeParamHandler::validateConfig(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates) +{ + return true; +} + +inline const char* ExtActivatedSludgeModelThreeParamHandler::identifier() CADET_NOEXCEPT { return "EXT_ACTIVATED_SLUDGE_MODEL3"; } + +inline bool ExtActivatedSludgeModelThreeParamHandler::validateConfig(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates) +{ + return true; +} + +namespace +{ + /** + * @brief Registers a matrix-valued parameter (row-major storage) with components as rows + * @details The matrix-valued parameter has as many rows as there are components in the system. + * @param [in,out] parameters Parameter map + * @param [in] unitOpIdx Unit operation id + * @param [in] parTypeIdx Particle type index + * @param [in] paramName Name of the parameter + * @param [in] mat Matrix to register + */ + inline void registerCompRowMatrix(std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, const std::string& paramName, cadet::linalg::ActiveDenseMatrix& mat) + { + const cadet::StringHash hashName = cadet::hashStringRuntime(paramName); + cadet::registerParam2DArray(parameters, mat.data(), mat.elements(), [=](bool multi, unsigned int row, unsigned int col) + { + return cadet::makeParamId(hashName, unitOpIdx, row, parTypeIdx, cadet::BoundStateIndep, col, cadet::SectionIndep); + }, + mat.columns() + ); + } +} + + +/** + * @brief Defines a Michaelis-Menten reaction kinetic with simple inhibition + * @details Implements the Michaelis-Menten kinetics: \f[ \begin{align} + * S \nu, + * \end{align} \f] + * where \f$ S \f$ is the stoichiometric matrix and the fluxes are given by + * \f[ \begin{align} + * \nu_i = \frac{\mu_{\mathrm{max},i} c_S}{k_{\mathrm{MM},i} + c_S}. + * \end{align} \f] + * The substrate component \f$ c_S \f$ is identified by the index of the + * first negative entry in the stoichiometry of this reaction. + * In addition, the reaction might be inhibited by other components. In this + * case, the flux has the form + * \f[ \begin{align} + * \nu_i = \frac{\mu_{\mathrm{max},i} c_S}{k_{\mathrm{MM},i} + c_S} \cdot \frac{1}{1 + \sum_i \frac{1+ k_{\mathrm{I},i,j}}{k_{\mathrm{I},i,j} + c_{\mathrm{I},j}}}. + * \end{align} \f] + * The value of \f$ k_{\mathrm{I},i,j} \f$ decides whether component \f$ j \f$ + * inhibits reaction \f$ i \f$. If \f$ k_{\mathrm{I},i,j} \leq 0 \f$, the component + * does not inhibit the reaction. + * Only reactions in liquid phase are supported (no solid phase or cross-phase reactions). + * @tparam ParamHandler_t Type that can add support for external function dependence + */ +template +class ActivatedSludgeModelThreeBase : public DynamicReactionModelBase +{ +public: + + ActivatedSludgeModelThreeBase() CADET_NOEXCEPT { } + virtual ~ActivatedSludgeModelThreeBase() CADET_NOEXCEPT { } + + static const char* identifier() { return ParamHandler_t::identifier(); } + virtual const char* name() const CADET_NOEXCEPT { return ParamHandler_t::identifier(); } + + virtual void setExternalFunctions(IExternalFunction** extFuns, unsigned int size) { _paramHandler.setExternalFunctions(extFuns, size); } + virtual bool dependsOnTime() const CADET_NOEXCEPT { return ParamHandler_t::dependsOnTime(); } + virtual bool requiresWorkspace() const CADET_NOEXCEPT { return true; } + virtual unsigned int workspaceSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nBoundStates) const CADET_NOEXCEPT + { + return _paramHandler.cacheSize(_stoichiometry.columns(), nComp, totalNumBoundStates) + std::max(_stoichiometry.columns() * sizeof(active), 2 * (_nComp + totalNumBoundStates) * sizeof(double)); + } + + virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBound, unsigned int const* boundOffset) + { + DynamicReactionModelBase::configureModelDiscretization(paramProvider, nComp, nBound, boundOffset); + _stoichiometry.resize(nComp, 13); + return true; + } + + virtual unsigned int numReactionsLiquid() const CADET_NOEXCEPT { return _stoichiometry.columns(); } + virtual unsigned int numReactionsCombined() const CADET_NOEXCEPT { return 0; } + + CADET_DYNAMICREACTIONMODEL_BOILERPLATE + +protected: + ParamHandler_t _paramHandler; //!< Handles parameters and their dependence on external functions + + linalg::ActiveDenseMatrix _stoichiometry; + + unsigned int _idxSO = 0; //!< SO component index, default 0 + unsigned int _idxSS_nad = 1; //!< SS / SS_ad component index, default 1 + unsigned int _idxSNH = 2; //!< SNH component index, default 2 + unsigned int _idxSNO = 3; //!< SNO component index, default 3 + unsigned int _idxSN2 = 4; //!< SN2 component index, default 4 + unsigned int _idxSALK = 5; //!< SALK component index, default 5 + unsigned int _idxSI_nad = 6; //!< SI component index, default 6 + unsigned int _idxXI = 7; //!< XI component index, default 7 + unsigned int _idxXS = 8; //!< XS component index, default 8 + unsigned int _idxXH = 9; //!< XH component index, default 9 + unsigned int _idxXSTO = 10; //!< XSTO component index, default 10 + unsigned int _idxXA = 11; //!< XA component index, default 11 + unsigned int _idxXMI = 12; //!< XMI component index, default 12 + unsigned int _idxSI_ad = 13; + unsigned int _idxSS_ad = 14; + + bool _fractionate = false; + bool _activeAeration = true; + + + virtual bool configureStoich(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) + { + return true; + } + + + virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) + { + _paramHandler.configure(paramProvider, _stoichiometry.columns(), _nComp, _nBoundStates); + _paramHandler.registerParameters(_parameters, unitOpIdx, parTypeIdx, _nComp, _nBoundStates); + + if (paramProvider.exists("ASM3_FRACTIONATE")) + _fractionate = paramProvider.getBool("ASM3_FRACTIONATE"); + + if (_nComp < 13 && !_fractionate) + throw InvalidParameterException("ASM3 configuration: To use the ASM3 model the number of components must be at least 13"); + else if (_nComp < 15 && _fractionate) + throw InvalidParameterException("ASM3 configuration: To use the ASM3 model with fractionation the number of components must be at least 15"); + + if (paramProvider.exists("ASM3_COMP_IDX")) + { + const std::vector compIdx = paramProvider.getUint64Array("ASM3_COMP_IDX"); + if (_fractionate) + { + if (compIdx.size() != 15) + throw InvalidParameterException("ASM3 configuration: ASM3_COMP_IDX must have 15 elements"); + + _idxSS_ad = compIdx[13]; + _idxSI_ad = compIdx[14]; + } + else if (compIdx.size() != 13) + { + throw InvalidParameterException("ASM3 configuration: ASM3_COMP_IDX must have 13 elements"); + } + + LOG(Debug) << "ASM3_COMP_IDX set: " << compIdx; + _idxSO = compIdx[0]; + _idxSS_nad = compIdx[1]; + _idxSNH = compIdx[2]; + _idxSNO = compIdx[3]; + _idxSN2 = compIdx[4]; + _idxSALK = compIdx[5]; + _idxSI_nad = compIdx[6]; + _idxXI = compIdx[7]; + _idxXS = compIdx[8]; + _idxXH = compIdx[9]; + _idxXSTO = compIdx[10]; + _idxXA = compIdx[11]; + _idxXMI = compIdx[12]; + } + else + { + LOG(Debug) << "ASM3_COMP_IDX not set, using defaults"; + } + + // handle optional Aeration + const double V = paramProvider.getDouble("ASM3_V"); + const double IO2 = paramProvider.getDouble("ASM3_IO2"); + + if ( V > 0 ) + { + LOG(Debug) << "Activating reaction aeration in ASM3"; + _activeAeration = true; + _stoichiometry.resize(_nComp, 13); + } + else + { + LOG(Debug) << "Deactivating reaction aeration in ASM3"; + _activeAeration = false; + _stoichiometry.resize(_nComp, 12); + } + + _stoichiometry.setAll(0); + + // parameter set ASM3h + const double iNSI = paramProvider.getDouble("ASM3_INSI"); + const double iNSS = paramProvider.getDouble("ASM3_INSS"); + const double iNXI = paramProvider.getDouble("ASM3_INXI"); + const double iNXS = paramProvider.getDouble("ASM3_INXS"); + const double iNBM = paramProvider.getDouble("ASM3_INBM"); + + const double fSI = paramProvider.getDouble("ASM3_FSI"); + const double fXI = paramProvider.getDouble("ASM3_FXI"); + + const double YH_aer = paramProvider.getDouble("ASM3_YH_AER"); + if (YH_aer < 1e-16) + throw InvalidParameterException("ASM3 configuration: YH_aer must be bigger than zero"); + const double YH_anox = paramProvider.getDouble("ASM3_YH_ANOX"); + const double YSTO_aer = paramProvider.getDouble("ASM3_YSTO_AER"); + if (YSTO_aer < 1e-16) + throw InvalidParameterException("ASM3 configuration: YSTO_aer must be bigger than zero"); + const double YSTO_anox = paramProvider.getDouble("ASM3_YSTO_ANOX"); + const double YA = paramProvider.getDouble("ASM3_YA"); + if (YA < 1e-16) + throw InvalidParameterException("ASM3 configuration: YA must be bigger than zero"); + + const double fiSS_BM_prod = paramProvider.getDouble("ASM3_FISS_BM_PROD"); + const double iVSS_BM = paramProvider.getDouble("ASM3_IVSS_BM"); + const double iTSS_VSS_BM = paramProvider.getDouble("ASM3_ITSS_VSS_BM"); + + // internal variables + const double fXMI_BM = fiSS_BM_prod * fXI * iVSS_BM * (iTSS_VSS_BM - 1); + const double c1n = iNXS - iNSI * fSI - (1 - fSI) * iNSS; + const double c2n = iNSS; + const double c3n = iNSS; + const double c4n = -iNBM; + const double c5n = c4n; + const double c6n = -fXI * iNXI + iNBM; + const double c7n = c6n; + const double c10n = -1 / YA - iNBM; + const double c11n = -fXI * iNXI + iNBM; + const double c12n = c11n; + const double c3no = (YSTO_anox - 1.0) / (40.0 / 14.0); + const double c5no = (1.0 - 1 / YH_anox) / (40.0 / 14.0); + const double c7no = (fXI - 1) / (40.0 / 14.0); + const double c9no = -14.0 / 40.0; + const double c10no = 1 / YA; + const double c12no = c7no; + const double c1a = c1n / 14.0; + const double c2a = c2n / 14.0; + const double c3a = (c3n - c3no) / 14.0; + const double c4a = c4n / 14.0; + const double c5a = (c5n - c5no) / 14.0; + const double c6a = c6n / 14.0; + const double c7a = (c7n - c7no) / 14.0; + const double c9a = 1 / 40.0; + const double c10a = (c10n - c10no) / 14.0; + const double c11a = c11n / 14.0; + const double c12a = (c12n - c12no) / 14.0; + + // SO + _stoichiometry.native(_idxSO, 1) = YSTO_aer - 1; + _stoichiometry.native(_idxSO, 3) = 1 - 1 / YH_aer; + _stoichiometry.native(_idxSO, 5) = -1 * (1 - fXI); + _stoichiometry.native(_idxSO, 7) = -1; + _stoichiometry.native(_idxSO, 9) = -(64.0 / 14.0) * 1 / YA + 1; + _stoichiometry.native(_idxSO, 10) = -1 * (1 - fXI); + if (_activeAeration) + _stoichiometry.native(_idxSO, 12) = 1; + + + // SS_nad + _stoichiometry.native(_idxSS_nad, 0) = (1 - fSI); + _stoichiometry.native(_idxSS_nad, 1) = -1; + _stoichiometry.native(_idxSS_nad, 2) = -1; + + // SS_ad + if (_fractionate) + { + _stoichiometry.native(_idxSS_ad, 0) = (1 - fSI); + _stoichiometry.native(_idxSS_ad, 1) = -1; + _stoichiometry.native(_idxSS_ad, 2) = -1; + } + // SNH + _stoichiometry.native(_idxSNH, 0) = c1n; + _stoichiometry.native(_idxSNH, 1) = c2n; + _stoichiometry.native(_idxSNH, 2) = c3n; + _stoichiometry.native(_idxSNH, 3) = c4n; + _stoichiometry.native(_idxSNH, 4) = c5n; + _stoichiometry.native(_idxSNH, 5) = c6n; + _stoichiometry.native(_idxSNH, 6) = c7n; + _stoichiometry.native(_idxSNH, 9) = c10n; + _stoichiometry.native(_idxSNH, 10) = c11n; + _stoichiometry.native(_idxSNH, 11) = c12n; + + // SNO + _stoichiometry.native(_idxSNO, 2) = c3no; + _stoichiometry.native(_idxSNO, 4) = c5no; + _stoichiometry.native(_idxSNO, 6) = c7no; + _stoichiometry.native(_idxSNO, 8) = c9no; + _stoichiometry.native(_idxSNO, 9) = c10no; + _stoichiometry.native(_idxSNO, 11) = c12no; + + // SN2 + _stoichiometry.native(_idxSN2, 2) = -c3no; + _stoichiometry.native(_idxSN2, 4) = -c5no; + _stoichiometry.native(_idxSN2, 6) = -c7no; + _stoichiometry.native(_idxSN2, 8) = -c9no; + _stoichiometry.native(_idxSN2, 11) = -c12no; + + // SALK + _stoichiometry.native(_idxSALK, 0) = c1a; + _stoichiometry.native(_idxSALK, 1) = c2a; + _stoichiometry.native(_idxSALK, 2) = c3a; + _stoichiometry.native(_idxSALK, 3) = c4a; + _stoichiometry.native(_idxSALK, 4) = c5a; + _stoichiometry.native(_idxSALK, 5) = c6a; + _stoichiometry.native(_idxSALK, 6) = c7a; + _stoichiometry.native(_idxSALK, 8) = c9a; + _stoichiometry.native(_idxSALK, 9) = c10a; + _stoichiometry.native(_idxSALK, 10) = c11a; + _stoichiometry.native(_idxSALK, 11) = c12a; + + // SI_ad + _stoichiometry.native(_idxSI_nad, 0) = fSI; + + //SI_nad + if(_fractionate) + _stoichiometry.native(_idxSI_ad, 0) = fSI; + + // XI + _stoichiometry.native(_idxXI, 5) = fXI; + _stoichiometry.native(_idxXI, 6) = fXI; + _stoichiometry.native(_idxXI, 10) = fXI; + _stoichiometry.native(_idxXI, 11) = fXI; + + // XS + _stoichiometry.native(_idxXS, 0) = -1; + + // XH + _stoichiometry.native(_idxXH, 3) = 1; + _stoichiometry.native(_idxXH, 4) = 1; + _stoichiometry.native(_idxXH, 5) = -1; + _stoichiometry.native(_idxXH, 6) = -1; + + // XSTO + _stoichiometry.native(_idxXSTO, 1) = YSTO_aer; + _stoichiometry.native(_idxXSTO, 2) = YSTO_anox; + _stoichiometry.native(_idxXSTO, 3) = -1 / YH_aer; + _stoichiometry.native(_idxXSTO, 4) = -1 / YH_anox; + _stoichiometry.native(_idxXSTO, 7) = -1; + _stoichiometry.native(_idxXSTO, 8) = -1; + + // XA + _stoichiometry.native(_idxXA, 9) = 1; + _stoichiometry.native(_idxXA, 10) = -1; + _stoichiometry.native(_idxXA, 11) = -1; + + // XMI + _stoichiometry.native(_idxXMI, 5) = fXMI_BM; + _stoichiometry.native(_idxXMI, 6) = fXMI_BM; + _stoichiometry.native(_idxXMI, 10) = fXMI_BM; + _stoichiometry.native(_idxXMI, 11) = fXMI_BM; + + return true; + } + + template + int residualLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, + StateType const* y, ResidualType* res, const FactorType& factor, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + // Calculate fluxes + typedef typename DoubleActivePromoter::type flux_t; + BufferedArray fluxes = workSpace.array(_stoichiometry.columns()); + + const flux_t kh20 = static_cast::type>(p->kh20); + const flux_t T = static_cast::type>(p->T); + if (T < 0) + throw InvalidParameterException("ASM3 configuration: Temperature T must be non-negative"); + + flux_t io2; + flux_t V; + if (_activeAeration) + { + io2 = static_cast::type>(p->io2); + V = static_cast::type>(p->V); + if(V < 0) + throw InvalidParameterException("ASM3 configuration: Aeration volume V must be non-negative"); + } + const flux_t k_sto20 = static_cast::type>(p->k_sto20); + const flux_t kx = static_cast::type>(p->kx); + const flux_t kho2 = static_cast::type>(p->kho2); + const flux_t khss = static_cast::type>(p->khss); + const flux_t khn03 = static_cast::type>(p->khn03); + const flux_t etahno3 = static_cast::type>(p->etahno3); + const flux_t khnh4 = static_cast::type>(p->khnh4); + const flux_t khalk = static_cast::type>(p->khalk); + const flux_t khsto = static_cast::type>(p->khsto); + const flux_t muh2o = static_cast::type>(p->muh2o); + const flux_t etahend = static_cast::type>(p->etahend); + const flux_t bh20 = static_cast::type>(p->bh20); + const flux_t muAUT20 = static_cast::type>(p->muAUT20); + const flux_t kno2 = static_cast::type>(p->kno2); + const flux_t knnh4 = static_cast::type>(p->knnh4); + const flux_t knalk = static_cast::type>(p->knalk); + const flux_t baut20 = static_cast::type>(p->baut20); + const flux_t etanend = static_cast::type>(p->etanend); + + // derived parameters + const double ft04 = exp(-0.04 * (20.0 - static_cast(T))); + const double ft07 = exp(-0.06952 * (20.0 - static_cast(T))); + const double ft105 = exp(-0.105 * (20.0 - static_cast(T))); + const double k_sto = static_cast(k_sto20) * ft07; + const double muH = static_cast(muh2o) * ft07; + const double bH = static_cast(bh20) * ft07; + const double muAUT = static_cast(muAUT20) * ft105; + const double bAUT = static_cast(baut20) * ft105; + + StateType SO = y[_idxSO]; + StateType SS_nad = y[_idxSS_nad]; + StateType SS_ad = 0.0; + if (_fractionate) + SS_ad = y[_idxSS_ad]; + StateType SNH = y[_idxSNH]; + StateType SNO = y[_idxSNO]; + StateType SN2 = y[_idxSN2]; + StateType SALK = y[_idxSALK]; + // StateType SI = y[_idxSI]; // unused + // StateType XI = y[_idxXI]; // unused + StateType XS = y[_idxXS]; + StateType XH = y[_idxXH]; + StateType XSTO = y[_idxXSTO]; + StateType XA = y[_idxXA]; + // StateType XMI = y[_idxXMI]; // unused + + + // p1: Hydrolysis of organic structures + fluxes[0] = kh20 * ft04 * (XS/XH) / ((XS/XH) + kx) * XH; + if (XH < 0.1) + fluxes[0] = kh20 * ft04 * (XS/0.1) / ((XS/0.1) + kx) * XH; + + + // p2: Aerobic storage of SS + fluxes[1] = k_sto * SO / (SO + kho2) * ( SS_ad + SS_nad ) / ( ( SS_ad + SS_nad ) + khss ) * XH; + + // p3: Anoxic storage of SS + fluxes[2] = k_sto * etahno3 * kho2 / (SO + kho2) * ( SS_ad + SS_nad ) / ( ( SS_ad + SS_nad ) + khss ) * SNO / (SNO + khn03) * XH; + + // p4: Aerobic growth of heterotrophic biomass (XH) + fluxes[3] = muH * SO / (SO + kho2) * SNH / (SNH + khnh4) * SALK / (SALK + khalk) * ((XSTO/XH)) / (((XSTO/XH)) + khsto) * XH; + if (XH < 0.1) + fluxes[3] = muH * SO / (SO + kho2) * SNH / (SNH + khnh4) * SALK / (SALK + khalk) * (XSTO / 0.1) / ((XSTO / 0.1) + khsto) * XH; + + // p5: Anoxic growth of heterotrophic biomass (XH, denitrification) + fluxes[4] = muH * etahno3 * kho2 / (kho2 + SO) * SNH / (khnh4 + SNH) * SALK / (khalk + SALK) * (XSTO / XH) * 1 / (khsto + (XSTO / XH)) * SNO / (khn03 + SNO) * XH; + if (XH < 0.1) + fluxes[4] = muH * etahno3 * kho2 / (kho2 + SO) * SNH / (khnh4 + SNH) * SALK / (khalk + SALK) * (XSTO / 0.1) * 1 / (khsto + (XSTO/0.1)) * SNO / (khn03 + SNO) * XH; + + // r6: Aerobic endogenous respiration of heterotroph microorganisms (XH) + fluxes[5] = bH * SO / (SO + kho2) * XH; + + // r7: Anoxic endogenous respiration of heterotroph microorganisms (XH) + fluxes[6] = bH * etahend * kho2 / (SO + kho2) * SNO / (SNO + khn03) * XH; + + // r8: Aerobic respiration of internal cell storage products + fluxes[7] = bH * SO / (SO + kho2) * XSTO; + + // r9: Anoxic respiration of internal cell storage products + fluxes[8] = bH * etahend * kho2 / (SO + kho2) * SNO / (SNO + khn03) * XSTO; + + // r10: Aerobic growth of autotrophic biomass (XAUT, nitrification) + fluxes[9] = muAUT * SO / (SO + kno2) * SNH / (SNH + knnh4) * SALK / (SALK + knalk) * XA; + + // r11: Aerobic endogenous respiration of autotrophic biomass (XAUT) + fluxes[10] = bAUT * SO / (SO + kho2) * XA; + + // r12: Anoxic endogenous respiration of autotrophic biomass (XAUT) + fluxes[11] = bAUT * etanend * SNO / (SNO + khn03) * kho2 / (SO + kho2) * XA; + + // r13: Aeration + if (_activeAeration) + fluxes[12] = io2 / V; + + // Add reaction terms to residual + _stoichiometry.multiplyVector(static_cast(fluxes), factor, res); + + return 0; + } + + template + int residualCombinedImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, + StateType const* yLiquid, StateType const* ySolid, ResidualType* resLiquid, ResidualType* resSolid, double factor, LinearBufferAllocator workSpace) const + { + std::fill_n(resLiquid, _nComp, 0.0); + + if (_nTotalBoundStates == 0) + return 0; + + std::fill_n(resSolid, _nTotalBoundStates, 0.0); + + return 0; + } + + template + void jacobianLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, const RowIterator& jac, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + //parameters + const double kh20 = static_cast(p->kh20); + const double T = static_cast(p->T); + const double io2 = static_cast(p->io2); + const double V = static_cast(p->V); + const double k_sto20 = static_cast(p->k_sto20); + const double kx = static_cast(p->kx); + const double kho2 = static_cast(p->kho2); + const double khss = static_cast(p->khss); + const double khn03 = static_cast(p->khn03); + const double etahno3 = static_cast(p->etahno3); + const double khnh4 = static_cast(p->khnh4); + const double khalk = static_cast(p->khalk); + const double khsto = static_cast(p->khsto); + const double muh2o = static_cast(p->muh2o); + const double etahend = static_cast(p->etahend); + const double bh20 = static_cast(p->bh20); + const double muAUT20 = static_cast(p->muAUT20); + const double kno2 = static_cast(p->kno2); + const double knnh4 = static_cast(p->knnh4); + const double knalk = static_cast(p->knalk); + const double baut20 = static_cast(p->baut20); + const double etanend = static_cast(p->etanend); + + // derived parameters + const double ft04 = exp(-0.04 * (20.0 - static_cast(T))); + const double ft07 = exp(-0.06952 * (20.0 - static_cast(T))); + const double ft105 = exp(-0.105 * (20.0 - static_cast(T))); + const double k_sto = static_cast(k_sto20) * ft07; + const double muH = static_cast(muh2o) * ft07; + const double bH = static_cast(bh20) * ft07; + const double muAUT = static_cast(muAUT20) * ft105; + const double bAUT = static_cast(baut20) * ft105; + + double SO = y[_idxSO]; + double SS_nad = y[_idxSS_nad]; + double SS_ad = 0.0; + double SNH = y[_idxSNH]; + double SNO = y[_idxSNO]; + double SN2 = y[_idxSN2]; + double SALK = y[_idxSALK]; + double SI_nad = y[_idxSI_nad]; + double SI_ad = 0.0; + double XS = y[_idxXS]; + double XH = y[_idxXH]; + double XSTO = y[_idxXSTO]; + double XA = y[_idxXA]; + + if (_fractionate) + { + SS_ad = y[_idxSS_ad]; + SI_ad = y[_idxSI_ad]; + } + + // initialize jacobian + double d[13][15] = {}; + // p1: Hydrolysis: kh20 * ft04 * XS/XH_S / (XS/XH_S + kx) * XH; + d[0][_idxXS] = kh20 * ft04 + * XH / ((XS + XH * kx) + * (XS + XH * kx)) * XH; + d[0][_idxXH] = kh20 * ft04 + * (XS * XS) / ((XS + kx * XH) + * (XS + kx * XH)); + if (XH < 0.1) + { + d[0][_idxXS] = kh20 * ft04 + * 0.1 / ((XS + 0.1 * kx) * (XS + 0.1 * kx)) * XH; + d[0][_idxXH] = 0.0; + } + + // p2: Aerobic storage of SS: k_sto * SO / (SO + kho2) * ( SS_ad + SS_nad ) / ( ( SS_ad + SS_nad ) + khss ) * XH; + d[1][_idxSO] = k_sto + * (SS_ad + SS_nad) / ((SS_ad + SS_nad) + khss) + * kho2 / ((SO + kho2) * (SO + kho2)) * XH; + if (_fractionate) + d[1][_idxSS_ad] = k_sto + * SO / (SO + kho2) + * khss / ((SS_ad + SS_nad + khss) * (SS_ad + SS_nad + khss)) * XH; + d[1][_idxSS_nad] = k_sto + * SO / (SO + kho2) + * khss / ((SS_ad + SS_nad + khss) * (SS_ad + SS_nad + khss)) * XH; + d[1][_idxXH] = k_sto + * SO / (SO + kho2) + * (SS_ad + SS_nad) / ((SS_ad + SS_nad) + khss); + + // p3: Anoxic storage of SS: k_sto * etahno3 * kho2 / (SO + kho2) * ( SS_ad + SS_nad ) / ( ( SS_ad + SS_nad ) + khss ) * SNO / (SNO + khn03) * XH; + d[2][_idxSO] = k_sto * etahno3 + * -kho2 / ((SO + kho2) * (SO + kho2)) + * (SS_ad + SS_nad) / ((SS_ad + SS_nad) + khss) + * SNO / (SNO + khn03) * XH; + if(_fractionate) + d[2][_idxSS_ad] = k_sto * etahno3 + * kho2 / (SO + kho2) + * khss / ((SS_ad + SS_nad + khss) * (SS_ad + SS_nad + khss)) + * SNO / (SNO + khn03) * XH; + d[2][_idxSS_nad] = k_sto * etahno3 + * kho2 / (SO + kho2) + * khss / ((SS_ad + SS_nad + khss) * (SS_ad + SS_nad + khss)) + * SNO / (SNO + khn03) * XH; + d[2][_idxSNO] = k_sto * etahno3 + * kho2 / (SO + kho2) + * (SS_ad + SS_nad) / ((SS_ad + SS_nad) + khss) + * khn03 / ((SNO + khn03) * (SNO + khn03)) * XH; + d[2][_idxXH] = k_sto * etahno3 + * kho2 / (SO + kho2) + * (SS_ad + SS_nad) / ((SS_ad + SS_nad) + khss) + * SNO / (SNO + khn03); + + // p4: Aerobic growth: muH * SO / (SO + kho2) * SNH / (SNH + khnh4) * SALK / (SALK + khalk) * (XSTO/XH_S) / ((XSTO/XH_S) + khsto) * XH; + d[3][_idxSO] = muH + * kho2 / ((SO + kho2) * (SO + kho2)) + * SNH / (SNH + khnh4) + * SALK / (SALK + khalk) + * (XSTO / XH) / ((XSTO / XH) + khsto) * XH; + d[3][_idxSNH] = muH + * SO / (SO + kho2) + * khnh4 / ((SNH + khnh4) * (SNH + khnh4)) + * SALK / (SALK + khalk) + * (XSTO / XH) / ((XSTO / XH) + khsto) * XH; + d[3][_idxSALK] = muH + * SO / (SO + kho2) + * SNH / (SNH + khnh4) + * khalk / ((SALK + khalk) * (SALK + khalk)) + * (XSTO / XH) / ((XSTO / XH) + khsto) * XH; + d[3][_idxXSTO] = muH + * SO / (SO + kho2) + * SNH / (SNH + khnh4) + * SALK / (SALK + khalk) + * (khsto * XH) / ((XSTO + khsto * XH) * (XSTO + khsto * XH)) * XH; + d[3][_idxXH] = muH + * SO / (SO + kho2) + * SNH / (SNH + khnh4) + * SALK / (SALK + khalk) + * (XSTO * XSTO) / ((XSTO + khsto * XH) * (XSTO + khsto * XH)); + + if (XH < 0.1) + { + d[3][_idxSO] = muH + * -kho2 / ((SO + kho2) * (SO + kho2)) + * SNH / (SNH + khnh4) + * SALK / (SALK + khalk) + * (XSTO / 0.1) / ((XSTO / 0.1) + khsto) * 0.1; + d[3][_idxSNH] = muH + * SO / (SO + kho2) + * khnh4 / ((SNH + khnh4) * (SNH + khnh4)) + * khnh4 / ((SNH + khnh4) * (SNH + khnh4)) + * SALK / (SALK + khalk) + * (XSTO / 0.1) / ((XSTO / 0.1) + khsto) * 0.1; + d[3][_idxSALK] = muH + * SO / (SO + kho2) + * SNH / (SNH + khnh4) + * khalk / ((SALK + khalk) * (SALK + khalk)) + * (XSTO / 0.1) / ((XSTO / 0.1) + khsto) * 0.1; + d[3][_idxXSTO] = muH + * SO / (SO + kho2) + * SNH / (SNH + khnh4) + * SALK / (SALK + khalk) + * (khsto * 0.1) / ((XSTO + khsto * 0.1) * (XSTO + khsto * 0.1)) * 0.1; + d[3][_idxXH] = 0.0; + } + + // p5: Anoxic growth: muH * etahno3 * kho2 / (kho2 + SO) * SNH / (khnh4 + SNH) * SALK / (khalk + SALK) * SNO / (khn03 + SNO)* (XSTO / XH_S) / (khsto + (XSTO/XH)_S) * XH; + d[4][_idxSO] = muH * etahno3 + * -kho2 / ((kho2 + SO) * (kho2 + SO)) + * SNH / (khnh4 + SNH) + * SALK / (khalk + SALK) + * (XSTO / XH) / (khsto + (XSTO / XH)) + * SNO / (khn03 + SNO) * XH; + d[4][_idxSNH] = muH * etahno3 + * kho2 / (kho2 + SO) + * khnh4 / ((khnh4 + SNH) * (khnh4 + SNH)) + * SALK / (khalk + SALK) + * (XSTO / XH) / (khsto + (XSTO / XH)) + * SNO / (khn03 + SNO) * XH; + d[4][_idxSALK] = muH * etahno3 + * kho2 / (kho2 + SO) + * SNH / (khnh4 + SNH) + * khalk / ((khalk + SALK) * (khalk + SALK)) + * (XSTO / XH) / (khsto + (XSTO / XH)) + * SNO / (khn03 + SNO) * XH; + d[4][_idxSNO] = muH * etahno3 + * kho2 / (kho2 + SO) + * SNH / (khnh4 + SNH) + * SALK / (khalk + SALK) + * (XSTO / XH) / (khsto + (XSTO / XH)) + * khn03 / ((khn03 + SNO) * (khn03 + SNO)) * XH; + d[4][_idxXSTO] = muH * etahno3 + * kho2 / (kho2 + SO) + * SNH / (khnh4 + SNH) + * SALK / (khalk + SALK) + * (1.0 / XH) + * SNO / (khn03 + SNO) + * (khsto * XH * XH) / ((XSTO + khsto * XH) * (XSTO + khsto * XH)) * XH; + d[4][_idxXH] = muH * etahno3 + * kho2 / (kho2 + SO) + * SNH / (khnh4 + SNH) + * SALK / (khalk + SALK) + * SNO / (khn03 + SNO) + * (XSTO * XSTO) / ((XSTO + khsto * XH) * (XSTO + khsto * XH)); + + + if (XH < 0.1) + { + d[4][_idxSO] = muH * etahno3 + * -kho2 / ((kho2 + SO) * (kho2 + SO)) + * SNH / (khnh4 + SNH) + * SALK / (khalk + SALK) + * (XSTO / 0.1) / (khsto + (XSTO / 0.1)) + * SNO / (khn03 + SNO) * 0.1; + d[4][_idxSNH] = muH * etahno3 + * kho2 / (kho2 + SO) + * khnh4 / ((khnh4 + SNH) * (khnh4 + SNH)) + * SALK / (khalk + SALK) + * (XSTO / 0.1) / (khsto + (XSTO / 0.1)) + * SNO / (khn03 + SNO) * 0.1; + d[4][_idxSALK] = muH * etahno3 + * kho2 / (kho2 + SO) + * SNH / (khnh4 + SNH) + * khalk / ((khalk + SALK) * (khalk + SALK)) + * (XSTO / 0.1) / (khsto + (XSTO / 0.1)) + * SNO / (khn03 + SNO) * 0.1; + d[4][_idxSNO] = muH * etahno3 + * kho2 / (kho2 + SO) + * SNH / (khnh4 + SNH) + * SALK / (khalk + SALK) + * (XSTO / 0.1) / (khsto + (XSTO / 0.1)) + * khn03 / ((khn03 + SNO) * (khn03 + SNO)) * 0.1; + d[4][_idxXSTO] = muH * etahno3 + * kho2 / (kho2 + SO) + * SNH / (khnh4 + SNH) + * SALK / (khalk + SALK) + * (1.0 / 0.1) + * SNO / (khn03 + SNO) + * (khsto * 0.1 * 0.1) / ((XSTO + khsto * 0.1) * (XSTO + khsto * 0.1)) * 0.1; + d[4][_idxXH] = 0.0; + } + + //reaction6: bH * SO / (SO + kho2) * XH; + d[5][_idxSO] = bH + * kho2 / ((SO + kho2) * (SO + kho2)) * XH; + d[5][_idxXH] = bH * SO / (SO + kho2); + + //reaction7: bH * etahend * kho2 / (SO + kho2) * SNO / (SNO + khn03) * XH; + d[6][_idxSO] = bH * etahend + * -kho2 / ((SO + kho2) * (SO + kho2)) + * SNO / (SNO + khn03) * XH; + d[6][_idxSNO] = bH * etahend + * kho2 / (SO + kho2) + * khn03 / ((SNO + khn03) * (SNO + khn03)) * XH; + d[6][_idxXH] = bH * etahend + * kho2 / (SO + kho2) + * SNO / (SNO + khn03); + + //reaction8: bH * SO / (SO + kho2) * XSTO; + d[7][_idxSO] = bH + * kho2 / ((SO + kho2) * (SO + kho2)) * XSTO; + d[7][_idxXSTO] = bH * SO / (SO + kho2); + + //reaction9: bH * etahend * kho2 / (SO + kho2) * SNO / (SNO + khn03) * XSTO; + d[8][_idxSO] = bH * etahend + * -kho2 / ((SO + kho2) * (SO + kho2)) + * SNO / (SNO + khn03) * XSTO; + d[8][_idxSNO] = bH * etahend + * kho2 / (SO + kho2) + * khn03 / ((SNO + khn03) * (SNO + khn03)) * XSTO; + d[8][_idxXSTO] = bH * etahend + * kho2 / (SO + kho2) + * SNO / (SNO + khn03); + + //reaction10: muAUT * SO / (SO + kno2) * SNH / (SNH + knnh4) * SALK / (SALK + knalk) * XA; + d[9][_idxSO] = muAUT + * kno2 / ((SO + kno2) * (SO + kno2)) + * SNH / (SNH + knnh4) + * SALK / (SALK + knalk) * XA; + d[9][_idxSALK] = muAUT + * SO / (SO + kno2) + * SNH / (SNH + knnh4) + * knalk / ((SALK + knalk) * (SALK + knalk)) * XA; + d[9][_idxSNH] = muAUT + * SO / (SO + kno2) + * SALK / (SALK + knalk) + * knnh4 / ((SNH + knnh4) * (SNH + knnh4)) * XA; + d[9][_idxXA] = muAUT + * SO / (SO + kno2) + * SNH / (SNH + knnh4) + * SALK / (SALK + knalk); + + //reaction11: bAUT * SO / (SO + kho2) * XA; + d[10][_idxSO] = bAUT + * kho2 / ((SO + kho2) * (SO + kho2)) * XA; + d[10][_idxXA] = bAUT + * SO / (SO + kho2); + + //reaction12: bAUT * etanend * SNO / (SNO + khn03) * kho2 / (SO + kho2) * XA; + d[11][_idxSO] = bAUT * etanend + * SNO / (SNO + khn03) + * -kho2 / ((SO + kho2) * (SO + kho2)) * XA; + d[11][_idxSNO] = bAUT * etanend + * kho2 / (SO + kho2) + * khn03 / ((SNO + khn03) * (SNO + khn03)) * XA; + d[11][_idxXA] = bAUT * etanend + * SNO / (SNO + khn03) + * kho2 / (SO + kho2); + + //reaction13: Aeration: io2 / V; + // no jacobian terms + + + RowIterator curJac = jac; + for (size_t rIdx = 0; rIdx < _stoichiometry.columns(); rIdx++) + { + RowIterator curJac = jac; + for (int row = 0; row < _stoichiometry.rows(); ++row, ++curJac) + { + double colFactor = static_cast(_stoichiometry.native(row, rIdx)); + for (size_t compIdx = 0; compIdx < _stoichiometry.rows(); compIdx++) + { + if (_fractionate) + { + if (compIdx == _idxSI_ad) + colFactor *= SI_ad / (SI_ad + SI_nad); + if (compIdx == _idxSI_nad) + colFactor *= SI_nad / (SI_ad + SI_nad); + if (compIdx == _idxSS_ad) + colFactor *= SS_ad / (SS_ad + SS_nad); + if (compIdx == _idxSS_nad) + colFactor *= SS_nad / (SS_ad + SS_nad); + } + curJac[compIdx - static_cast(row)] += colFactor * d[rIdx][compIdx]; + } + } + } + + } + + template + void jacobianCombinedImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yLiquid, double const* ySolid, double factor, const RowIteratorLiquid& jacLiquid, const RowIteratorSolid& jacSolid, LinearBufferAllocator workSpace) const + { + } +}; + +typedef ActivatedSludgeModelThreeBase ActivatedSludgeModelThreeReaction; +typedef ActivatedSludgeModelThreeBase ExternalActivatedSludgeModelThreeReaction; + +namespace reaction +{ + void registerActivatedSludgeModelThreeReaction(std::unordered_map>& reactions) + { + reactions[ActivatedSludgeModelThreeReaction::identifier()] = []() { return new ActivatedSludgeModelThreeReaction(); }; + reactions[ExternalActivatedSludgeModelThreeReaction::identifier()] = []() { return new ExternalActivatedSludgeModelThreeReaction(); }; + } +} // namespace reaction + +} // namespace model + +} // namespace cadet diff --git a/test/MatrixHelper.hpp b/test/MatrixHelper.hpp index edb6c15b4..9dacdeea3 100644 --- a/test/MatrixHelper.hpp +++ b/test/MatrixHelper.hpp @@ -56,9 +56,9 @@ Matrix_t createBandMatrix(unsigned int rows, unsigned int lower, unsigned int up double val = 1.0; for (int row = 0; row < bm.rows(); ++row) { - const int lower = std::max(-static_cast(bm.lowerBandwidth()), -static_cast(row)); - const int upper = std::min(static_cast(bm.upperBandwidth()), static_cast(bm.rows() - row) - 1); - for (int col = lower; col <= upper; ++col) + const int low = std::max(-static_cast(bm.lowerBandwidth()), -static_cast(row)); + const int up = std::min(static_cast(bm.upperBandwidth()), static_cast(bm.rows() - row) - 1); + for (int col = low; col <= up; ++col) { bm.centered(row, col) = val; val += 1.0; diff --git a/test/ReactionModelTests.cpp b/test/ReactionModelTests.cpp index 5164bb941..dceb317e2 100644 --- a/test/ReactionModelTests.cpp +++ b/test/ReactionModelTests.cpp @@ -357,7 +357,7 @@ namespace reaction ad::prepareAdVectorSeedsForDenseMatrix(adY, 0, numDofs); ad::copyToAd(yState.data(), adY, numDofs); ad::resetAd(adRes, numDofs); - crm.model().residualCombinedAdd(1.0, 0u, ColumnPosition{0.0, 0.0, 0.0}, adY, adY + crm.nComp(), adRes, adRes + crm.nComp(), 1.0, crm.buffer()); + crm.model().residualCombinedAdd(1.0, 0u, ColumnPosition{ 0.0, 0.0, 0.0 }, adY, adY + crm.nComp(), adRes, adRes + crm.nComp(), 1.0, crm.buffer()); // Extract Jacobian cadet::linalg::DenseMatrix jacAD; @@ -367,30 +367,30 @@ namespace reaction // Calculate analytic Jacobian cadet::linalg::DenseMatrix jacAna; jacAna.resize(numDofs, numDofs); - crm.model().analyticJacobianCombinedAdd(1.0, 0u, ColumnPosition{0.0, 0.0, 0.0}, yState.data(), yState.data() + crm.nComp(), 1.0, jacAna.row(0), jacAna.row(crm.nComp()), crm.buffer()); + crm.model().analyticJacobianCombinedAdd(1.0, 0u, ColumnPosition{ 0.0, 0.0, 0.0 }, yState.data(), yState.data() + crm.nComp(), 1.0, jacAna.row(0), jacAna.row(crm.nComp()), crm.buffer()); cadet::test::checkJacobianPatternFD( [&](double const* lDir, double* res) -> void - { - std::fill_n(res, numDofs, 0.0); - crm.model().residualCombinedAdd(1.0, 0u, ColumnPosition{0.0, 0.0, 0.0}, lDir, lDir + crm.nComp(), res, res + crm.nComp(), 1.0, crm.buffer()); - }, - [&](double const* lDir, double* res) -> void - { - jacAna.multiplyVector(lDir, res); - }, + { + std::fill_n(res, numDofs, 0.0); + crm.model().residualCombinedAdd(1.0, 0u, ColumnPosition{ 0.0, 0.0, 0.0 }, lDir, lDir + crm.nComp(), res, res + crm.nComp(), 1.0, crm.buffer()); + }, + [&](double const* lDir, double* res) -> void + { + jacAna.multiplyVector(lDir, res); + }, yState.data(), dir.data(), colA.data(), colB.data(), numDofs, numDofs); cadet::test::checkJacobianPatternFD( [&](double const* lDir, double* res) -> void - { - std::fill_n(res, numDofs, 0.0); - crm.model().residualCombinedAdd(1.0, 0u, ColumnPosition{0.0, 0.0, 0.0}, lDir, lDir + crm.nComp(), res, res + crm.nComp(), 1.0, crm.buffer()); - }, - [&](double const* lDir, double* res) -> void - { - jacAD.multiplyVector(lDir, res); - }, + { + std::fill_n(res, numDofs, 0.0); + crm.model().residualCombinedAdd(1.0, 0u, ColumnPosition{ 0.0, 0.0, 0.0 }, lDir, lDir + crm.nComp(), res, res + crm.nComp(), 1.0, crm.buffer()); + }, + [&](double const* lDir, double* res) -> void + { + jacAD.multiplyVector(lDir, res); + }, yState.data(), dir.data(), colA.data(), colB.data(), numDofs, numDofs); // Check Jacobians against each other @@ -408,7 +408,7 @@ namespace reaction // Evaluate with AD ad::resetAd(adRes, numDofs); - crm.model().residualLiquidAdd(1.0, 0u, ColumnPosition{0.0, 0.0, 0.0}, adY, adRes, 1.0, crm.buffer()); + crm.model().residualLiquidAdd(1.0, 0u, ColumnPosition{ 0.0, 0.0, 0.0 }, adY, adRes, 1.0, crm.buffer()); // Extract Jacobian jacAD.setAll(0.0); @@ -416,33 +416,115 @@ namespace reaction // Calculate analytic Jacobian jacAna.setAll(0.0); - crm.model().analyticJacobianLiquidAdd(1.0, 0u, ColumnPosition{0.0, 0.0, 0.0}, yState.data(), 1.0, jacAna.row(0), crm.buffer()); + crm.model().analyticJacobianLiquidAdd(1.0, 0u, ColumnPosition{ 0.0, 0.0, 0.0 }, yState.data(), 1.0, jacAna.row(0), crm.buffer()); delete[] adY; delete[] adRes; cadet::test::checkJacobianPatternFD( [&](double const* lDir, double* res) -> void - { - std::fill_n(res, nComp, 0.0); - crm.model().residualLiquidAdd(1.0, 0u, ColumnPosition{0.0, 0.0, 0.0}, lDir, res, 1.0, crm.buffer()); - }, - [&](double const* lDir, double* res) -> void - { - jacAna.submatrixMultiplyVector(lDir, 0, 0, nComp, nComp, res); - }, + { + std::fill_n(res, nComp, 0.0); + crm.model().residualLiquidAdd(1.0, 0u, ColumnPosition{ 0.0, 0.0, 0.0 }, lDir, res, 1.0, crm.buffer()); + }, + [&](double const* lDir, double* res) -> void + { + jacAna.submatrixMultiplyVector(lDir, 0, 0, nComp, nComp, res); + }, yState.data(), dir.data(), colA.data(), colB.data(), nComp, nComp); cadet::test::checkJacobianPatternFD( [&](double const* lDir, double* res) -> void - { - std::fill_n(res, nComp, 0.0); - crm.model().residualLiquidAdd(1.0, 0u, ColumnPosition{0.0, 0.0, 0.0}, lDir, res, 1.0, crm.buffer()); - }, - [&](double const* lDir, double* res) -> void - { - jacAD.submatrixMultiplyVector(lDir, 0, 0, nComp, nComp, res); - }, + { + std::fill_n(res, nComp, 0.0); + crm.model().residualLiquidAdd(1.0, 0u, ColumnPosition{ 0.0, 0.0, 0.0 }, lDir, res, 1.0, crm.buffer()); + }, + [&](double const* lDir, double* res) -> void + { + jacAD.submatrixMultiplyVector(lDir, 0, 0, nComp, nComp, res); + }, + yState.data(), dir.data(), colA.data(), colB.data(), nComp, nComp); + + // Check Jacobians against each other + for (unsigned int row = 0; row < nComp; ++row) + { + for (unsigned int col = 0; col < nComp; ++col) + { + CAPTURE(row); + CAPTURE(col); + CHECK(jacAna.native(row, col) == makeApprox(jacAD.native(row, col), absTol, relTol)); + } + } + } + + void testLiquidReactionJacobianAD(const char* modelName, unsigned int nComp, unsigned int const* nBound, const char* config, double const* point, double absTol, double relTol) + { + ConfiguredDynamicReactionModel crm = ConfiguredDynamicReactionModel::create(modelName, nComp, nBound, config); + + const unsigned int numDofs = crm.nComp() + crm.numBoundStates(); + std::vector yState(numDofs, 0.0); + std::copy_n(point, numDofs, yState.data()); + + std::vector dir(numDofs, 0.0); + std::vector colA(numDofs, 0.0); + std::vector colB(numDofs, 0.0); + + // Enable AD + cadet::ad::setDirections(cadet::ad::getMaxDirections()); + cadet::active* adRes = new cadet::active[numDofs]; + cadet::active* adY = new cadet::active[numDofs]; + + // Evaluate with AD + ad::prepareAdVectorSeedsForDenseMatrix(adY, 0, numDofs); + ad::copyToAd(yState.data(), adY, numDofs); + ad::resetAd(adRes, numDofs); + + // Extract Jacobian + cadet::linalg::DenseMatrix jacAD; + jacAD.resize(numDofs, numDofs); + ad::extractDenseJacobianFromAd(adRes, 0, jacAD); + + // Calculate analytic Jacobian + cadet::linalg::DenseMatrix jacAna; + jacAna.resize(numDofs, numDofs); + + // Evaluate with AD + ad::resetAd(adRes, numDofs); + crm.model().residualLiquidAdd(1.0, 0u, ColumnPosition{ 0.0, 0.0, 0.0 }, adY, adRes, 1.0, crm.buffer()); + + // Extract Jacobian + jacAD.setAll(0.0); + ad::extractDenseJacobianFromAd(adRes, 0, jacAD); + + // Calculate analytic Jacobian + jacAna.setAll(0.0); + crm.model().analyticJacobianLiquidAdd(1.0, 0u, ColumnPosition{ 0.0, 0.0, 0.0 }, yState.data(), 1.0, jacAna.row(0), crm.buffer()); + + delete[] adY; + delete[] adRes; + + cadet::test::checkJacobianPatternFD( + [&](double const* lDir, double* res) -> void + { + std::fill_n(res, nComp, 0.0); + crm.model().residualLiquidAdd(1.0, 0u, ColumnPosition{ 0.0, 0.0, 0.0 }, lDir, res, 1.0, crm.buffer()); + }, + [&](double const* lDir, double* res) -> void + { + jacAna.submatrixMultiplyVector(lDir, 0, 0, nComp, nComp, res); + }, + yState.data(), dir.data(), colA.data(), colB.data(), nComp, nComp); + + cadet::test::checkJacobianPatternFD( + [&](double const* lDir, double* res) -> void + { + std::fill_n(res, nComp, 0.0); + crm.model().residualLiquidAdd(1.0, 0u, ColumnPosition{ 0.0, 0.0, 0.0 }, lDir, res, 1.0, crm.buffer()); + }, + [&](double const* lDir, double* res) -> void + { + jacAD.submatrixMultiplyVector(lDir, 0, 0, nComp, nComp, res); + }, yState.data(), dir.data(), colA.data(), colB.data(), nComp, nComp); // Check Jacobians against each other diff --git a/test/ReactionModelTests.hpp b/test/ReactionModelTests.hpp index aa0060144..0b4fed543 100644 --- a/test/ReactionModelTests.hpp +++ b/test/ReactionModelTests.hpp @@ -123,6 +123,18 @@ namespace reaction */ void testMichaelisMentenToSMAInhibitionMicroKinetic(const std::string configFilePathMM, const std::string configFilePathSMA, const double absTol, const double relTol); + /** + * @brief Checks the analytic Jacobians of the dynamic reaction model against AD + * @param [in] modelName Name of the reaction model + * @param [in] nComp Number of components + * @param [in] nBound Array with number of bound states for each component + * @param [in] config JSON string with reaction model parameters + * @param [in] point Liquid phase and solid phase values to check Jacobian at + * @param [in] absTol Absolute error tolerance + * @param [in] relTol Relative error tolerance + */ + void testDynamicJacobianAD(const char* modelName, unsigned int nComp, unsigned int const* nBound, const char* config, double const* point, double absTol = 0.0, double relTol = std::numeric_limits::epsilon() * 100.0); + /** * @brief Checks the analytic Jacobians of the dynamic reaction model against AD * @param [in] modelName Name of the reaction model @@ -133,7 +145,7 @@ namespace reaction * @param [in] absTol Absolute error tolerance * @param [in] relTol Relative error tolerance */ - void testDynamicJacobianAD(const char* modelName, unsigned int nComp, unsigned int const* nBound, const char* config, double const* point, double absTol = 0.0, double relTol = std::numeric_limits::epsilon() * 100.0); + void testLiquidReactionJacobianAD(const char* modelName, unsigned int nComp, unsigned int const* nBound, const char* config, double const* point, double absTol = 0.0, double relTol = std::numeric_limits::epsilon() * 100.0); /** * @brief Extends a model with dynamic reactions in each phase and particle type diff --git a/test/ReactionModels.cpp b/test/ReactionModels.cpp index 9c37134a8..7f477986d 100644 --- a/test/ReactionModels.cpp +++ b/test/ReactionModels.cpp @@ -180,6 +180,58 @@ TEST_CASE("MichaelisMenten kinetic analytic Jacobian vs AD with inhibition", "[M ); } +TEST_CASE("ASM3 analytic Jacobian vs AD", "[ASM3],[ReactionModel],[Jacobian],[AD]") +{ + const unsigned int nBound[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ,0, 0, 0 }; + unsigned int ncomp = 13; + const double point[] = { 1.0, 2.0, 1.4, 2.1, 0.2, 1.1, 1.8, 1.5, 1.0, 4.2, 1.4, 0.3, 1.4 }; + cadet::test::reaction::testLiquidReactionJacobianAD("ACTIVATED_SLUDGE_MODEL3", ncomp, nBound, + R"json({ + "ASM3_FISS_BM_PROD": 1.0, + "ASM3_FSI": 0.0, + "ASM3_YH_AER": 0.8, + "ASM3_YH_ANOX": 0.65, + "ASM3_YSTO_AER": 0.8375, + "ASM3_YSTO_ANOX": 0.7, + "ASM3_FXI": 0.2, + "ASM3_YA": 0.24, + "ASM3_KH20": 9.0, + "ASM3_KX": 1.0, + "ASM3_KSTO20": 12.0, + "ASM3_MU_H20": 3.0, + "ASM3_BH20": 0.33, + "ASM3_ETA_HNO3": 0.5, + "ASM3_KHO2": 0.2, + "ASM3_KHSS": 10.0, + "ASM3_KHNO3": 0.5, + "ASM3_KHNH4": 0.01, + "ASM3_KHALK": 0.1, + "ASM3_KHSTO": 0.1, + "ASM3_MU_AUT20": 1.12, + "ASM3_BAUT20": 0.18, + "ASM3_ETAH_END": 0.5, + "ASM3_ETAN_END": 0.5, + "ASM3_KNO2": 0.5, + "ASM3_KNNH4": 0.7, + "ASM3_KNALK": 0.5, + "ASM3_T": 12.0, + "ASM3_V": 1000.0, + "ASM3_IO2": 0.0, + "ASM3_INSI": 0.01, + "ASM3_INSS": 0.03, + "ASM3_INXI": 0.04, + "ASM3_INXS": 0.03, + "ASM3_INBM": 0.07, + "ASM3_IVSS_XI": 0.751879699, + "ASM3_IVSS_XS": 0.555555556, + "ASM3_IVSS_STO": 0.6, + "ASM3_IVSS_BM": 0.704225352, + "ASM3_ITSS_VSS_BM": 1.086956522 + })json", + point, 1e-15, 1e-15 + ); +} + TEST_CASE("MassActionLaw old interface vs. two separate reactions", "[MassActionLaw],[ReactionModel],[Simulation],[CI]") { std::string modelFilePath = std::string("/data/model_CSTR_reacMAL_3comp_nreac_2.json");