From a6a58dd513cfe271726b724e9d90302fca622d1b Mon Sep 17 00:00:00 2001 From: David Orme Date: Mon, 4 Sep 2023 12:33:51 +0100 Subject: [PATCH 1/8] Added Chen water density coefficients to PModelConst --- pyrealm/constants/pmodel_const.py | 42 +++++++++++++++++++++++++++++-- 1 file changed, 40 insertions(+), 2 deletions(-) diff --git a/pyrealm/constants/pmodel_const.py b/pyrealm/constants/pmodel_const.py index 51f96457..d0abe594 100644 --- a/pyrealm/constants/pmodel_const.py +++ b/pyrealm/constants/pmodel_const.py @@ -17,12 +17,18 @@ class PModelConst(ConstantsClass): value and units shown in brackets and the sources for default parameterisations are given below: - * **Density of water**. Values for the Tumlirz equation taken from Table 5 of - :cite:t:`Fisher:1975tm`: + * **Density of water using Fisher**. Values for the Tumlirz equation taken from + Table 5 of :cite:t:`Fisher:1975tm`: (:attr:`~pyrealm.constants.pmodel_const.PModelConst.fisher_dial_lambda`, :attr:`~pyrealm.constants.pmodel_const.PModelConst.fisher_dial_Po`, :attr:`~pyrealm.constants.pmodel_const.PModelConst.fisher_dial_Vinf`) + * **Density of water using Chen**. Values taken from :cite:t:`chen:2008a`: + (:attr:`~pyrealm.constants.pmodel_const.PModelConst.chen_po`, + :attr:`~pyrealm.constants.pmodel_const.PModelConst.chen_ko`, + :attr:`~pyrealm.constants.pmodel_const.PModelConst.chen_ca`, + :attr:`~pyrealm.constants.pmodel_const.PModelConst.chen_cb`) + * **Viscosity of water**. Values for the parameterisation taken from Table 2 and 3 of :cite:t:`Huber:2009fy`: (:attr:`~pyrealm.constants.pmodel_const.PModelConst.huber_tk_ast`, @@ -154,6 +160,38 @@ class PModelConst(ConstantsClass): r"""Temperature dependent Vinf parameterisation of the Tumlirz equation (:math:`V_{\infty}`).""" + # Chen water density + chen_po: NDArray[np.float32] = np.array( + [ + 0.99983952, + 6.788260e-5, + -9.08659e-6, + 1.022130e-7, + -1.35439e-9, + 1.471150e-11, + -1.11663e-13, + 5.044070e-16, + -1.00659e-18, + ] + ) + r"""Polynomial relationship of water density at 1 atm (kg/m^3) with temperature.""" + + chen_ko: NDArray[np.float32] = np.array( + [19652.17, 148.1830, -2.29995, 0.01281, -4.91564e-5, 1.035530e-7] + ) + r"""Polynomial relationship of bulk modulus of water at 1 atm (kg/m^3) with + temperature.""" + + chen_ca: NDArray[np.float32] = np.array( + [3.26138, 5.223e-4, 1.324e-4, -7.655e-7, 8.584e-10] + ) + r"""Polynomial temperature dependent coefficient :math:`c_{a}.""" + + chen_cb: NDArray[np.float32] = np.array( + [7.2061e-5, -5.8948e-6, 8.69900e-8, -1.0100e-9, 4.3220e-12] + ) + r"""Polynomial temperature dependent coefficient :math:`c_{b}.""" + # Huber simple_viscosity: bool = False """Boolean setting for use of simple viscosity calculations""" From 897771c783be509cba2861c0f10b10865d86793c Mon Sep 17 00:00:00 2001 From: David Orme Date: Mon, 4 Sep 2023 12:35:47 +0100 Subject: [PATCH 2/8] Added functions for chen water density --- pyrealm/pmodel/functions.py | 98 +++++++++++++++++++++++++++++++++++++ 1 file changed, 98 insertions(+) diff --git a/pyrealm/pmodel/functions.py b/pyrealm/pmodel/functions.py index ccaea951..43270b7a 100644 --- a/pyrealm/pmodel/functions.py +++ b/pyrealm/pmodel/functions.py @@ -11,6 +11,104 @@ from pyrealm.utilities import check_input_shapes +def calc_density_h2o_chen( + tc: NDArray, p: NDArray, const: PModelConst = PModelConst() +) -> NDArray: + """Calculate the density of water using Chen et al 2008. + + This function calculates the density of water at a given temperature and pressure + (kg/m^3) following :cite:t:`chen:2008a`. + + Args: + tc: Air temperature (°C) + p: Atmospheric pressure (Pa) + + Returns: + The calculated density of water + """ + + # Calculate density at 1 atm (kg/m^3): + po_coef = const.chen_po + po = po_coef[0] + po_coef[1] * tc + po += po_coef[2] * tc * tc + po += po_coef[3] * tc * tc * tc + po += po_coef[4] * tc * tc * tc * tc + po += po_coef[5] * tc * tc * tc * tc * tc + po += po_coef[6] * tc * tc * tc * tc * tc * tc + po += po_coef[7] * tc * tc * tc * tc * tc * tc * tc + po += po_coef[8] * tc * tc * tc * tc * tc * tc * tc * tc + + # Calculate bulk modulus at 1 atm (bar): + ko_coef = const.chen_ko + ko = ko_coef[0] + ko_coef[1] * tc + ko += ko_coef[2] * tc * tc + ko += ko_coef[3] * tc * tc * tc + ko += ko_coef[4] * tc * tc * tc * tc + ko += ko_coef[5] * tc * tc * tc * tc * tc + + # Calculate temperature dependent coefficients: + ca_coef = const.chen_ca + ca = ca_coef[0] + ca_coef[1] * tc + ca += ca_coef[2] * tc * tc + ca += ca_coef[3] * tc * tc * tc + ca += ca_coef[4] * tc * tc * tc * tc + + cb_coef = const.chen_cb + cb = cb_coef[0] + cb_coef[1] * tc + cb += cb_coef[2] * tc * tc + cb += cb_coef[3] * tc * tc * tc + cb += cb_coef[4] * tc * tc * tc * tc + + # Convert atmospheric pressure to bar (1 bar = 100000 Pa) + pbar = (1.0e-5) * p + + pw = ko + ca * pbar + cb * pbar**2.0 + pw /= ko + ca * pbar + cb * pbar**2.0 - pbar + pw *= (1e3) * po + return pw + + +def calc_density_h2o_chen_matrix( + tc: NDArray, + patm: NDArray, + const: PModelConst = PModelConst(), +) -> NDArray: + """Calculate the density of water using Chen et al 2008. + + This function calculates the density of water at a given temperature and pressure + (kg/m^3) following :cite:t:`chen:2008a`. + + Args: + tc: Air temperature (°C) + p: Atmospheric pressure (Pa) + + Returns: + The calculated density of water + """ + + # TODO - merge + tc_pow = np.power.outer(tc, np.arange(0, 9)) + + # Calculate density at 1 atm (kg/m^3): + rho_ref = np.sum(const.chen_po * tc_pow, axis=-1) + + # Calculate bulk modulus at 1 atm (bar): + bulk_mod_ref = np.sum(const.chen_ko * tc_pow[..., :6], axis=-1) + + # Calculate temperature dependent coefficients: + ca = np.sum(const.chen_ca * tc_pow[..., :5], axis=-1) + cb = np.sum(const.chen_cb * tc_pow[..., :5], axis=-1) + + # Convert atmospheric pressure to bar (1 bar = 100000 Pa) + pbar = (1.0e-5) * patm + + return ( + (bulk_mod_ref + ca * pbar + cb * pbar**2.0) + / (bulk_mod_ref + ca * pbar + cb * pbar**2.0 - pbar) + * ((1e3) * rho_ref) + ) + + def calc_density_h2o( tc: NDArray, patm: NDArray, From a717e143958115574746856ff23fcb9a5dcf4f4a Mon Sep 17 00:00:00 2001 From: David Orme Date: Mon, 4 Sep 2023 13:19:42 +0100 Subject: [PATCH 3/8] Added separate functions for existing fisher dial water density --- pyrealm/pmodel/functions.py | 158 ++++++++++++++++++++++++++++++++++++ 1 file changed, 158 insertions(+) diff --git a/pyrealm/pmodel/functions.py b/pyrealm/pmodel/functions.py index 43270b7a..ffa394aa 100644 --- a/pyrealm/pmodel/functions.py +++ b/pyrealm/pmodel/functions.py @@ -109,6 +109,164 @@ def calc_density_h2o_chen_matrix( ) +def calc_density_h2o_fisher( + tc: NDArray, + patm: NDArray, + const: PModelConst = PModelConst(), + safe: bool = True, +) -> NDArray: + """Calculate water density. + + Calculates the density of water as a function of temperature and atmospheric + pressure, using the Tumlirz Equation and coefficients calculated by + :cite:t:`Fisher:1975tm`. + + Args: + tc: air temperature, °C + patm: atmospheric pressure, Pa + const: Instance of :class:`~pyrealm.constants.pmodel_const.PModelConst`. + safe: Prevents the function from estimating density below -30°C, where the + function behaves poorly + + PModel Parameters: + lambda_: polynomial coefficients of Tumlirz equation (``fisher_dial_lambda``). + Po: polynomial coefficients of Tumlirz equation (``fisher_dial_Po``). + Vinf: polynomial coefficients of Tumlirz equation (``fisher_dial_Vinf``). + + Returns: + Water density as a float in (g cm^-3) + + Raises: + ValueError: if ``tc`` is less than -30°C and ``safe`` is True, or if the inputs + have incompatible shapes. + + Examples: + >>> round(calc_density_h2o(20, 101325), 3) + 998.206 + """ + + # It doesn't make sense to use this function for tc < 0, but in particular + # the calculation shows wild numeric instability between -44 and -46 that + # leads to numerous downstream issues - see the extreme values documentation. + if safe and np.nanmin(tc) < -30: + raise ValueError( + "Water density calculations below about -30°C are " + "unstable. See argument safe to calc_density_h2o" + ) + + # Check input shapes, shape not used + _ = check_input_shapes(tc, patm) + + # Calculate lambda, (bar cm^3)/g: + lambda_coef = const.fisher_dial_lambda + lambda_val = lambda_coef[0] + lambda_coef[1] * tc + lambda_val += lambda_coef[2] * tc * tc + lambda_val += lambda_coef[3] * tc * tc * tc + lambda_val += lambda_coef[4] * tc * tc * tc * tc + + # Calculate po, bar + po_coef = const.fisher_dial_Po + po_val = po_coef[0] + po_coef[1] * tc + po_val += po_coef[2] * tc * tc + po_val += po_coef[3] * tc * tc * tc + po_val += po_coef[4] * tc * tc * tc * tc + + # Calculate vinf, cm^3/g + vinf_coef = const.fisher_dial_Vinf + vinf_val = vinf_coef[0] + vinf_coef[1] * tc + vinf_val += vinf_coef[2] * tc * tc + vinf_val += vinf_coef[3] * tc * tc * tc + vinf_val += vinf_coef[4] * tc * tc * tc * tc + vinf_val += vinf_coef[5] * tc * tc * tc * tc * tc + vinf_val += vinf_coef[6] * tc * tc * tc * tc * tc * tc + vinf_val += vinf_coef[7] * tc * tc * tc * tc * tc * tc * tc + vinf_val += vinf_coef[8] * tc * tc * tc * tc * tc * tc * tc * tc + vinf_val += vinf_coef[9] * tc * tc * tc * tc * tc * tc * tc * tc * tc + + # Convert pressure to bars (1 bar <- 100000 Pa) + pbar = 1e-5 * patm + + # Calculate the specific volume (cm^3 g^-1): + spec_vol = vinf_val + lambda_val / (po_val + pbar) + + # Convert to density (g cm^-3) -> 1000 g/kg; 1000000 cm^3/m^3 -> kg/m^3: + rho = 1e3 / spec_vol + + return rho + + +def calc_density_h2o_fisher_matrix( + tc: NDArray, + patm: NDArray, + const: PModelConst = PModelConst(), + safe: bool = True, +) -> NDArray: + """Calculate water density. + + Calculates the density of water as a function of temperature and atmospheric + pressure, using the Tumlirz Equation and coefficients calculated by + :cite:t:`Fisher:1975tm`. + + Args: + tc: air temperature, °C + patm: atmospheric pressure, Pa + const: Instance of :class:`~pyrealm.constants.pmodel_const.PModelConst`. + safe: Prevents the function from estimating density below -30°C, where the + function behaves poorly + + PModel Parameters: + lambda_: polynomial coefficients of Tumlirz equation (``fisher_dial_lambda``). + Po: polynomial coefficients of Tumlirz equation (``fisher_dial_Po``). + Vinf: polynomial coefficients of Tumlirz equation (``fisher_dial_Vinf``). + + Returns: + Water density as a float in (g cm^-3) + + Raises: + ValueError: if ``tc`` is less than -30°C and ``safe`` is True, or if the inputs + have incompatible shapes. + + Examples: + >>> round(calc_density_h2o(20, 101325), 3) + 998.206 + """ + + # It doesn't make sense to use this function for tc < 0, but in particular + # the calculation shows wild numeric instability between -44 and -46 that + # leads to numerous downstream issues - see the extreme values documentation. + if safe and np.nanmin(tc) < np.array([-30]): + raise ValueError( + "Water density calculations below about -30°C are " + "unstable. See argument safe to calc_density_h2o" + ) + + # Check input shapes, shape not used + _ = check_input_shapes(tc, patm) + + # Get powers of tc, including tc^0 = 1 for constant terms + tc_pow = np.power.outer(tc, np.arange(0, 10)) + + # Calculate lambda, (bar cm^3)/g: + lambda_val = np.sum(const.fisher_dial_lambda * tc_pow[..., :5], axis=-1) + + # Calculate po, bar + po_val = np.sum(const.fisher_dial_Po * tc_pow[..., :5], axis=-1) + + # Calculate vinf, cm^3/g + vinf_val = np.sum(const.fisher_dial_Vinf * tc_pow, axis=-1) + + # Convert pressure to bars (1 bar <- 100000 Pa) + pbar = 1e-5 * patm + + # Calculate the specific volume (cm^3 g^-1): + spec_vol = vinf_val + lambda_val / (po_val + pbar) + + # Convert to density (g cm^-3) -> 1000 g/kg; 1000000 cm^3/m^3 -> kg/m^3: + rho = 1e3 / spec_vol + + return rho + + def calc_density_h2o( tc: NDArray, patm: NDArray, From 5de25e953950d13f1c0c3b4a55f42ba45f49f15a Mon Sep 17 00:00:00 2001 From: David Orme Date: Mon, 4 Sep 2023 17:28:21 +0100 Subject: [PATCH 4/8] Remove slow matrix implementations of calc_density_h2o methods and update docs for remaining methods --- pyrealm/pmodel/functions.py | 152 ++++++------------------------------ 1 file changed, 23 insertions(+), 129 deletions(-) diff --git a/pyrealm/pmodel/functions.py b/pyrealm/pmodel/functions.py index ffa394aa..90854129 100644 --- a/pyrealm/pmodel/functions.py +++ b/pyrealm/pmodel/functions.py @@ -19,12 +19,29 @@ def calc_density_h2o_chen( This function calculates the density of water at a given temperature and pressure (kg/m^3) following :cite:t:`chen:2008a`. + Warning: + The predictions from this function are numerically unstable around -58°C. + Args: tc: Air temperature (°C) p: Atmospheric pressure (Pa) + const: Instance of :class:`~pyrealm.constants.pmodel_const.PModelConst`. + + PModel Parameters: + chen_po: polynomial coefficients of water density equation at 1 atm. + chen_ko: polynomial coefficients of bulk modulus at 1 atm. + chen_ca: polynomial coefficients of temperature coefficient :math:`c_a`. + chen_cb: polynomial coefficients of temperature coefficient :math:`c_b`. Returns: - The calculated density of water + Water density as a float in (g cm^-3) + + Raises: + ValueError: if the inputs have incompatible shapes. + + Examples: + >>> round(calc_density_h2o_chen(20, 101325), 3) + 998.25 """ # Calculate density at 1 atm (kg/m^3): @@ -68,52 +85,10 @@ def calc_density_h2o_chen( return pw -def calc_density_h2o_chen_matrix( - tc: NDArray, - patm: NDArray, - const: PModelConst = PModelConst(), -) -> NDArray: - """Calculate the density of water using Chen et al 2008. - - This function calculates the density of water at a given temperature and pressure - (kg/m^3) following :cite:t:`chen:2008a`. - - Args: - tc: Air temperature (°C) - p: Atmospheric pressure (Pa) - - Returns: - The calculated density of water - """ - - # TODO - merge - tc_pow = np.power.outer(tc, np.arange(0, 9)) - - # Calculate density at 1 atm (kg/m^3): - rho_ref = np.sum(const.chen_po * tc_pow, axis=-1) - - # Calculate bulk modulus at 1 atm (bar): - bulk_mod_ref = np.sum(const.chen_ko * tc_pow[..., :6], axis=-1) - - # Calculate temperature dependent coefficients: - ca = np.sum(const.chen_ca * tc_pow[..., :5], axis=-1) - cb = np.sum(const.chen_cb * tc_pow[..., :5], axis=-1) - - # Convert atmospheric pressure to bar (1 bar = 100000 Pa) - pbar = (1.0e-5) * patm - - return ( - (bulk_mod_ref + ca * pbar + cb * pbar**2.0) - / (bulk_mod_ref + ca * pbar + cb * pbar**2.0 - pbar) - * ((1e3) * rho_ref) - ) - - def calc_density_h2o_fisher( tc: NDArray, patm: NDArray, const: PModelConst = PModelConst(), - safe: bool = True, ) -> NDArray: """Calculate water density. @@ -121,12 +96,13 @@ def calc_density_h2o_fisher( pressure, using the Tumlirz Equation and coefficients calculated by :cite:t:`Fisher:1975tm`. + Warning: + The predictions from this function are unstable around -45°C. + Args: tc: air temperature, °C patm: atmospheric pressure, Pa const: Instance of :class:`~pyrealm.constants.pmodel_const.PModelConst`. - safe: Prevents the function from estimating density below -30°C, where the - function behaves poorly PModel Parameters: lambda_: polynomial coefficients of Tumlirz equation (``fisher_dial_lambda``). @@ -137,23 +113,13 @@ def calc_density_h2o_fisher( Water density as a float in (g cm^-3) Raises: - ValueError: if ``tc`` is less than -30°C and ``safe`` is True, or if the inputs - have incompatible shapes. + ValueError: if the inputs have incompatible shapes. Examples: - >>> round(calc_density_h2o(20, 101325), 3) + >>> round(calc_density_h2o_fisher(20, 101325), 3) 998.206 """ - # It doesn't make sense to use this function for tc < 0, but in particular - # the calculation shows wild numeric instability between -44 and -46 that - # leads to numerous downstream issues - see the extreme values documentation. - if safe and np.nanmin(tc) < -30: - raise ValueError( - "Water density calculations below about -30°C are " - "unstable. See argument safe to calc_density_h2o" - ) - # Check input shapes, shape not used _ = check_input_shapes(tc, patm) @@ -195,78 +161,6 @@ def calc_density_h2o_fisher( return rho -def calc_density_h2o_fisher_matrix( - tc: NDArray, - patm: NDArray, - const: PModelConst = PModelConst(), - safe: bool = True, -) -> NDArray: - """Calculate water density. - - Calculates the density of water as a function of temperature and atmospheric - pressure, using the Tumlirz Equation and coefficients calculated by - :cite:t:`Fisher:1975tm`. - - Args: - tc: air temperature, °C - patm: atmospheric pressure, Pa - const: Instance of :class:`~pyrealm.constants.pmodel_const.PModelConst`. - safe: Prevents the function from estimating density below -30°C, where the - function behaves poorly - - PModel Parameters: - lambda_: polynomial coefficients of Tumlirz equation (``fisher_dial_lambda``). - Po: polynomial coefficients of Tumlirz equation (``fisher_dial_Po``). - Vinf: polynomial coefficients of Tumlirz equation (``fisher_dial_Vinf``). - - Returns: - Water density as a float in (g cm^-3) - - Raises: - ValueError: if ``tc`` is less than -30°C and ``safe`` is True, or if the inputs - have incompatible shapes. - - Examples: - >>> round(calc_density_h2o(20, 101325), 3) - 998.206 - """ - - # It doesn't make sense to use this function for tc < 0, but in particular - # the calculation shows wild numeric instability between -44 and -46 that - # leads to numerous downstream issues - see the extreme values documentation. - if safe and np.nanmin(tc) < np.array([-30]): - raise ValueError( - "Water density calculations below about -30°C are " - "unstable. See argument safe to calc_density_h2o" - ) - - # Check input shapes, shape not used - _ = check_input_shapes(tc, patm) - - # Get powers of tc, including tc^0 = 1 for constant terms - tc_pow = np.power.outer(tc, np.arange(0, 10)) - - # Calculate lambda, (bar cm^3)/g: - lambda_val = np.sum(const.fisher_dial_lambda * tc_pow[..., :5], axis=-1) - - # Calculate po, bar - po_val = np.sum(const.fisher_dial_Po * tc_pow[..., :5], axis=-1) - - # Calculate vinf, cm^3/g - vinf_val = np.sum(const.fisher_dial_Vinf * tc_pow, axis=-1) - - # Convert pressure to bars (1 bar <- 100000 Pa) - pbar = 1e-5 * patm - - # Calculate the specific volume (cm^3 g^-1): - spec_vol = vinf_val + lambda_val / (po_val + pbar) - - # Convert to density (g cm^-3) -> 1000 g/kg; 1000000 cm^3/m^3 -> kg/m^3: - rho = 1e3 / spec_vol - - return rho - - def calc_density_h2o( tc: NDArray, patm: NDArray, From 80c4c6248c1f707625d9c00236b6e41288080f4e Mon Sep 17 00:00:00 2001 From: David Orme Date: Mon, 4 Sep 2023 17:29:37 +0100 Subject: [PATCH 5/8] Update calc_density_h2o to wrap the chen and fisher methods --- pyrealm/pmodel/functions.py | 51 +++++++++++++------------------------ 1 file changed, 18 insertions(+), 33 deletions(-) diff --git a/pyrealm/pmodel/functions.py b/pyrealm/pmodel/functions.py index 90854129..685177ea 100644 --- a/pyrealm/pmodel/functions.py +++ b/pyrealm/pmodel/functions.py @@ -165,41 +165,41 @@ def calc_density_h2o( tc: NDArray, patm: NDArray, const: PModelConst = PModelConst(), + method: str = "fisher", safe: bool = True, ) -> NDArray: """Calculate water density. Calculates the density of water as a function of temperature and atmospheric - pressure, using the Tumlirz Equation and coefficients calculated by - :cite:t:`Fisher:1975tm`. + pressure. This function uses either the method provided by :cite:t:`Fisher:1975tm` + (:func:`~pyrealm.pmodel.functions.calc_density_h20_fisher`) or :cite:t:`chen:2008a` + (:func:`~pyrealm.pmodel.functions.calc_density_h20_chen`). + + The ``method`` argument can be used to switch between the ``fisher`` and ``chen`` + methods. This can also be set via the `water_density_method` argument to + :class:`~pyrealm.constants.PModelConsts`, which takes priority over the ``method`` + argument. Args: tc: air temperature, °C patm: atmospheric pressure, Pa const: Instance of :class:`~pyrealm.constants.pmodel_const.PModelConst`. safe: Prevents the function from estimating density below -30°C, where the - function behaves poorly - - PModel Parameters: - lambda_: polynomial coefficients of Tumlirz equation (``fisher_dial_lambda``). - Po: polynomial coefficients of Tumlirz equation (``fisher_dial_Po``). - Vinf: polynomial coefficients of Tumlirz equation (``fisher_dial_Vinf``). + functions are numerically unstable. Returns: Water density as a float in (g cm^-3) Raises: - ValueError: if ``tc`` is less than -30°C and ``safe`` is True, or if the inputs - have incompatible shapes. + ValueError: if ``tc`` contains values below -30°C and ``safe`` is True, or if + the inputs have incompatible shapes. Examples: >>> round(calc_density_h2o(20, 101325), 3) 998.206 """ - # It doesn't make sense to use this function for tc < 0, but in particular - # the calculation shows wild numeric instability between -44 and -46 that - # leads to numerous downstream issues - see the extreme values documentation. + # Safe guard against instability in functions at low temperature. if safe and np.nanmin(tc) < np.array([-30]): raise ValueError( "Water density calculations below about -30°C are " @@ -209,28 +209,13 @@ def calc_density_h2o( # Check input shapes, shape not used _ = check_input_shapes(tc, patm) - # Get powers of tc, including tc^0 = 1 for constant terms - tc_pow = np.power.outer(tc, np.arange(0, 10)) - - # Calculate lambda, (bar cm^3)/g: - lambda_val = np.sum(const.fisher_dial_lambda * tc_pow[..., :5], axis=-1) + if method == "fisher": + return calc_density_h2o_fisher(tc, patm, const) - # Calculate po, bar - po_val = np.sum(const.fisher_dial_Po * tc_pow[..., :5], axis=-1) + if method == "chen": + return calc_density_h2o_chen(tc, patm, const) - # Calculate vinf, cm^3/g - vinf_val = np.sum(const.fisher_dial_Vinf * tc_pow, axis=-1) - - # Convert pressure to bars (1 bar <- 100000 Pa) - pbar = 1e-5 * patm - - # Calculate the specific volume (cm^3 g^-1): - spec_vol = vinf_val + lambda_val / (po_val + pbar) - - # Convert to density (g cm^-3) -> 1000 g/kg; 1000000 cm^3/m^3 -> kg/m^3: - rho = 1e3 / spec_vol - - return rho + raise ValueError("Unknown method provided to calc_density_h2o") def calc_ftemp_arrh( From 4ff4354f2ddfd1aa00804b190a734002dc32fcb4 Mon Sep 17 00:00:00 2001 From: David Orme Date: Mon, 4 Sep 2023 17:52:28 +0100 Subject: [PATCH 6/8] Implemented calc_density_h20 as wrapper around methods, set via PModelConst --- pyrealm/constants/pmodel_const.py | 3 +++ pyrealm/pmodel/functions.py | 11 ++++------- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/pyrealm/constants/pmodel_const.py b/pyrealm/constants/pmodel_const.py index d0abe594..4643978e 100644 --- a/pyrealm/constants/pmodel_const.py +++ b/pyrealm/constants/pmodel_const.py @@ -130,6 +130,9 @@ class PModelConst(ConstantsClass): k_CtoK: float = 273.15 """Conversion from °C to K (:math:`CtoK` , 273.15, -)""" + water_density_method: str = "fisher" + """Set the method used for calculating water density ('fisher' or 'chen'.""" + # Fisher Dial fisher_dial_lambda: NDArray[np.float32] = np.array( [1788.316, 21.55053, -0.4695911, 0.003096363, -7.341182e-06] diff --git a/pyrealm/pmodel/functions.py b/pyrealm/pmodel/functions.py index 685177ea..d821e972 100644 --- a/pyrealm/pmodel/functions.py +++ b/pyrealm/pmodel/functions.py @@ -165,7 +165,6 @@ def calc_density_h2o( tc: NDArray, patm: NDArray, const: PModelConst = PModelConst(), - method: str = "fisher", safe: bool = True, ) -> NDArray: """Calculate water density. @@ -175,10 +174,8 @@ def calc_density_h2o( (:func:`~pyrealm.pmodel.functions.calc_density_h20_fisher`) or :cite:t:`chen:2008a` (:func:`~pyrealm.pmodel.functions.calc_density_h20_chen`). - The ``method`` argument can be used to switch between the ``fisher`` and ``chen`` - methods. This can also be set via the `water_density_method` argument to - :class:`~pyrealm.constants.PModelConsts`, which takes priority over the ``method`` - argument. + The `water_density_method` argument to :class:`~pyrealm.constants.PModelConsts` is + used to set which of the ``fisher`` or ``chen`` methods is used. Args: tc: air temperature, °C @@ -209,10 +206,10 @@ def calc_density_h2o( # Check input shapes, shape not used _ = check_input_shapes(tc, patm) - if method == "fisher": + if const.water_density_method == "fisher": return calc_density_h2o_fisher(tc, patm, const) - if method == "chen": + if const.water_density_method == "chen": return calc_density_h2o_chen(tc, patm, const) raise ValueError("Unknown method provided to calc_density_h2o") From 0283d75e070fc125c7413a41c988834c5de3e860 Mon Sep 17 00:00:00 2001 From: David Orme Date: Tue, 5 Sep 2023 15:35:49 +0100 Subject: [PATCH 7/8] Alternate version of calc_viscosity_h20; added some basic tests for pmodel.functions separate from rpmodel benchmarking --- pyrealm/pmodel/functions.py | 68 +++++++++++++++++++++++++++++ tests/pmodel/test_functions.py | 79 ++++++++++++++++++++++++++++++++++ 2 files changed, 147 insertions(+) create mode 100644 tests/pmodel/test_functions.py diff --git a/pyrealm/pmodel/functions.py b/pyrealm/pmodel/functions.py index d821e972..c5439cd8 100644 --- a/pyrealm/pmodel/functions.py +++ b/pyrealm/pmodel/functions.py @@ -845,6 +845,74 @@ def calc_viscosity_h2o( tbar = (tc + const.k_CtoK) / const.huber_tk_ast rbar = rho / const.huber_rho_ast + # Calculate mu0 (Eq. 11 & Table 2, Huber et al., 2009): + mu0 = const.huber_H_i[0] + const.huber_H_i[1] / tbar + mu0 += const.huber_H_i[2] / (tbar * tbar) + mu0 += const.huber_H_i[3] / (tbar * tbar * tbar) + mu0 = (1e2 * np.sqrt(tbar)) / mu0 + + # Calculate mu1 (Eq. 12 & Table 3, Huber et al., 2009): + ctbar = (1.0 / tbar) - 1.0 + mu1 = 0.0 + + # Iterate over the rows of the H_ij constants matrix + for row_idx in np.arange(const.huber_H_ij.shape[1]): + cf1 = ctbar**row_idx + cf2 = 0.0 + for col_idx in np.arange(const.huber_H_ij.shape[0]): + cf2 += const.huber_H_ij[col_idx, row_idx] * (rbar - 1.0) ** col_idx + mu1 += cf1 * cf2 + + mu1 = np.exp(rbar * mu1) + + # Calculate mu_bar (Eq. 2, Huber et al., 2009), assumes mu2 = 1 + mu_bar = mu0 * mu1 + + # Calculate mu (Eq. 1, Huber et al., 2009) + return mu_bar * const.huber_mu_ast # Pa s + + +def calc_viscosity_h2o_matrix( + tc: NDArray, + patm: NDArray, + const: PModelConst = PModelConst(), + simple: bool = False, +) -> NDArray: + r"""Calculate the viscosity of water. + + Calculates the viscosity of water (:math:`\eta`) as a function of temperature and + atmospheric pressure :cite:p:`Huber:2009fy`. + + Args: + tc: air temperature (°C) + patm: atmospheric pressure (Pa) + const: Instance of :class:`~pyrealm.constants.pmodel_const.PModelConst`. + simple: Use the simple formulation. + + Returns: + A float giving the viscosity of water (mu, Pa s) + + Examples: + >>> # Density of water at 20 degrees C and standard atmospheric pressure: + >>> round(calc_viscosity_h2o(20, 101325), 7) + 0.0010016 + """ + + # Check inputs, return shape not used + _ = check_input_shapes(tc, patm) + + if simple or const.simple_viscosity: + # The reference for this is unknown, but is used in some implementations + # so is included here to allow intercomparison. + return np.exp(-3.719 + 580 / ((tc + 273) - 138)) + + # Get the density of water, kg/m^3 + rho = calc_density_h2o(tc, patm, const=const) + + # Calculate dimensionless parameters: + tbar = (tc + const.k_CtoK) / const.huber_tk_ast + rbar = rho / const.huber_rho_ast + # Calculate mu0 (Eq. 11 & Table 2, Huber et al., 2009): tbar_pow = np.power.outer(tbar, np.arange(0, 4)) mu0 = (1e2 * np.sqrt(tbar)) / np.sum(np.array(const.huber_H_i) / tbar_pow, axis=-1) diff --git a/tests/pmodel/test_functions.py b/tests/pmodel/test_functions.py new file mode 100644 index 00000000..2e2f23cd --- /dev/null +++ b/tests/pmodel/test_functions.py @@ -0,0 +1,79 @@ +"""Test the pmodel functions + +TODO - note that there are parallel tests in test_pmodel that benchmark against the +rpmodel outputs and test a wider range of inputs. Those could be moved here. These tests +check the size of outputs and that the results meet a simple benchmark value. +""" +import numpy as np +import pytest + + +@pytest.mark.parametrize(argnames="shape", argvalues=[(1,), (6, 9), (4, 7, 3)]) +def test_calc_density_h20_fisher(shape): + """Test the fisher method""" + from pyrealm.pmodel.functions import calc_density_h2o_fisher + + rho = calc_density_h2o_fisher( + np.full(shape, fill_value=20), np.full(shape, fill_value=101325) + ) + + assert np.allclose(rho.round(3), np.full(shape, fill_value=998.206)) + + +@pytest.mark.parametrize(argnames="shape", argvalues=[(1,), (6, 9), (4, 7, 3)]) +def test_calc_density_h20_chen(shape): + """Test the chen method""" + from pyrealm.pmodel.functions import calc_density_h2o_chen + + rho = calc_density_h2o_chen( + np.full(shape, fill_value=20), np.full(shape, fill_value=101325) + ) + + assert np.allclose(rho.round(3), np.full(shape, fill_value=998.25)) + + +@pytest.mark.parametrize( + argnames="const_args, exp", + argvalues=[ + pytest.param(None, 998.206, id="defaults"), + pytest.param("fisher", 998.206, id="fisher_via_const"), + pytest.param("chen", 998.25, id="chen_via_const"), + ], +) +def test_calc_density_h20(const_args, exp): + """Test the wrapper method dispatches as expected""" + from pyrealm.constants import PModelConst + from pyrealm.pmodel.functions import calc_density_h2o + + args = {} + + if const_args is not None: + args["const"] = PModelConst(water_density_method=const_args) + + rho = calc_density_h2o(np.array([20]), np.array([101325]), **args) + + assert rho.round(3) == exp + + +@pytest.mark.parametrize(argnames="shape", argvalues=[(1,), (6, 9), (4, 7, 3)]) +def test_calc_viscosity_h20(shape): + """Test the viscosity calculation""" + from pyrealm.pmodel.functions import calc_viscosity_h2o + + eta = calc_viscosity_h2o( + np.full(shape, fill_value=20), np.full(shape, fill_value=101325) + ) + + assert np.allclose(eta.round(7), np.full(shape, fill_value=0.0010016)) + + +@pytest.mark.parametrize(argnames="shape", argvalues=[(1,), (6, 9), (4, 7, 3)]) +def test_calc_viscosity_h20_matrix(shape): + """Test the viscosity calculation""" + from pyrealm.pmodel.functions import calc_viscosity_h2o_matrix + + eta = calc_viscosity_h2o_matrix( + np.full(shape, fill_value=20), np.full(shape, fill_value=101325) + ) + + assert np.allclose(eta.round(7), np.full(shape, fill_value=0.0010016)) From 968761b36a3874bcec0aa3dc2e54a722f33292f5 Mon Sep 17 00:00:00 2001 From: David Orme Date: Tue, 5 Sep 2023 15:53:33 +0100 Subject: [PATCH 8/8] Docstring and sphinx fixes --- pyrealm/constants/pmodel_const.py | 4 ++-- pyrealm/pmodel/functions.py | 9 +++++---- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/pyrealm/constants/pmodel_const.py b/pyrealm/constants/pmodel_const.py index 4643978e..d5d16765 100644 --- a/pyrealm/constants/pmodel_const.py +++ b/pyrealm/constants/pmodel_const.py @@ -188,12 +188,12 @@ class PModelConst(ConstantsClass): chen_ca: NDArray[np.float32] = np.array( [3.26138, 5.223e-4, 1.324e-4, -7.655e-7, 8.584e-10] ) - r"""Polynomial temperature dependent coefficient :math:`c_{a}.""" + r"""Polynomial temperature dependent coefficient :math:`c_{a}`.""" chen_cb: NDArray[np.float32] = np.array( [7.2061e-5, -5.8948e-6, 8.69900e-8, -1.0100e-9, 4.3220e-12] ) - r"""Polynomial temperature dependent coefficient :math:`c_{b}.""" + r"""Polynomial temperature dependent coefficient :math:`c_{b}`.""" # Huber simple_viscosity: bool = False diff --git a/pyrealm/pmodel/functions.py b/pyrealm/pmodel/functions.py index c5439cd8..b5793841 100644 --- a/pyrealm/pmodel/functions.py +++ b/pyrealm/pmodel/functions.py @@ -171,11 +171,12 @@ def calc_density_h2o( Calculates the density of water as a function of temperature and atmospheric pressure. This function uses either the method provided by :cite:t:`Fisher:1975tm` - (:func:`~pyrealm.pmodel.functions.calc_density_h20_fisher`) or :cite:t:`chen:2008a` - (:func:`~pyrealm.pmodel.functions.calc_density_h20_chen`). + (:func:`~pyrealm.pmodel.functions.calc_density_h2o_fisher`) or :cite:t:`chen:2008a` + (:func:`~pyrealm.pmodel.functions.calc_density_h2o_chen`). - The `water_density_method` argument to :class:`~pyrealm.constants.PModelConsts` is - used to set which of the ``fisher`` or ``chen`` methods is used. + The `water_density_method` argument to + :class:`~pyrealm.constants.pmodel_const.PModelConst` is used to set which of the + ``fisher`` or ``chen`` methods is used. Args: tc: air temperature, °C