From a18999c1804286066c27b68f6ea557b4f4c21e14 Mon Sep 17 00:00:00 2001 From: vgro Date: Fri, 27 Oct 2023 10:44:23 +0100 Subject: [PATCH 01/18] first draft wind added --- docs/source/api/abiotic/wind.md | 24 + docs/source/index.md | 1 + tests/models/abiotic/test_wind.py | 318 ++++++++++ virtual_rainforest/models/abiotic/__init__.py | 5 + .../models/abiotic/constants.py | 22 + virtual_rainforest/models/abiotic/wind.py | 590 ++++++++++++++++++ 6 files changed, 960 insertions(+) create mode 100644 docs/source/api/abiotic/wind.md create mode 100644 tests/models/abiotic/test_wind.py create mode 100644 virtual_rainforest/models/abiotic/wind.py diff --git a/docs/source/api/abiotic/wind.md b/docs/source/api/abiotic/wind.md new file mode 100644 index 000000000..aae9c0aad --- /dev/null +++ b/docs/source/api/abiotic/wind.md @@ -0,0 +1,24 @@ +--- +jupytext: + cell_metadata_filter: -all + formats: md:myst + main_language: python + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.13.8 +kernelspec: + display_name: vr_python3 + language: python + name: vr_python3 +--- + +#  API for the {mod}`~virtual_rainforest.models.abiotic.wind` module + +```{eval-rst} +.. automodule:: virtual_rainforest.models.abiotic.wind + :autosummary: + :members: + :special-members: __init__ +``` diff --git a/docs/source/index.md b/docs/source/index.md index c9d7f9578..65ad75613 100644 --- a/docs/source/index.md +++ b/docs/source/index.md @@ -108,6 +108,7 @@ team. Abiotic Mechanistic Model Abiotic Mechanistic Constants Abiotic Mechanistic Tools + Abiotic Mechanistic Wind Hydrology Overview Hydrology Model Hydrology Above-ground diff --git a/tests/models/abiotic/test_wind.py b/tests/models/abiotic/test_wind.py new file mode 100644 index 000000000..2939652a5 --- /dev/null +++ b/tests/models/abiotic/test_wind.py @@ -0,0 +1,318 @@ +"""Test module for wind.py.""" + + +import numpy as np +import pytest +from xarray import DataArray + +from virtual_rainforest.models.abiotic.constants import AbioticConsts + + +@pytest.fixture +def dummy_data(): + """Creates a dummy data object for use in wind tests. + + One grid cell has no vegetation, two grid cells represent a range of values. + """ + + from virtual_rainforest.core.data import Data + from virtual_rainforest.core.grid import Grid + + # Setup the data object with two cells. + grid = Grid(cell_nx=3, cell_ny=1) + data = Data(grid) + + # Add the required data. + data["canopy_height"] = DataArray([0, 10, 50], dims=["cell_id"]) + data["wind_layer_heights"] = DataArray( + [[0.1, 5, 15], [1, 15, 25], [5, 25, 55]], dims=["cell_id", "layers"] + ) + data["leaf_area_index"] = DataArray( + [ + [0, 1, 5], + [0, 2, 5], + [0, 3, 5], + ], + dims=["cell_id", "layers"], + ) + data["air_temperature"] = DataArray([30, 20, 30], dims=["cell_id"]) + data["atmospheric_pressure"] = DataArray([101, 102, 103], dims=["cell_id"]) + data["sensible_heat_flux_topofcanopy"] = DataArray([100, 50, 10], dims=["cell_id"]) + data["friction_velocity"] = DataArray([12, 5, 2], dims=["cell_id"]) + data["wind_speed_ref"] = DataArray([0, 5, 10], dims=["cell_id"]) + + return data + + +def test_calculate_zero_plane_displacement(dummy_data): + """Test if calculated correctly and set to zero without vegetation.""" + + from virtual_rainforest.models.abiotic.wind import calculate_zero_plane_displacement + + plant_area_index = dummy_data.data["leaf_area_index"].sum(dim="layers") + result = calculate_zero_plane_displacement( + canopy_height=dummy_data.data["canopy_height"].to_numpy(), + plant_area_index=plant_area_index.to_numpy(), + zero_plane_scaling_parameter=AbioticConsts.zero_plane_scaling_parameter, + ) + + np.testing.assert_allclose(result, np.array([0.0, 8.620853, 43.547819])) + + +def test_calculate_roughness_length_momentum(dummy_data): + """Test roughness length governing momentum transfer.""" + + from virtual_rainforest.models.abiotic.wind import ( + calculate_roughness_length_momentum, + ) + + plant_area_index = dummy_data.data["leaf_area_index"].sum(dim="layers") + result = calculate_roughness_length_momentum( + canopy_height=dummy_data.data["canopy_height"].to_numpy(), + plant_area_index=plant_area_index.to_numpy(), + zero_plane_displacement=np.array([0, 8, 43]), + substrate_surface_drag_coefficient=( + AbioticConsts.substrate_surface_drag_coefficient + ), + roughness_element_drag_coefficient=( + AbioticConsts.roughness_element_drag_coefficient + ), + roughness_sublayer_depth_parameter=( + AbioticConsts.roughness_sublayer_depth_parameter + ), + max_ratio_wind_to_friction_velocity=( + AbioticConsts.max_ratio_wind_to_friction_velocity + ), + von_karman_constant=AbioticConsts.von_karmans_constant, + ) + + np.testing.assert_allclose( + result, np.array([0.003, 0.434662, 1.521318]), rtol=1e-3, atol=1e-3 + ) + + +def test_calculate_diabatic_correction_above(dummy_data): + """Test diabatic correction factors for heat and momentum.""" + + from virtual_rainforest.models.abiotic.wind import ( + calculate_diabatic_correction_above, + ) + + result = calculate_diabatic_correction_above( + temperature=dummy_data["air_temperature"].to_numpy(), + atmospheric_pressure=dummy_data["atmospheric_pressure"].to_numpy(), + sensible_heat_flux=dummy_data["sensible_heat_flux_topofcanopy"].to_numpy(), + friction_velocity=dummy_data["friction_velocity"].to_numpy(), + wind_heights_above_canopy=np.array([1, 15, 50]), + zero_plane_displacement=np.array([0, 8, 43]), + celsius_to_kelvin=AbioticConsts.celsius_to_kelvin, + standard_pressure=AbioticConsts.standard_pressure, + standard_mole=AbioticConsts.standard_mole, + molar_heat_capacity_air=AbioticConsts.molar_heat_capacity_air, + specific_heat_equ_factor_1=AbioticConsts.specific_heat_equ_factor_1, + specific_heat_equ_factor_2=AbioticConsts.specific_heat_equ_factor_2, + von_karmans_constant=AbioticConsts.von_karmans_constant, + ) + + exp_result_h = np.array([7.5154e-05, 6.2562e-04, 3.0957e-04]) + exp_result_m = np.array([4.5093e-05, 3.7534e-04, 1.857e-04]) + np.testing.assert_allclose(result["psi_h"], exp_result_h, rtol=1e-4, atol=1e-4) + np.testing.assert_allclose(result["psi_m"], exp_result_m, rtol=1e-4, atol=1e-4) + + +def test_calculate_wind_attenuation_coefficient(dummy_data): + """Test wind attenuation coefficient with and without vegetation.""" + + from virtual_rainforest.models.abiotic.wind import ( + calculate_wind_attenuation_coefficient, + ) + + plant_area_index = dummy_data.data["leaf_area_index"].sum(dim="layers") + result = calculate_wind_attenuation_coefficient( + canopy_height=dummy_data.data["canopy_height"].to_numpy(), + plant_area_index=plant_area_index, + mixing_length=np.array([0, 0.4, 1.5]), + drag_coefficient=AbioticConsts.drag_coefficient, + relative_turbulence_intensity=AbioticConsts.relative_turbulence_intensity, + diabatic_correction_factor_below=AbioticConsts.diabatic_correction_factor_below, + ) + + np.testing.assert_allclose(result, np.array([0.0, 5.91608, 7.302967])) + + +def test_calculate_mixing_length(dummy_data): + """Test mixing length with and without vegetation.""" + + from virtual_rainforest.models.abiotic.wind import calculate_mixing_length + + result = calculate_mixing_length( + canopy_height=dummy_data.data["canopy_height"].to_numpy(), + zero_plane_displacement=np.array([0, 8, 43]), + roughness_length_momentum=np.array([0.003, 0.435, 1.521]), + mixing_length_factor=AbioticConsts.mixing_length_factor, + ) + + np.testing.assert_allclose( + result, np.array([0, 0.419, 1.467]), rtol=1e-3, atol=1e-3 + ) + + +def test_wind_log_profile(): + """Test log wind profile.""" + + from virtual_rainforest.models.abiotic.wind import wind_log_profile + + result = wind_log_profile( + height=np.array([[1, 15, 50], [5, 20, 60]]), + zeroplane_displacement=np.array([0, 8, 43]), + roughness_length_momentum=np.array([0.003, 0.435, 1.521]), + diabatic_correction_momentum=np.array([-0.1, 0.0, 0.1]), + ) + np.testing.assert_allclose( + result, + np.array([[5.709, 2.778, 1.627], [7.319, 3.317, 2.514]]), + rtol=1e-3, + atol=1e-3, + ) + + +def test_calculate_wind_profile(dummy_data): + """Test if wind profile is calculated correctly for different height.""" + + from virtual_rainforest.models.abiotic.wind import calculate_wind_profile + + plant_area_index = dummy_data.data["leaf_area_index"].sum(dim="layers") + result = calculate_wind_profile( + wind_speed_ref=dummy_data["wind_speed_ref"].to_numpy(), + wind_layer_heights=np.array([[5, 12, 60], [0.1, 8, 10]]), + reference_height=10.0, + attenuation_coefficient=np.array([0, 4, 7]), + plant_area_index=plant_area_index, + canopy_height=dummy_data["canopy_height"].to_numpy(), + diabatic_correction_momentum=np.array([0, 0, 0]), + ground_layer_vegetation_height=0.1, + roughness_sublayer_depth_parameter=( + AbioticConsts.roughness_sublayer_depth_parameter + ), + zero_plane_scaling_parameter=AbioticConsts.zero_plane_scaling_parameter, + substrate_surface_drag_coefficient=( + AbioticConsts.substrate_surface_drag_coefficient + ), + roughness_element_drag_coefficient=( + AbioticConsts.roughness_element_drag_coefficient + ), + max_ratio_wind_to_friction_velocity=( + AbioticConsts.max_ratio_wind_to_friction_velocity + ), + von_karmans_constant=AbioticConsts.von_karmans_constant, + ) + + np.testing.assert_allclose( + result, + np.array([[0, 7.935657, 24.623732], [0, 1.471923, 0.036979]]), + rtol=1e-3, + atol=1e-3, + ) + + +# def test_calculate_wind_above_canopy(dummy_data): +# """Test wind above canopy.""" + +# from virtual_rainforest.models.abiotic import wind + +# result = wind.calculate_wind_above_canopy( +# wind_heights=DataArray( +# [[10, 20, 30], [20, 30, 40], [40, 50, 60]], dims=["heights", "cell_id"] +# ), +# zero_plane_displacement=DataArray([0, 10, 10], dims="cell_id"), +# roughness_length_momentum=DataArray([0.003, 0.4, 1.5], dims="cell_id"), +# diabatic_correction_momentum_above=DataArray( +# [[1, 1, 1], [1, 1, 1], [1, 1, 1]], dims=["heights", "cell_id"] +# ), +# friction_velocity=dummy_data["friction_velocity"], +# ) + +# xr.testing.assert_allclose( +# result, +# DataArray( +# [ +# [244.3518425, 265.14625792, 285.94067333], +# [41.23594781, 49.90028757, 58.56462732], +# [13.95133583, 15.97866137, 18.53278949], +# ], +# dims=["cell_id", "heights"], +# ), +# ) + + +# def test_calculate_wind_below_canopy(dummy_data): +# """Test wind within canopy.""" + +# from virtual_rainforest.models.abiotic import wind + +# result = wind.calculate_wind_below_canopy( +# canopy_node_heights=dummy_data.data["canopy_node_heights"], +# wind_profile_above=DataArray( +# [ +# [244.0, 265.0, 285.0], +# [41.0, 49.0, 58.0], +# [13.0, 15.0, 18.0], +# ], +# dims=["cell_id", "heights"], +# ), +# wind_attenuation_coefficient=DataArray([0, 0.1, 0.3], dims=["cell_id"]), +# canopy_height=dummy_data.data["canopy_height"], +# ) + +# xr.testing.assert_allclose( +# result, +# DataArray( +# [ +# [0.0, 0.0, 0.0], +# [52.48057, 60.973724, 67.386386], +# [13.334728, 15.492744, 16.450761], +# ], +# dims=["cell_id", "canopy_layers"], +# ), +# ) + + +# def test_calculate_wind_profile(dummy_data): +# """Test wind profile above and within canopy.""" + +# from virtual_rainforest.models.abiotic import wind + +# result = wind.calculate_wind_profile( +# wind_heights=DataArray( +# [[10, 20, 30], [20, 30, 40], [40, 50, 60]], dims=["heights", "cell_id"] +# ), +# canopy_node_heights=dummy_data.data["canopy_node_heights"], +# data=dummy_data, +# ) + +# # check wind above canopy +# xr.testing.assert_allclose( +# result[0], +# DataArray( +# DataArray( +# [ +# [244.351842, 265.146258, 285.940673], +# [46.458135, 54.341054, 62.595567], +# [0.0, 0.0, 13.311866], +# ], +# dims=["cell_id", "heights"], +# ), +# ), +# ) +# # check wind below canopy +# xr.testing.assert_allclose( +# result[1], +# DataArray( +# [ +# [0.0, 0.0, 0.0], +# [5.950523e-02, 2.030195e03, 2.135631e06], +# [6.086929e-03, 2.846548e-01, 1.325216e00], +# ], +# dims=["cell_id", "canopy_layers"], +# ), +# ) diff --git a/virtual_rainforest/models/abiotic/__init__.py b/virtual_rainforest/models/abiotic/__init__.py index f5e12b00d..b8795fc16 100644 --- a/virtual_rainforest/models/abiotic/__init__.py +++ b/virtual_rainforest/models/abiotic/__init__.py @@ -15,6 +15,11 @@ * The :mod:`~virtual_rainforest.models.abiotic.abiotic_tools` submodule contains a set of general functions that are shared across submodules in the :mod:`~virtual_rainforest.models.abiotic` model. + +* The :mod:`~virtual_rainforest.models.abiotic.wind` submodule calculates the + above- and within-canopy wind profiles for the Virtual Rainforest. These profiles will + determine the exchange of heat, water, and :math:`CO_{2}` between soil and atmosphere + below the canopy as well as the exchange with the atmsophere above the canopy. """ # noqa: D205, D415 from virtual_rainforest.models.abiotic.abiotic_model import AbioticModel # noqa: F401 diff --git a/virtual_rainforest/models/abiotic/constants.py b/virtual_rainforest/models/abiotic/constants.py index b915a79af..a4bab440d 100644 --- a/virtual_rainforest/models/abiotic/constants.py +++ b/virtual_rainforest/models/abiotic/constants.py @@ -38,3 +38,25 @@ class AbioticConsts(ConstantsDataclass): latent_heat_vap_equ_factor_2: float = 33.91 """Factor in calculation of latent heat of vaporisation, (Henderson-Sellers, 1984). """ + + zero_plane_scaling_parameter: float = 7.5 + """Control parameter for scaling zero displacement/height + :cite:p:`raupach_simplified_1994`.""" + substrate_surface_drag_coefficient: float = 0.003 + """Substrate-surface drag coefficient :cite:p:`maclean_microclimc_2021`.""" + roughness_element_drag_coefficient: float = 0.3 + """Roughness-element drag coefficient :cite:p:`maclean_microclimc_2021`.""" + roughness_sublayer_depth_parameter: float = 0.193 + """Parameter characterizes the roughness sublayer depth + :cite:p:`maclean_microclimc_2021`.""" + max_ratio_wind_to_friction_velocity: float = 0.3 + """Maximum ratio of wind velocity to friction velocity + :cite:p:`maclean_microclimc_2021`.""" + drag_coefficient: float = 0.2 + """Drag coefficient :cite:p:`maclean_microclimc_2021`.""" + relative_turbulence_intensity: float = 0.5 + """Relative turbulence intensity :cite:p:`maclean_microclimc_2021`.""" + diabatic_correction_factor_below: float = 1 + "Diabatic correction factor below canopy." + mixing_length_factor: float = 0.32 + """Factor in calculation of mixing length :cite:p:`maclean_microclimc_2021`.""" diff --git a/virtual_rainforest/models/abiotic/wind.py b/virtual_rainforest/models/abiotic/wind.py new file mode 100644 index 000000000..3836069a8 --- /dev/null +++ b/virtual_rainforest/models/abiotic/wind.py @@ -0,0 +1,590 @@ +r"""The wind module calculates the above- and within-canopy wind profiles for the +Virtual Rainforest. These profiles will determine the exchange of heat, water, and +:math:`CO_{2}` between soil and atmosphere below the canopy as well as the exchange with +the atmsophere above the canopy. +TODO: add sanity checks, errors and logging +TODO check equations for above profile, and howdoes it need adjustment for below +TODO add function that alls all the steps and returns info for data object +TODO select wind layer heights +""" # noqa: D205, D415 + +from typing import Union + +import numpy as np +from numpy.typing import NDArray + +from virtual_rainforest.models.abiotic.abiotic_tools import ( + calculate_molar_density_air, + calculate_specific_heat_air, +) + + +def calculate_zero_plane_displacement( + canopy_height: NDArray[np.float32], + plant_area_index: NDArray[np.float32], + zero_plane_scaling_parameter: float, +) -> NDArray[np.float32]: + """Calculate zero plane displacement. + + The zero plane displacement height of a vegetated surface is the height at which the + wind speed would go to zero if the logarithmic wind profile was maintained from the + outer flow all the way down to the surface (that is, in the absence of the + vegetation when the value is set to zero). Implementation after + :cite:t:`maclean_microclimc_2021`. + + Args: + canopy_height: canopy height, [m] + plant_area_index: leaf area index, [m m-1] + zero_plane_scaling_parameter: Control parameter for scaling d/h + :cite:p:`raupach_simplified_1994` + + Returns: + zero place displacement height, [m] + """ + + # Calculate in presence of vegetation + plant_area_index_displacement = np.where( + plant_area_index > 0, plant_area_index, np.nan + ) + + scale_displacement = np.sqrt( + zero_plane_scaling_parameter * plant_area_index_displacement + ) + + zero_place_displacement = ( + (1 - (1 - np.exp(-scale_displacement)) / scale_displacement) * canopy_height, + ) + + # no displacement in absence of vegetation + return np.nan_to_num(zero_place_displacement, nan=0.0).squeeze() + + +def calculate_roughness_length_momentum( + canopy_height: NDArray[np.float32], + plant_area_index: NDArray[np.float32], + zero_plane_displacement: NDArray[np.float32], + substrate_surface_drag_coefficient: float, + roughness_element_drag_coefficient: float, + roughness_sublayer_depth_parameter: float, + max_ratio_wind_to_friction_velocity: float, + von_karman_constant: float, +) -> NDArray[np.float32]: + """Calculate roughness length governing momentum transfer. + + Roughness length is defined as the height at which the mean velocity is zero due to + substrate roughness. Real surfaces such as the ground or vegetation are not smooth + and often have varying degrees of roughness. Roughness length accounts for that + effect. Implementation after :cite:t:`maclean_microclimc_2021`. + + Args: + canopy_height: canopy height, [m] + plant_area_index: plant area index, [m m-1] + substrate_surface_drag_coefficient: substrate-surface drag coefficient + zero_plane_displacement: height above ground within the canopy where the wind + profile extrapolates to zero, [m] + roughness_element_drag_coefficient: roughness-element drag coefficient + roughness_sublayer_depth_parameter: parameter characterizes the roughness + sublayer depth + max_ratio_wind_to_friction_velocity: Maximum ratio of wind velocity to friction + velocity + + Returns: + momentum roughness length, [m] + """ + + # calculate ratio of wind velocity to friction velocity + ratio_wind_to_friction_velocity = np.sqrt( + substrate_surface_drag_coefficient + + (roughness_element_drag_coefficient * plant_area_index) / 2 + ) + + # if the ratio of wind velocity to friction velocity is larger than the set maximum, + # set the value to set maximum + set_maximum_ratio = np.where( + ratio_wind_to_friction_velocity > max_ratio_wind_to_friction_velocity, + max_ratio_wind_to_friction_velocity, + ratio_wind_to_friction_velocity, + ) + + # calculate initial roughness length + initial_roughness_length = (canopy_height - zero_plane_displacement) * np.exp( + -von_karman_constant * (1 / set_maximum_ratio) + - roughness_sublayer_depth_parameter + ) + + # if roughness smaller than the substrate surface drag coefficient, set to value to + # the substrate surface drag coefficient + roughness_length = np.where( + initial_roughness_length < substrate_surface_drag_coefficient, + substrate_surface_drag_coefficient, + initial_roughness_length, + ) + + return np.where(roughness_length <= 0, 0.01, roughness_length) + + +def calculate_diabatic_correction_above( + temperature: NDArray[np.float32], + atmospheric_pressure: NDArray[np.float32], + sensible_heat_flux: NDArray[np.float32], + friction_velocity: NDArray[np.float32], + wind_heights_above_canopy: NDArray[np.float32], + zero_plane_displacement: NDArray[np.float32], + celsius_to_kelvin: float, + standard_mole: float, + standard_pressure: float, + molar_heat_capacity_air: float, + specific_heat_equ_factor_1: float, + specific_heat_equ_factor_2: float, + von_karmans_constant: float, +) -> dict[str, NDArray[np.float32]]: + r"""Calculates the diabatic correction factors for momentum and heat above canopy. + + Diabatic correction factor for heat and momentum are used to adjust wind profiles + for surface heating and cooling :cite:p:`maclean_microclimc_2021`. When the surface + is strongly heated, the diabatic correction factor for momemtum :math:`psi_{m}` + becomes negative and drops to values of around -1.5. In contrast, when the surface + is much cooler than the air above it, it increases to values around 4. + + Args: + temperature: 2m temperature + atmospheric_pressure: atmospheric pressure, [kPa] + sensible_heat_flux: sensible heat flux from canopy to atmosphere above, + [W m-2], # TODO: could be the top entry of the general sensible heat flux + friction_velocity: friction velocity + wind_heights: vector of heights for which wind speed is calculated, [m] + zero_plane_displacement: height above ground within the canopy where the wind + profile extrapolates to zero, [m] + celsius_to_kelvin: Factor to convert temperature in Celsius to absolute + temperature in Kelvin + standard_mole: Moles of ideal gas in 1 m^3 air at standard atmosphere + standard_pressure: Standard atmospheric pressure, [kPa] + molar_heat_capacity_air: molar heat capacity of air, [J mol-1 C-1] + specific_heat_equ_factor_1: Factor in calculation of molar specific heat of air + specific_heat_equ_factor_2: Factor in calculation of molar specific heat of air + von_karmans_constant + + Returns: + diabatic correction factors for heat (psi_h) and momentum (psi_m) transfer + """ + + molar_density_air = calculate_molar_density_air( + temperature=temperature, + atmospheric_pressure=atmospheric_pressure, + standard_mole=standard_mole, + standard_pressure=standard_pressure, + celsius_to_kelvin=celsius_to_kelvin, + ) + specific_heat_air = calculate_specific_heat_air( + temperature=temperature, + molar_heat_capacity_air=molar_heat_capacity_air, + specific_heat_equ_factor_1=specific_heat_equ_factor_1, + specific_heat_equ_factor_2=specific_heat_equ_factor_2, + ) + + # calculate atmospheric stability + stability = ( + von_karmans_constant + * (wind_heights_above_canopy - zero_plane_displacement) + * sensible_heat_flux + ) / ( + molar_density_air + * specific_heat_air + * (temperature + celsius_to_kelvin) + * friction_velocity + ) + + stable_condition = 6 * np.log(1 - stability) + unstable_condition = -2 * np.log((1 + np.sqrt(1 - 16 * stability)) / 2) + + diabatic_correction_heat = np.where( + stability < 0, stable_condition, unstable_condition + ) + + diabatic_correction_momentum = np.where( + stability < 0, diabatic_correction_heat, 0.6 * diabatic_correction_heat + ) + + # Apply upper threshold + diabatic_correction_momentum = np.minimum(diabatic_correction_momentum, 5) + diabatic_correction_heat = np.minimum(diabatic_correction_heat, 5) + + return {"psi_m": diabatic_correction_momentum, "psi_h": diabatic_correction_heat} + + +def calculate_mixing_length( + canopy_height: NDArray[np.float32], + zero_plane_displacement: NDArray[np.float32], + roughness_length_momentum: NDArray[np.float32], + mixing_length_factor: float, +) -> NDArray[np.float32]: + """Calculate mixing length for canopy air transport. + + The mixing length is used to calculate turbulent air transport inside vegetated + canopies. It is made equivalent to the above canopy value at the canopy surface. In + absence of vegetation, it is set to zero. + Implementation after :cite:t:`maclean_microclimc_2021`. + + Args: + canopy_height: canopy height, [m] + leaf_area_index: leaf area index, [m m-1] + zero_plane_displacement: height above ground within the canopy where the wind + profile extrapolates to zero, [m] + roughness_length_momentum: momentum roughness length, [m] + mixing_length_factor: Factor in calculation of mixing length + + Returns: + mixing length for canopy air transport, [m] + """ + + mixing_length = ( + mixing_length_factor * (canopy_height - zero_plane_displacement) + ) / np.log((canopy_height - zero_plane_displacement) / roughness_length_momentum) + + return np.nan_to_num(mixing_length, nan=0) + + +def calculate_wind_attenuation_coefficient( + canopy_height: NDArray[np.float32], + plant_area_index: NDArray[np.float32], + mixing_length: NDArray[np.float32], + drag_coefficient: float, + relative_turbulence_intensity: float, + diabatic_correction_factor_below: float, +) -> NDArray[np.float32]: + """Calculate wind attenuation coefficient. + + The wind attenuation coefficient describes how wind is slowed down by the presence + of vegetation. In absence of vegetation, the coefficient is set to zero. + Implementation after :cite:t:`maclean_microclimc_2021`. + + Args: + canopy_height: canopy height, [m] + plant_area_index: plant area index, [m m-1] + mixing_length: mixing length for canopy air transport, [m] + drag_coefficient: drag coefficient + relative_turbulence_intensity: relative turbulence intensity + diabatic_correction_factor_below: diabatic correction factor below canopy + + Returns: + wind attenuation coefficient + """ + + intermediate_coefficient1 = ( + (drag_coefficient * plant_area_index * canopy_height) + / (2 * mixing_length * relative_turbulence_intensity), + ) + + intermediate_coefficient2 = np.power(intermediate_coefficient1, 0.5) + + atten_coeff = intermediate_coefficient2 * np.power( + diabatic_correction_factor_below, 0.5 + ) + + return np.nan_to_num(atten_coeff, nan=0).squeeze() + + +def wind_log_profile( + height: Union[float, NDArray[np.float32]], + zeroplane_displacement: Union[float, NDArray[np.float32]], + roughness_length_momentum: Union[float, NDArray[np.float32]], + diabatic_correction_momentum: Union[float, NDArray[np.float32]], +) -> NDArray[np.float32]: + """Calculate logarithmic wind profile. + + Args: + height: array of heights for which wind speed is calculated, [m] + zeroplane_displacement: height above ground within the canopy where the wind + profile extrapolates to zero, [m] + roughness_length_momentum: momentum roughness length, [m] + diabatic_correction_momentum: diabaric correction factor for momentum + + Returns: + logarithmic wind profile, [m s-1] + """ + + wind_profile = ( + np.log((height - zeroplane_displacement) / roughness_length_momentum) + + diabatic_correction_momentum, + ) + + return np.nan_to_num(wind_profile, nan=1).squeeze() + + +def calculate_wind_profile( + wind_speed_ref: NDArray[np.float32], + reference_height: float, + wind_layer_heights: NDArray[np.float32], + attenuation_coefficient: NDArray[np.float32], + plant_area_index: NDArray[np.float32], + canopy_height: NDArray[np.float32], + diabatic_correction_momentum: NDArray[np.float32], + ground_layer_vegetation_height: Union[float, NDArray[np.float32]], + roughness_sublayer_depth_parameter: float, + zero_plane_scaling_parameter: float, + substrate_surface_drag_coefficient: float, + roughness_element_drag_coefficient: float, + max_ratio_wind_to_friction_velocity: float, + von_karmans_constant: float, +) -> NDArray[np.float32]: + """Calculate wind speed at any given height from wind speed at reference height. + + Wind profile above the canopy dictates heat and vapor exchange between the canopy + and air above it, and therefore ultimately determine temperature and vapor profiles. + The wind profile above canopy typically follows a logarithmic height profile, which + extrapolates to zero roughly two thirds of the way to the top of the canopy. The + profile itself is thus dependent on the height of the canopy, but also on the + roughness of the vegetation layer, which causes wind shear. We follow the + implementation by :cite:t:`campbell_introduction_1998` as described in + :cite:t:`maclean_microclimc_2021`. + + Args: + wind_speed_ref: Wind speed at reference height, [m s-1] + reference_height: Height of wind measurement, [m] + wind_layer_heights: Heights for which wind speed is required, [m] + attenuation_coefficient: Attenuation coefficient as returned by + :func:`~virtual_rainforest.models.abiotic.wind.calculate_wind_attenuation_coefficient` + plant_area_index: Plant area index, [m m-1] + canopy_height: Canopy height, [m] + diabatic_correction_momentum: Diabatic correction factor for momentum as + returned by + :func:`~virtual_rainforest.models.abiotic.wind.calculate_diabatic_correction_above` + ground_layer_vegetation_height: Height of ground vegetation layer below canopy, + [m] + roughness_sublayer_depth_parameter: parameter characterizes the roughness + sublayer depth, [m] + zero_plane_scaling_parameter: Control parameter for scaling zero + displacement/height + substrate_surface_drag_coefficient: Substrate-surface drag coefficient + roughness_element_drag_coefficient: Roughness-element drag coefficient + max_ratio_wind_to_friction_velocity: Maximum ratio of wind velocity to friction + velocity + von_karmans_constant: Von Karman's constant, unitless, describes the logarithmic + velocity profile of a turbulent fluid near a no-slip boundary + + Returns: + wind speed at required heights above and below canopy, [m s-1] + """ + + # Calculate zero-plane displacement + zeroplane_displacement = calculate_zero_plane_displacement( + canopy_height=canopy_height, + plant_area_index=plant_area_index, + zero_plane_scaling_parameter=zero_plane_scaling_parameter, + ) + + # Calculate roughness length for momentum, [m] + roughness_length = calculate_roughness_length_momentum( + canopy_height=canopy_height, + plant_area_index=plant_area_index, + zero_plane_displacement=zeroplane_displacement, + substrate_surface_drag_coefficient=substrate_surface_drag_coefficient, + roughness_element_drag_coefficient=roughness_element_drag_coefficient, + roughness_sublayer_depth_parameter=roughness_sublayer_depth_parameter, + max_ratio_wind_to_friction_velocity=max_ratio_wind_to_friction_velocity, + von_karman_constant=von_karmans_constant, + ) + + roughness_length_ground_vegetation = 0.1 * ground_layer_vegetation_height + + # reference wind profile + wind_profile_reference = wind_log_profile( + height=reference_height, + zeroplane_displacement=zeroplane_displacement, + roughness_length_momentum=roughness_length, + diabatic_correction_momentum=diabatic_correction_momentum, + ) + + friction_velocity = von_karmans_constant * (wind_speed_ref / wind_profile_reference) + + wind_profile_above = wind_log_profile( + height=wind_layer_heights, + zeroplane_displacement=zeroplane_displacement, + roughness_length_momentum=roughness_length, + diabatic_correction_momentum=diabatic_correction_momentum, + ) + wind_profile_above[wind_profile_above < 0.55] = 0.55 + + wind_profile_below = wind_log_profile( + height=wind_layer_heights, + zeroplane_displacement=zeroplane_displacement, + roughness_length_momentum=roughness_length, + diabatic_correction_momentum=diabatic_correction_momentum, + ) + + wind_profile = np.where( + wind_layer_heights > canopy_height, wind_profile_above, wind_profile_below + ) + wind_speed = (friction_velocity / von_karmans_constant) * wind_profile + + # Required height above 10% of canopy hgt, but below top of canopy + wind_10percent_above = wind_speed * np.exp( + attenuation_coefficient * ((wind_layer_heights / canopy_height) - 1) + ) + + wind_speed = np.where( + (wind_layer_heights > (0.1 * canopy_height)) + & (wind_layer_heights < canopy_height), + wind_10percent_above, + wind_speed, + ) + + # Required height below 10% of canopy hgt + wind_10percent_below = wind_speed * np.exp( + attenuation_coefficient * (((0.1 * canopy_height) / canopy_height) - 1) + ) + + wind_profile_reference_below = wind_log_profile( + height=0.1 * canopy_height, + zeroplane_displacement=0.0, + roughness_length_momentum=roughness_length_ground_vegetation, + diabatic_correction_momentum=diabatic_correction_momentum, + ) + + friction_velocity_below = von_karmans_constant * ( + wind_10percent_below / wind_profile_reference_below + ) + + wind_profile_10percent_below = wind_log_profile( + height=wind_layer_heights, + zeroplane_displacement=zeroplane_displacement, + roughness_length_momentum=roughness_length_ground_vegetation, + diabatic_correction_momentum=diabatic_correction_momentum, + ) + + wind_speed_10percent_below = ( + friction_velocity_below / von_karmans_constant + ) * wind_profile_10percent_below + + complete_wind_profile = np.where( + wind_layer_heights <= (0.1 * canopy_height), + wind_speed_10percent_below, + wind_speed, + ) + + return complete_wind_profile + + +# def update_wind_profile( TODO this is outdated +# wind_heights_above_canopy: NDArray[np.float32], +# wind_heights_below_canopy: NDArray[np.float32], +# temperature: NDArray[np.float32], +# atmospheric_pressure: NDArray[np.float32], +# friction_velocity: NDArray[np.float32], +# canopy_height: NDArray[np.float32], +# leaf_area_index: NDArray[np.float32], +# sensible_heat_topofcanopy: NDArray[np.float32], +# abiotic_const: AbioticConsts, +# ) -> Tuple[DataArray, DataArray]: +# r"""Calculate wind profile above and below canopy. + +# The wind profile above the canopy is described as follows +# (based on :cite:p:`campbell_introduction_1998` as implemented in +# :cite:t:`maclean_microclimc_2021`): + +# :math:`u_z = \frac{u^{*}}{0.4} ln \frac{z-d}{z_M} + \phi_M` + +# where :math:`u_z` is wind speed at height :math:`z` above the canopy, :math:`d` is +# the height above ground within the canopy where the wind profile extrapolates to +# zero, :math:`z_m` the roughness length for momentum, :math:`\phi_M` is a diabatic +# correction for momentum and :math:`u^{*}` is the friction velocity, which gives the +# wind speed at height :math:`d + z_m`. + +# The wind profile below canopy is derived as follows: + +# :math:`u_z = u_h exp(a(\frac{z}{h} - 1))` + +# where :math:`u_z` is wind speed at height :math:`z` within the canopy, :math:`u_h` +# is wind speed at the top of the canopy at height :math:`h`, and :math:`a` is a wind +# attenuation coefficient given by :math:`a = 2 l_m i_w`, where :math:`c_d` is a drag +# coefficient that varies with leaf inclination and shape, :math:`i_w` is a +# coefficient describing relative turbulence intensity and :math:`l_m` is the mean +# mixing length, equivalent to the free space between the leaves and stems. For +# details, see :cite:t:`maclean_microclimc_2021`. + +# Args: +# wind_heights_above_canopy: array of heights above canopy for which wind speed +# is calculated, [m] +# wind_heights_below_canopy: array of heights below canopy for which wind speed +# is calculated, [m], base on 'layer_heights' variable +# temperature: Air temperature, [C] +# atmospheric_pressure: Atmospheric pressure, [kPa] +# friction_velocity: friction velocity, TODO +# canopy_height: Canopy height, [m] +# leaf_area_index: leaf area index, [m m-1] +# sensible_heat_topofcanopy: Sensible heat flux at the top of the canopy, TODO +# abiotic_const: An AbioticConsts instance + +# Returns: +# vertical wind profile above canopy, [m s-1] +# vertical wind profile within canopy, [m s-1] +# """ + +# # preparatory calculations for wind profiles +# zero_plane_displacement = calculate_zero_plane_displacement( +# canopy_height=canopy_height, +# leaf_area_index=leaf_area_index, +# zero_plane_scaling_parameter=abiotic_const.zero_plane_scaling_parameter, +# ) +# roughness_length_momentum = calculate_roughness_length_momentum( +# canopy_height=canopy_height, +# leaf_area_index=leaf_area_index, +# zero_plane_displacement=zero_plane_displacement, +# substrate_surface_drag_coefficient=( +# abiotic_const.substrate_surface_drag_coefficient +# ), +# roughness_element_drag_coefficient=( +# abiotic_const.roughness_element_drag_coefficient +# ), +# roughness_sublayer_depth_parameter=( +# abiotic_const.roughness_sublayer_depth_parameter +# ), +# max_ratio_wind_to_friction_velocity=( +# abiotic_const.max_ratio_wind_to_friction_velocity +# ), +# ) + +# diabatic_correction_momentum_above = calculate_diabatic_correction_momentum_above( +# temperature=temperature, +# atmospheric_pressure=atmospheric_pressure, +# sensible_heat_flux=sensible_heat_topofcanopy, +# friction_velocity=friction_velocity, +# wind_heights=wind_heights_above_canopy, +# zero_plane_displacement=zero_plane_displacement, +# celsius_to_kelvin=abiotic_const.celsius_to_kelvin, +# gravity=abiotic_const.gravity, +# ) + +# mixing_length = calculate_mixing_length( +# canopy_height=canopy_height, +# zero_plane_displacement=zero_plane_displacement, +# roughness_length_momentum=roughness_length_momentum, +# mixing_length_factor=abiotic_const.mixing_length_factor, +# ) + +# wind_attenuation_coefficient = calculate_wind_attenuation_coefficient( +# canopy_height=canopy_height, +# leaf_area_index=leaf_area_index, +# mixing_length=mixing_length, +# drag_coefficient=abiotic_const.drag_coefficient, +# relative_turbulence_intensity=abiotic_const.relative_turbulence_intensity, +# diabatic_correction_factor_below=abiotic_const.diabatic_correction_factor_below, +# ) + +# # Calculate wind profiles and return as Tuple +# wind_profile_above = calculate_wind_above_canopy( +# wind_heights_above_canopy=wind_heights_above_canopy, +# zero_plane_displacement=zero_plane_displacement, +# roughness_length_momentum=roughness_length_momentum, +# diabatic_correction_momentum_above=diabatic_correction_momentum_above, +# friction_velocity=friction_velocity, +# ) + +# wind_profile_canopy = calculate_wind_below_canopy( +# canopy_node_heights=wind_heights_below_canopy, +# wind_profile_above=wind_profile_above, +# wind_attenuation_coefficient=wind_attenuation_coefficient, +# canopy_height=canopy_height, +# ) + +# return wind_profile_above, wind_profile_canopy From eef1daad8a4bb49f41b5c07844a9e92aa0c2a443 Mon Sep 17 00:00:00 2001 From: vgro Date: Wed, 22 Nov 2023 11:57:12 +0000 Subject: [PATCH 02/18] revisited all functions --- tests/models/abiotic/test_wind.py | 276 +++++------- .../models/abiotic/constants.py | 10 + virtual_rainforest/models/abiotic/wind.py | 422 ++++++++---------- 3 files changed, 319 insertions(+), 389 deletions(-) diff --git a/tests/models/abiotic/test_wind.py b/tests/models/abiotic/test_wind.py index 2939652a5..7cade9145 100644 --- a/tests/models/abiotic/test_wind.py +++ b/tests/models/abiotic/test_wind.py @@ -49,10 +49,10 @@ def test_calculate_zero_plane_displacement(dummy_data): from virtual_rainforest.models.abiotic.wind import calculate_zero_plane_displacement - plant_area_index = dummy_data.data["leaf_area_index"].sum(dim="layers") + leaf_area_index = dummy_data.data["leaf_area_index"].sum(dim="layers") result = calculate_zero_plane_displacement( canopy_height=dummy_data.data["canopy_height"].to_numpy(), - plant_area_index=plant_area_index.to_numpy(), + leaf_area_index=leaf_area_index.to_numpy(), zero_plane_scaling_parameter=AbioticConsts.zero_plane_scaling_parameter, ) @@ -66,10 +66,10 @@ def test_calculate_roughness_length_momentum(dummy_data): calculate_roughness_length_momentum, ) - plant_area_index = dummy_data.data["leaf_area_index"].sum(dim="layers") + leaf_area_index = dummy_data.data["leaf_area_index"].sum(dim="layers") result = calculate_roughness_length_momentum( canopy_height=dummy_data.data["canopy_height"].to_numpy(), - plant_area_index=plant_area_index.to_numpy(), + leaf_area_index=leaf_area_index.to_numpy(), zero_plane_displacement=np.array([0, 8, 43]), substrate_surface_drag_coefficient=( AbioticConsts.substrate_surface_drag_coefficient @@ -83,6 +83,7 @@ def test_calculate_roughness_length_momentum(dummy_data): max_ratio_wind_to_friction_velocity=( AbioticConsts.max_ratio_wind_to_friction_velocity ), + min_roughness_length=AbioticConsts.min_roughness_length, von_karman_constant=AbioticConsts.von_karmans_constant, ) @@ -99,62 +100,92 @@ def test_calculate_diabatic_correction_above(dummy_data): ) result = calculate_diabatic_correction_above( + molar_density_air=np.repeat(28.96, 3), + specific_heat_air=np.repeat(1, 3), temperature=dummy_data["air_temperature"].to_numpy(), - atmospheric_pressure=dummy_data["atmospheric_pressure"].to_numpy(), sensible_heat_flux=dummy_data["sensible_heat_flux_topofcanopy"].to_numpy(), friction_velocity=dummy_data["friction_velocity"].to_numpy(), - wind_heights_above_canopy=np.array([1, 15, 50]), + wind_heights=np.array([1, 15, 50]), zero_plane_displacement=np.array([0, 8, 43]), celsius_to_kelvin=AbioticConsts.celsius_to_kelvin, - standard_pressure=AbioticConsts.standard_pressure, - standard_mole=AbioticConsts.standard_mole, - molar_heat_capacity_air=AbioticConsts.molar_heat_capacity_air, - specific_heat_equ_factor_1=AbioticConsts.specific_heat_equ_factor_1, - specific_heat_equ_factor_2=AbioticConsts.specific_heat_equ_factor_2, von_karmans_constant=AbioticConsts.von_karmans_constant, ) - exp_result_h = np.array([7.5154e-05, 6.2562e-04, 3.0957e-04]) - exp_result_m = np.array([4.5093e-05, 3.7534e-04, 1.857e-04]) + exp_result_h = np.array([0.003044, 0.026923, 0.012881]) + exp_result_m = np.array([0.001827, 0.016154, 0.007729]) np.testing.assert_allclose(result["psi_h"], exp_result_h, rtol=1e-4, atol=1e-4) np.testing.assert_allclose(result["psi_m"], exp_result_m, rtol=1e-4, atol=1e-4) -def test_calculate_wind_attenuation_coefficient(dummy_data): - """Test wind attenuation coefficient with and without vegetation.""" +def test_calculate_mean_mixing_length(dummy_data): + """Test mixing length with and without vegetation.""" + + from virtual_rainforest.models.abiotic.wind import calculate_mean_mixing_length + + result = calculate_mean_mixing_length( + canopy_height=dummy_data.data["canopy_height"].to_numpy(), + zero_plane_displacement=np.array([0, 8, 43]), + roughness_length_momentum=np.array([0.003, 0.435, 1.521]), + mixing_length_factor=AbioticConsts.mixing_length_factor, + ) + + np.testing.assert_allclose( + result, np.array([0, 0.419, 1.467]), rtol=1e-3, atol=1e-3 + ) + + +def test_generate_relative_turbulence_intensity(dummy_data): + """Test relative turbulence intensity.""" from virtual_rainforest.models.abiotic.wind import ( - calculate_wind_attenuation_coefficient, + generate_relative_turbulence_intensity, ) - plant_area_index = dummy_data.data["leaf_area_index"].sum(dim="layers") - result = calculate_wind_attenuation_coefficient( - canopy_height=dummy_data.data["canopy_height"].to_numpy(), - plant_area_index=plant_area_index, - mixing_length=np.array([0, 0.4, 1.5]), - drag_coefficient=AbioticConsts.drag_coefficient, - relative_turbulence_intensity=AbioticConsts.relative_turbulence_intensity, - diabatic_correction_factor_below=AbioticConsts.diabatic_correction_factor_below, + result = generate_relative_turbulence_intensity( + layer_heights=dummy_data["wind_layer_heights"].to_numpy(), + min_relative_turbulence_intensity=( + AbioticConsts.min_relative_turbulence_intensity + ), + max_relative_turbulence_intensity=( + AbioticConsts.max_relative_turbulence_intensity + ), + increasing_with_height=True, ) - np.testing.assert_allclose(result, np.array([0.0, 5.91608, 7.302967])) + exp_result = np.array( + [[0.414, 3.06, 8.46], [0.9, 8.46, 13.86], [3.06, 13.86, 30.06]] + ) + np.testing.assert_allclose(result, exp_result) -def test_calculate_mixing_length(dummy_data): - """Test mixing length with and without vegetation.""" +def test_calculate_wind_attenuation_coefficient(dummy_data): + """Test wind attenuation coefficient with and without vegetation.""" - from virtual_rainforest.models.abiotic.wind import calculate_mixing_length + from virtual_rainforest.models.abiotic.wind import ( + calculate_wind_attenuation_coefficient, + ) - result = calculate_mixing_length( + leaf_area_index = dummy_data.data["leaf_area_index"].to_numpy() + relative_turbulence_intensity = np.array( + [[0.414, 3.06, 8.46], [0.9, 8.46, 13.86], [3.06, 13.86, 30.06]] + ) + result = calculate_wind_attenuation_coefficient( canopy_height=dummy_data.data["canopy_height"].to_numpy(), - zero_plane_displacement=np.array([0, 8, 43]), - roughness_length_momentum=np.array([0.003, 0.435, 1.521]), - mixing_length_factor=AbioticConsts.mixing_length_factor, + leaf_area_index=leaf_area_index, + mean_mixing_length=np.array([0, 0.4, 1.5]), + drag_coefficient=AbioticConsts.drag_coefficient, + relative_turbulence_intensity=relative_turbulence_intensity, + diabatic_correction_factor_below=AbioticConsts.diabatic_correction_factor_below, ) - np.testing.assert_allclose( - result, np.array([0, 0.419, 1.467]), rtol=1e-3, atol=1e-3 + exp_result = np.array( + [ + [0.0, 0.9039, 1.4036], + [0.0, 0.7688, 1.0966], + [0.0, 0.7356, 0.7446], + ] ) + np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) def test_wind_log_profile(): @@ -176,143 +207,58 @@ def test_wind_log_profile(): ) -def test_calculate_wind_profile(dummy_data): - """Test if wind profile is calculated correctly for different height.""" +def test_calculate_friction_velocity(dummy_data): + """Calculate friction velocity.""" - from virtual_rainforest.models.abiotic.wind import calculate_wind_profile + from virtual_rainforest.models.abiotic.wind import calculate_fricition_velocity - plant_area_index = dummy_data.data["leaf_area_index"].sum(dim="layers") - result = calculate_wind_profile( - wind_speed_ref=dummy_data["wind_speed_ref"].to_numpy(), - wind_layer_heights=np.array([[5, 12, 60], [0.1, 8, 10]]), + result = calculate_fricition_velocity( + wind_speed_ref=dummy_data["wind_speed_ref"], reference_height=10.0, - attenuation_coefficient=np.array([0, 4, 7]), - plant_area_index=plant_area_index, - canopy_height=dummy_data["canopy_height"].to_numpy(), - diabatic_correction_momentum=np.array([0, 0, 0]), - ground_layer_vegetation_height=0.1, - roughness_sublayer_depth_parameter=( - AbioticConsts.roughness_sublayer_depth_parameter - ), - zero_plane_scaling_parameter=AbioticConsts.zero_plane_scaling_parameter, - substrate_surface_drag_coefficient=( - AbioticConsts.substrate_surface_drag_coefficient - ), - roughness_element_drag_coefficient=( - AbioticConsts.roughness_element_drag_coefficient - ), - max_ratio_wind_to_friction_velocity=( - AbioticConsts.max_ratio_wind_to_friction_velocity - ), + zeroplane_displacement=np.array([0, 8, 43]), + roughness_length_momentum=np.array([0.003, 0.435, 1.521]), + diabatic_correction_momentum=np.array([-0.1, 0.0, 0.1]), von_karmans_constant=AbioticConsts.von_karmans_constant, ) + exp_result = np.array([0.0, 1.310997, 4.0]) + np.testing.assert_allclose(result, exp_result) - np.testing.assert_allclose( - result, - np.array([[0, 7.935657, 24.623732], [0, 1.471923, 0.036979]]), - rtol=1e-3, - atol=1e-3, + +def test_calculate_wind_above_canopy(): + """Wind speed above canopy.""" + + from virtual_rainforest.models.abiotic.wind import calculate_wind_above_canopy + + result = calculate_wind_above_canopy( + friction_velocity=np.array([0, 1.3, 4]), + wind_layer_heights=np.array([2, 12, 52]), + zeroplane_displacement=np.array([0, 8, 43]), + roughness_length_momentum=np.array([0.003, 0.435, 1.521]), + diabatic_correction_momentum=np.array([-0.1, 0.0, 0.1]), + von_karmans_constant=AbioticConsts.von_karmans_constant, + min_wind_speed_above_canopy=AbioticConsts.min_wind_speed_above_canopy, ) + exp_result = np.array([0.0, 7.211, 18.779]) + np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) + + +def test_calculate_wind_canopy(dummy_data): + """Test below canopy wind profile.""" -# def test_calculate_wind_above_canopy(dummy_data): -# """Test wind above canopy.""" - -# from virtual_rainforest.models.abiotic import wind - -# result = wind.calculate_wind_above_canopy( -# wind_heights=DataArray( -# [[10, 20, 30], [20, 30, 40], [40, 50, 60]], dims=["heights", "cell_id"] -# ), -# zero_plane_displacement=DataArray([0, 10, 10], dims="cell_id"), -# roughness_length_momentum=DataArray([0.003, 0.4, 1.5], dims="cell_id"), -# diabatic_correction_momentum_above=DataArray( -# [[1, 1, 1], [1, 1, 1], [1, 1, 1]], dims=["heights", "cell_id"] -# ), -# friction_velocity=dummy_data["friction_velocity"], -# ) - -# xr.testing.assert_allclose( -# result, -# DataArray( -# [ -# [244.3518425, 265.14625792, 285.94067333], -# [41.23594781, 49.90028757, 58.56462732], -# [13.95133583, 15.97866137, 18.53278949], -# ], -# dims=["cell_id", "heights"], -# ), -# ) - - -# def test_calculate_wind_below_canopy(dummy_data): -# """Test wind within canopy.""" - -# from virtual_rainforest.models.abiotic import wind - -# result = wind.calculate_wind_below_canopy( -# canopy_node_heights=dummy_data.data["canopy_node_heights"], -# wind_profile_above=DataArray( -# [ -# [244.0, 265.0, 285.0], -# [41.0, 49.0, 58.0], -# [13.0, 15.0, 18.0], -# ], -# dims=["cell_id", "heights"], -# ), -# wind_attenuation_coefficient=DataArray([0, 0.1, 0.3], dims=["cell_id"]), -# canopy_height=dummy_data.data["canopy_height"], -# ) - -# xr.testing.assert_allclose( -# result, -# DataArray( -# [ -# [0.0, 0.0, 0.0], -# [52.48057, 60.973724, 67.386386], -# [13.334728, 15.492744, 16.450761], -# ], -# dims=["cell_id", "canopy_layers"], -# ), -# ) - - -# def test_calculate_wind_profile(dummy_data): -# """Test wind profile above and within canopy.""" - -# from virtual_rainforest.models.abiotic import wind - -# result = wind.calculate_wind_profile( -# wind_heights=DataArray( -# [[10, 20, 30], [20, 30, 40], [40, 50, 60]], dims=["heights", "cell_id"] -# ), -# canopy_node_heights=dummy_data.data["canopy_node_heights"], -# data=dummy_data, -# ) - -# # check wind above canopy -# xr.testing.assert_allclose( -# result[0], -# DataArray( -# DataArray( -# [ -# [244.351842, 265.146258, 285.940673], -# [46.458135, 54.341054, 62.595567], -# [0.0, 0.0, 13.311866], -# ], -# dims=["cell_id", "heights"], -# ), -# ), -# ) -# # check wind below canopy -# xr.testing.assert_allclose( -# result[1], -# DataArray( -# [ -# [0.0, 0.0, 0.0], -# [5.950523e-02, 2.030195e03, 2.135631e06], -# [6.086929e-03, 2.846548e-01, 1.325216e00], -# ], -# dims=["cell_id", "canopy_layers"], -# ), -# ) + from virtual_rainforest.models.abiotic.wind import calculate_wind_canopy + + result = calculate_wind_canopy( + top_of_canopy_wind_speed=np.array([2, 5, 7]), + wind_layer_heights=dummy_data["wind_layer_heights"], + canopy_height=dummy_data["canopy_height"], + attenuation_coefficient=np.array([0, 0.9, 1.4]), + ) + exp_result = np.array( + [ + [0.0, 1.797693e308, 1.797693e308], + [2.0, 7.841561, 57.16319], + [2.0, 3.188141, 8.051917], + ] + ) + np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) diff --git a/virtual_rainforest/models/abiotic/constants.py b/virtual_rainforest/models/abiotic/constants.py index a4bab440d..25f2d476f 100644 --- a/virtual_rainforest/models/abiotic/constants.py +++ b/virtual_rainforest/models/abiotic/constants.py @@ -60,3 +60,13 @@ class AbioticConsts(ConstantsDataclass): "Diabatic correction factor below canopy." mixing_length_factor: float = 0.32 """Factor in calculation of mixing length :cite:p:`maclean_microclimc_2021`.""" + min_relative_turbulence_intensity: float = 0.36 + """minimum relative turbulence intensity, default value from Shaw et al (1974) + Agricultural Meteorology, 13: 419-425. """ + max_relative_turbulence_intensity: float = 0.9 + """maximum relative turbulence intensity, default value from Shaw et al (1974) + Agricultural Meteorology, 13: 419-425.""" + min_wind_speed_above_canopy: float = 0.55 + """Minimum wind speed above the canopy, [m/s].""" + min_roughness_length: float = 0.05 + """Minimum roughness length, [m].""" diff --git a/virtual_rainforest/models/abiotic/wind.py b/virtual_rainforest/models/abiotic/wind.py index 3836069a8..b7ce33547 100644 --- a/virtual_rainforest/models/abiotic/wind.py +++ b/virtual_rainforest/models/abiotic/wind.py @@ -1,11 +1,34 @@ -r"""The wind module calculates the above- and within-canopy wind profiles for the -Virtual Rainforest. These profiles will determine the exchange of heat, water, and +r"""The wind module calculates the above- and within-canopy wind profile for the +Virtual Rainforest. The wind profile determines the exchange of heat, water, and :math:`CO_{2}` between soil and atmosphere below the canopy as well as the exchange with -the atmsophere above the canopy. +the atmosphere above the canopy. + +The wind profile above the canopy is described as follows (based on +:cite:p:`campbell_introduction_1998` as implemented in +:cite:t:`maclean_microclimc_2021`): + +:math:`u_z = \frac{u^{*}}{0.4} ln \frac{z-d}{z_M} + \phi_M` + +where :math:`u_z` is wind speed at height :math:`z` above the canopy, :math:`d` is +the height above ground within the canopy where the wind profile extrapolates to +zero, :math:`z_m` the roughness length for momentum, :math:`\phi_M` is a diabatic +correction for momentum and :math:`u^{*}` is the friction velocity, which gives the +wind speed at height :math:`d + z_m`. + +The wind profile below canopy is derived as follows: + +:math:`u_z = u_h exp(a(\frac{z}{h} - 1))` + +where :math:`u_z` is wind speed at height :math:`z` within the canopy, :math:`u_h` +is wind speed at the top of the canopy at height :math:`h`, and :math:`a` is a wind +attenuation coefficient given by :math:`a = 2 l_m i_w`, where :math:`c_d` is a drag +coefficient that varies with leaf inclination and shape, :math:`i_w` is a +coefficient describing relative turbulence intensity and :math:`l_m` is the mean +mixing length, equivalent to the free space between the leaves and stems. For +details, see :cite:t:`maclean_microclimc_2021`. + TODO: add sanity checks, errors and logging -TODO check equations for above profile, and howdoes it need adjustment for below -TODO add function that alls all the steps and returns info for data object -TODO select wind layer heights +TODO add function that calls all the steps and returns info for data object """ # noqa: D205, D415 from typing import Union @@ -13,15 +36,10 @@ import numpy as np from numpy.typing import NDArray -from virtual_rainforest.models.abiotic.abiotic_tools import ( - calculate_molar_density_air, - calculate_specific_heat_air, -) - def calculate_zero_plane_displacement( canopy_height: NDArray[np.float32], - plant_area_index: NDArray[np.float32], + leaf_area_index: NDArray[np.float32], zero_plane_scaling_parameter: float, ) -> NDArray[np.float32]: """Calculate zero plane displacement. @@ -33,8 +51,8 @@ def calculate_zero_plane_displacement( :cite:t:`maclean_microclimc_2021`. Args: - canopy_height: canopy height, [m] - plant_area_index: leaf area index, [m m-1] + canopy_height: Canopy height, [m] + leaf_area_index: Total leaf area index, [m m-1] zero_plane_scaling_parameter: Control parameter for scaling d/h :cite:p:`raupach_simplified_1994` @@ -43,12 +61,12 @@ def calculate_zero_plane_displacement( """ # Calculate in presence of vegetation - plant_area_index_displacement = np.where( - plant_area_index > 0, plant_area_index, np.nan + leaf_area_index_displacement = np.where( + leaf_area_index > 0, leaf_area_index, np.nan ) scale_displacement = np.sqrt( - zero_plane_scaling_parameter * plant_area_index_displacement + zero_plane_scaling_parameter * leaf_area_index_displacement ) zero_place_displacement = ( @@ -61,12 +79,13 @@ def calculate_zero_plane_displacement( def calculate_roughness_length_momentum( canopy_height: NDArray[np.float32], - plant_area_index: NDArray[np.float32], + leaf_area_index: NDArray[np.float32], zero_plane_displacement: NDArray[np.float32], substrate_surface_drag_coefficient: float, roughness_element_drag_coefficient: float, roughness_sublayer_depth_parameter: float, max_ratio_wind_to_friction_velocity: float, + min_roughness_length: float, von_karman_constant: float, ) -> NDArray[np.float32]: """Calculate roughness length governing momentum transfer. @@ -77,16 +96,19 @@ def calculate_roughness_length_momentum( effect. Implementation after :cite:t:`maclean_microclimc_2021`. Args: - canopy_height: canopy height, [m] - plant_area_index: plant area index, [m m-1] - substrate_surface_drag_coefficient: substrate-surface drag coefficient - zero_plane_displacement: height above ground within the canopy where the wind + canopy_height: Canopy height, [m] + leaf_area_index: Total leaf area index, [m m-1] + zero_plane_displacement: Height above ground within the canopy where the wind profile extrapolates to zero, [m] - roughness_element_drag_coefficient: roughness-element drag coefficient - roughness_sublayer_depth_parameter: parameter characterizes the roughness + substrate_surface_drag_coefficient: Substrate-surface drag coefficient + roughness_element_drag_coefficient: Roughness-element drag coefficient + roughness_sublayer_depth_parameter: Parameter that characterizes the roughness sublayer depth max_ratio_wind_to_friction_velocity: Maximum ratio of wind velocity to friction velocity + min_roughness_length: Minimum roughness length, [m] + von_karman_constant: Von Karman's constant, unitless constant describing the + logarithmic velocity profile of a turbulent fluid near a no-slip boundary. Returns: momentum roughness length, [m] @@ -95,7 +117,7 @@ def calculate_roughness_length_momentum( # calculate ratio of wind velocity to friction velocity ratio_wind_to_friction_velocity = np.sqrt( substrate_surface_drag_coefficient - + (roughness_element_drag_coefficient * plant_area_index) / 2 + + (roughness_element_drag_coefficient * leaf_area_index) / 2 ) # if the ratio of wind velocity to friction velocity is larger than the set maximum, @@ -120,22 +142,18 @@ def calculate_roughness_length_momentum( initial_roughness_length, ) - return np.where(roughness_length <= 0, 0.01, roughness_length) + return np.where(roughness_length <= 0, min_roughness_length, roughness_length) def calculate_diabatic_correction_above( + molar_density_air: Union[float, NDArray[np.float32]], + specific_heat_air: Union[float, NDArray[np.float32]], temperature: NDArray[np.float32], - atmospheric_pressure: NDArray[np.float32], sensible_heat_flux: NDArray[np.float32], friction_velocity: NDArray[np.float32], - wind_heights_above_canopy: NDArray[np.float32], + wind_heights: NDArray[np.float32], zero_plane_displacement: NDArray[np.float32], celsius_to_kelvin: float, - standard_mole: float, - standard_pressure: float, - molar_heat_capacity_air: float, - specific_heat_equ_factor_1: float, - specific_heat_equ_factor_2: float, von_karmans_constant: float, ) -> dict[str, NDArray[np.float32]]: r"""Calculates the diabatic correction factors for momentum and heat above canopy. @@ -147,45 +165,27 @@ def calculate_diabatic_correction_above( is much cooler than the air above it, it increases to values around 4. Args: + molar_density_air: molar density of air, [mol m-3] + specific_heat_air: specific heat of air, [J mol-1 K-1] temperature: 2m temperature - atmospheric_pressure: atmospheric pressure, [kPa] - sensible_heat_flux: sensible heat flux from canopy to atmosphere above, + sensible_heat_flux: Sensible heat flux from canopy to atmosphere above, [W m-2], # TODO: could be the top entry of the general sensible heat flux - friction_velocity: friction velocity - wind_heights: vector of heights for which wind speed is calculated, [m] - zero_plane_displacement: height above ground within the canopy where the wind + friction_velocity: Friction velocity + wind_heights: Vector of heights for which wind speed is calculated, [m] + zero_plane_displacement: Height above ground within the canopy where the wind profile extrapolates to zero, [m] celsius_to_kelvin: Factor to convert temperature in Celsius to absolute temperature in Kelvin - standard_mole: Moles of ideal gas in 1 m^3 air at standard atmosphere - standard_pressure: Standard atmospheric pressure, [kPa] - molar_heat_capacity_air: molar heat capacity of air, [J mol-1 C-1] - specific_heat_equ_factor_1: Factor in calculation of molar specific heat of air - specific_heat_equ_factor_2: Factor in calculation of molar specific heat of air von_karmans_constant Returns: diabatic correction factors for heat (psi_h) and momentum (psi_m) transfer """ - molar_density_air = calculate_molar_density_air( - temperature=temperature, - atmospheric_pressure=atmospheric_pressure, - standard_mole=standard_mole, - standard_pressure=standard_pressure, - celsius_to_kelvin=celsius_to_kelvin, - ) - specific_heat_air = calculate_specific_heat_air( - temperature=temperature, - molar_heat_capacity_air=molar_heat_capacity_air, - specific_heat_equ_factor_1=specific_heat_equ_factor_1, - specific_heat_equ_factor_2=specific_heat_equ_factor_2, - ) - # calculate atmospheric stability stability = ( von_karmans_constant - * (wind_heights_above_canopy - zero_plane_displacement) + * (wind_heights - zero_plane_displacement) * sensible_heat_flux ) / ( molar_density_air @@ -212,7 +212,7 @@ def calculate_diabatic_correction_above( return {"psi_m": diabatic_correction_momentum, "psi_h": diabatic_correction_heat} -def calculate_mixing_length( +def calculate_mean_mixing_length( canopy_height: NDArray[np.float32], zero_plane_displacement: NDArray[np.float32], roughness_length_momentum: NDArray[np.float32], @@ -220,36 +220,72 @@ def calculate_mixing_length( ) -> NDArray[np.float32]: """Calculate mixing length for canopy air transport. - The mixing length is used to calculate turbulent air transport inside vegetated + The mean mixing length is used to calculate turbulent air transport inside vegetated canopies. It is made equivalent to the above canopy value at the canopy surface. In - absence of vegetation, it is set to zero. - Implementation after :cite:t:`maclean_microclimc_2021`. + absence of vegetation, it is set to zero. Implementation after + :cite:t:`maclean_microclimc_2021`. Args: - canopy_height: canopy height, [m] - leaf_area_index: leaf area index, [m m-1] - zero_plane_displacement: height above ground within the canopy where the wind + canopy_height: Canopy height, [m] + zero_plane_displacement: Height above ground within the canopy where the wind profile extrapolates to zero, [m] - roughness_length_momentum: momentum roughness length, [m] - mixing_length_factor: Factor in calculation of mixing length + roughness_length_momentum: Momentum roughness length, [m] + mixing_length_factor: Factor in calculation of mean mixing length Returns: mixing length for canopy air transport, [m] """ - mixing_length = ( + mean_mixing_length = ( mixing_length_factor * (canopy_height - zero_plane_displacement) ) / np.log((canopy_height - zero_plane_displacement) / roughness_length_momentum) - return np.nan_to_num(mixing_length, nan=0) + return np.nan_to_num(mean_mixing_length, nan=0) + + +def generate_relative_turbulence_intensity( + layer_heights: NDArray[np.float32], + min_relative_turbulence_intensity: float, + max_relative_turbulence_intensity: float, + increasing_with_height: bool, +) -> NDArray[np.float32]: + """Generate relative turbulence intensity profile. + + At the moment, default values are for a maize crop Shaw et al (1974) + Agricultural Meteorology, 13: 419-425. TODO adjust to environment + + Args: + layer_heights: heights of above ground layers, m + min_relative_turbulence_intensity: minimum relative turbulence intensity + max_relative_turbulence_intensity: maximum relative turbulence intensity + increasing_with_height: increasing logical indicating whether turbulence + intensity increases (TRUE) or decreases (FALSE) with height + + Returns: + relative turbulence intensity for each node + """ + + if increasing_with_height: + turbulence_intensity = ( + min_relative_turbulence_intensity + + (max_relative_turbulence_intensity - min_relative_turbulence_intensity) + * layer_heights + ) + else: + turbulence_intensity = ( + max_relative_turbulence_intensity + - (max_relative_turbulence_intensity - min_relative_turbulence_intensity) + * layer_heights + ) + return turbulence_intensity def calculate_wind_attenuation_coefficient( canopy_height: NDArray[np.float32], - plant_area_index: NDArray[np.float32], - mixing_length: NDArray[np.float32], + leaf_area_index: NDArray[np.float32], + mean_mixing_length: NDArray[np.float32], drag_coefficient: float, - relative_turbulence_intensity: float, + relative_turbulence_intensity: NDArray[np.float32], diabatic_correction_factor_below: float, ) -> NDArray[np.float32]: """Calculate wind attenuation coefficient. @@ -259,20 +295,20 @@ def calculate_wind_attenuation_coefficient( Implementation after :cite:t:`maclean_microclimc_2021`. Args: - canopy_height: canopy height, [m] - plant_area_index: plant area index, [m m-1] - mixing_length: mixing length for canopy air transport, [m] - drag_coefficient: drag coefficient - relative_turbulence_intensity: relative turbulence intensity - diabatic_correction_factor_below: diabatic correction factor below canopy + canopy_height: Canopy height, [m] + leaf_area_index: Total leaf area index, [m m-1] + mean_mixing_length: Mixing length for canopy air transport, [m] + drag_coefficient: Drag coefficient + relative_turbulence_intensity: Relative turbulence intensity + diabatic_correction_factor_below: Diabatic correction factor for momentum Returns: wind attenuation coefficient """ intermediate_coefficient1 = ( - (drag_coefficient * plant_area_index * canopy_height) - / (2 * mixing_length * relative_turbulence_intensity), + (drag_coefficient * leaf_area_index * canopy_height) + / (2 * mean_mixing_length * relative_turbulence_intensity), ) intermediate_coefficient2 = np.power(intermediate_coefficient1, 0.5) @@ -293,11 +329,11 @@ def wind_log_profile( """Calculate logarithmic wind profile. Args: - height: array of heights for which wind speed is calculated, [m] - zeroplane_displacement: height above ground within the canopy where the wind + height: Array of heights for which wind speed is calculated, [m] + zeroplane_displacement: Height above ground within the canopy where the wind profile extrapolates to zero, [m] - roughness_length_momentum: momentum roughness length, [m] - diabatic_correction_momentum: diabaric correction factor for momentum + roughness_length_momentum: Momentum roughness length, [m] + diabatic_correction_momentum: Diabatic correction factor for momentum Returns: logarithmic wind profile, [m s-1] @@ -311,158 +347,120 @@ def wind_log_profile( return np.nan_to_num(wind_profile, nan=1).squeeze() -def calculate_wind_profile( +def calculate_fricition_velocity( wind_speed_ref: NDArray[np.float32], reference_height: float, - wind_layer_heights: NDArray[np.float32], - attenuation_coefficient: NDArray[np.float32], - plant_area_index: NDArray[np.float32], - canopy_height: NDArray[np.float32], + zeroplane_displacement: NDArray[np.float32], + roughness_length_momentum: NDArray[np.float32], diabatic_correction_momentum: NDArray[np.float32], - ground_layer_vegetation_height: Union[float, NDArray[np.float32]], - roughness_sublayer_depth_parameter: float, - zero_plane_scaling_parameter: float, - substrate_surface_drag_coefficient: float, - roughness_element_drag_coefficient: float, - max_ratio_wind_to_friction_velocity: float, von_karmans_constant: float, ) -> NDArray[np.float32]: - """Calculate wind speed at any given height from wind speed at reference height. - - Wind profile above the canopy dictates heat and vapor exchange between the canopy - and air above it, and therefore ultimately determine temperature and vapor profiles. - The wind profile above canopy typically follows a logarithmic height profile, which - extrapolates to zero roughly two thirds of the way to the top of the canopy. The - profile itself is thus dependent on the height of the canopy, but also on the - roughness of the vegetation layer, which causes wind shear. We follow the - implementation by :cite:t:`campbell_introduction_1998` as described in - :cite:t:`maclean_microclimc_2021`. + """Calculate friction velocity from wind speed at reference height. Args: wind_speed_ref: Wind speed at reference height, [m s-1] reference_height: Height of wind measurement, [m] - wind_layer_heights: Heights for which wind speed is required, [m] - attenuation_coefficient: Attenuation coefficient as returned by - :func:`~virtual_rainforest.models.abiotic.wind.calculate_wind_attenuation_coefficient` - plant_area_index: Plant area index, [m m-1] - canopy_height: Canopy height, [m] + zero_plane_displacement: Height above ground within the canopy where the wind + profile extrapolates to zero, [m] + roughness_length_momentum: Momentum roughness length, [m] diabatic_correction_momentum: Diabatic correction factor for momentum as returned by :func:`~virtual_rainforest.models.abiotic.wind.calculate_diabatic_correction_above` - ground_layer_vegetation_height: Height of ground vegetation layer below canopy, - [m] - roughness_sublayer_depth_parameter: parameter characterizes the roughness - sublayer depth, [m] - zero_plane_scaling_parameter: Control parameter for scaling zero - displacement/height - substrate_surface_drag_coefficient: Substrate-surface drag coefficient - roughness_element_drag_coefficient: Roughness-element drag coefficient - max_ratio_wind_to_friction_velocity: Maximum ratio of wind velocity to friction - velocity von_karmans_constant: Von Karman's constant, unitless, describes the logarithmic velocity profile of a turbulent fluid near a no-slip boundary Returns: - wind speed at required heights above and below canopy, [m s-1] + friction velocity """ - # Calculate zero-plane displacement - zeroplane_displacement = calculate_zero_plane_displacement( - canopy_height=canopy_height, - plant_area_index=plant_area_index, - zero_plane_scaling_parameter=zero_plane_scaling_parameter, - ) - - # Calculate roughness length for momentum, [m] - roughness_length = calculate_roughness_length_momentum( - canopy_height=canopy_height, - plant_area_index=plant_area_index, - zero_plane_displacement=zeroplane_displacement, - substrate_surface_drag_coefficient=substrate_surface_drag_coefficient, - roughness_element_drag_coefficient=roughness_element_drag_coefficient, - roughness_sublayer_depth_parameter=roughness_sublayer_depth_parameter, - max_ratio_wind_to_friction_velocity=max_ratio_wind_to_friction_velocity, - von_karman_constant=von_karmans_constant, - ) - - roughness_length_ground_vegetation = 0.1 * ground_layer_vegetation_height - - # reference wind profile wind_profile_reference = wind_log_profile( height=reference_height, zeroplane_displacement=zeroplane_displacement, - roughness_length_momentum=roughness_length, + roughness_length_momentum=roughness_length_momentum, diabatic_correction_momentum=diabatic_correction_momentum, ) - friction_velocity = von_karmans_constant * (wind_speed_ref / wind_profile_reference) + return von_karmans_constant * (wind_speed_ref / wind_profile_reference) - wind_profile_above = wind_log_profile( - height=wind_layer_heights, - zeroplane_displacement=zeroplane_displacement, - roughness_length_momentum=roughness_length, - diabatic_correction_momentum=diabatic_correction_momentum, - ) - wind_profile_above[wind_profile_above < 0.55] = 0.55 - wind_profile_below = wind_log_profile( +def calculate_wind_above_canopy( + friction_velocity: NDArray[np.float32], + wind_layer_heights: NDArray[np.float32], + zeroplane_displacement: NDArray[np.float32], + roughness_length_momentum: NDArray[np.float32], + diabatic_correction_momentum: NDArray[np.float32], + von_karmans_constant: float, + min_wind_speed_above_canopy: float, +) -> NDArray[np.float32]: + """Calculate wind speed above canopy from wind speed at reference height. + + Wind speed above the canopy dictates heat and vapor exchange between the canopy + and the air above it, and therefore ultimately determines temperature and vapor + profiles. + The wind profile above canopy typically follows a logarithmic height profile, which + extrapolates to zero roughly two thirds of the way to the top of the canopy. The + profile itself is thus dependent on the height of the canopy, but also on the + roughness of the vegetation layer, which causes wind shear. We follow the + implementation by :cite:t:`campbell_introduction_1998` as described in + :cite:t:`maclean_microclimc_2021`. + + Args: + friction_velocity: friction velocity + wind_layer_heights: Heights for which wind speed is required, [m] + zero_plane_displacement: Height above ground within the canopy where the wind + profile extrapolates to zero, [m] + roughness_length_momentum: Momentum roughness length, [m] + diabatic_correction_momentum: Diabatic correction factor for momentum as + returned by + :func:`~virtual_rainforest.models.abiotic.wind.calculate_diabatic_correction_above` + von_karmans_constant: Von Karman's constant, unitless, describes the logarithmic + velocity profile of a turbulent fluid near a no-slip boundary + + Returns: + wind speed at required heights above canopy, [m s-1] + """ + + wind_profile_above = wind_log_profile( height=wind_layer_heights, zeroplane_displacement=zeroplane_displacement, - roughness_length_momentum=roughness_length, + roughness_length_momentum=roughness_length_momentum, diabatic_correction_momentum=diabatic_correction_momentum, ) - - wind_profile = np.where( - wind_layer_heights > canopy_height, wind_profile_above, wind_profile_below - ) - wind_speed = (friction_velocity / von_karmans_constant) * wind_profile - - # Required height above 10% of canopy hgt, but below top of canopy - wind_10percent_above = wind_speed * np.exp( - attenuation_coefficient * ((wind_layer_heights / canopy_height) - 1) + np.where( + wind_profile_above < min_wind_speed_above_canopy, + min_wind_speed_above_canopy, + wind_profile_above, ) - wind_speed = np.where( - (wind_layer_heights > (0.1 * canopy_height)) - & (wind_layer_heights < canopy_height), - wind_10percent_above, - wind_speed, - ) + return (friction_velocity / von_karmans_constant) * wind_profile_above - # Required height below 10% of canopy hgt - wind_10percent_below = wind_speed * np.exp( - attenuation_coefficient * (((0.1 * canopy_height) / canopy_height) - 1) - ) - wind_profile_reference_below = wind_log_profile( - height=0.1 * canopy_height, - zeroplane_displacement=0.0, - roughness_length_momentum=roughness_length_ground_vegetation, - diabatic_correction_momentum=diabatic_correction_momentum, - ) - - friction_velocity_below = von_karmans_constant * ( - wind_10percent_below / wind_profile_reference_below - ) +def calculate_wind_canopy( + top_of_canopy_wind_speed: NDArray[np.float32], + wind_layer_heights: NDArray[np.float32], + canopy_height: NDArray[np.float32], + attenuation_coefficient: NDArray[np.float32], +) -> NDArray[np.float32]: + """Calculate wind profile for individual canopy layers with variable leaf area. - wind_profile_10percent_below = wind_log_profile( - height=wind_layer_heights, - zeroplane_displacement=zeroplane_displacement, - roughness_length_momentum=roughness_length_ground_vegetation, - diabatic_correction_momentum=diabatic_correction_momentum, - ) + This function can be extended to account for edge distance effects. - wind_speed_10percent_below = ( - friction_velocity_below / von_karmans_constant - ) * wind_profile_10percent_below + Args: + top_of_canopy_wind_speed: Wind speed at top of canopy layer, [m /s] + wind_layer_heights: Height of canopy layer node, [m] + canopy_height: Height to top of canopy layer, [m] + attenuation_coefficient: Attenuation coefficient as returned by + :func:`~virtual_rainforest.models.abiotic.wind.calculate_wind_attenuation_coefficient` - complete_wind_profile = np.where( - wind_layer_heights <= (0.1 * canopy_height), - wind_speed_10percent_below, - wind_speed, - ) + Returns: + wind speed at height of canopy node, [m/s] + """ - return complete_wind_profile + return np.nan_to_num( + top_of_canopy_wind_speed + * np.exp(attenuation_coefficient * ((wind_layer_heights / canopy_height) - 1)), + nan=0, + ).squeeze() # def update_wind_profile( TODO this is outdated @@ -478,30 +476,6 @@ def calculate_wind_profile( # ) -> Tuple[DataArray, DataArray]: # r"""Calculate wind profile above and below canopy. -# The wind profile above the canopy is described as follows -# (based on :cite:p:`campbell_introduction_1998` as implemented in -# :cite:t:`maclean_microclimc_2021`): - -# :math:`u_z = \frac{u^{*}}{0.4} ln \frac{z-d}{z_M} + \phi_M` - -# where :math:`u_z` is wind speed at height :math:`z` above the canopy, :math:`d` is -# the height above ground within the canopy where the wind profile extrapolates to -# zero, :math:`z_m` the roughness length for momentum, :math:`\phi_M` is a diabatic -# correction for momentum and :math:`u^{*}` is the friction velocity, which gives the -# wind speed at height :math:`d + z_m`. - -# The wind profile below canopy is derived as follows: - -# :math:`u_z = u_h exp(a(\frac{z}{h} - 1))` - -# where :math:`u_z` is wind speed at height :math:`z` within the canopy, :math:`u_h` -# is wind speed at the top of the canopy at height :math:`h`, and :math:`a` is a wind -# attenuation coefficient given by :math:`a = 2 l_m i_w`, where :math:`c_d` is a drag -# coefficient that varies with leaf inclination and shape, :math:`i_w` is a -# coefficient describing relative turbulence intensity and :math:`l_m` is the mean -# mixing length, equivalent to the free space between the leaves and stems. For -# details, see :cite:t:`maclean_microclimc_2021`. - # Args: # wind_heights_above_canopy: array of heights above canopy for which wind speed # is calculated, [m] @@ -555,7 +529,7 @@ def calculate_wind_profile( # gravity=abiotic_const.gravity, # ) -# mixing_length = calculate_mixing_length( +# mean_mixing_length = calculate_mean_mixing_length( # canopy_height=canopy_height, # zero_plane_displacement=zero_plane_displacement, # roughness_length_momentum=roughness_length_momentum, @@ -565,7 +539,7 @@ def calculate_wind_profile( # wind_attenuation_coefficient = calculate_wind_attenuation_coefficient( # canopy_height=canopy_height, # leaf_area_index=leaf_area_index, -# mixing_length=mixing_length, +# mean_mixing_length=mean_mixing_length, # drag_coefficient=abiotic_const.drag_coefficient, # relative_turbulence_intensity=abiotic_const.relative_turbulence_intensity, # diabatic_correction_factor_below=abiotic_const.diabatic_correction_factor_below, From 9376d80fa11460fbfb7914b920fcfafddf6728b0 Mon Sep 17 00:00:00 2001 From: vgro Date: Fri, 24 Nov 2023 11:24:35 +0000 Subject: [PATCH 03/18] old todos deleted --- tests/models/abiotic/test_wind.py | 12 +-- virtual_rainforest/models/abiotic/wind.py | 114 +--------------------- 2 files changed, 7 insertions(+), 119 deletions(-) diff --git a/tests/models/abiotic/test_wind.py b/tests/models/abiotic/test_wind.py index 7cade9145..b9a7fa1ce 100644 --- a/tests/models/abiotic/test_wind.py +++ b/tests/models/abiotic/test_wind.py @@ -175,15 +175,10 @@ def test_calculate_wind_attenuation_coefficient(dummy_data): mean_mixing_length=np.array([0, 0.4, 1.5]), drag_coefficient=AbioticConsts.drag_coefficient, relative_turbulence_intensity=relative_turbulence_intensity, - diabatic_correction_factor_below=AbioticConsts.diabatic_correction_factor_below, ) exp_result = np.array( - [ - [0.0, 0.9039, 1.4036], - [0.0, 0.7688, 1.0966], - [0.0, 0.7356, 0.7446], - ] + [[0.0, 0.816, 1.970], [0.0, 0.591, 1.202], [0.0, 0.541, 0.554]] ) np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) @@ -248,11 +243,14 @@ def test_calculate_wind_canopy(dummy_data): from virtual_rainforest.models.abiotic.wind import calculate_wind_canopy + atten_coeff = np.array( + [[0.0, 0.816, 1.970], [0.0, 0.591, 1.202], [0.0, 0.541, 0.554]] + ) result = calculate_wind_canopy( top_of_canopy_wind_speed=np.array([2, 5, 7]), wind_layer_heights=dummy_data["wind_layer_heights"], canopy_height=dummy_data["canopy_height"], - attenuation_coefficient=np.array([0, 0.9, 1.4]), + attenuation_coefficient=atten_coeff, ) exp_result = np.array( [ diff --git a/virtual_rainforest/models/abiotic/wind.py b/virtual_rainforest/models/abiotic/wind.py index b7ce33547..875e11fea 100644 --- a/virtual_rainforest/models/abiotic/wind.py +++ b/virtual_rainforest/models/abiotic/wind.py @@ -28,7 +28,6 @@ details, see :cite:t:`maclean_microclimc_2021`. TODO: add sanity checks, errors and logging -TODO add function that calls all the steps and returns info for data object """ # noqa: D205, D415 from typing import Union @@ -286,7 +285,6 @@ def calculate_wind_attenuation_coefficient( mean_mixing_length: NDArray[np.float32], drag_coefficient: float, relative_turbulence_intensity: NDArray[np.float32], - diabatic_correction_factor_below: float, ) -> NDArray[np.float32]: """Calculate wind attenuation coefficient. @@ -300,24 +298,17 @@ def calculate_wind_attenuation_coefficient( mean_mixing_length: Mixing length for canopy air transport, [m] drag_coefficient: Drag coefficient relative_turbulence_intensity: Relative turbulence intensity - diabatic_correction_factor_below: Diabatic correction factor for momentum Returns: wind attenuation coefficient """ - intermediate_coefficient1 = ( + intermediate_coefficient = ( (drag_coefficient * leaf_area_index * canopy_height) / (2 * mean_mixing_length * relative_turbulence_intensity), ) - intermediate_coefficient2 = np.power(intermediate_coefficient1, 0.5) - - atten_coeff = intermediate_coefficient2 * np.power( - diabatic_correction_factor_below, 0.5 - ) - - return np.nan_to_num(atten_coeff, nan=0).squeeze() + return np.nan_to_num(intermediate_coefficient, nan=0).squeeze() def wind_log_profile( @@ -461,104 +452,3 @@ def calculate_wind_canopy( * np.exp(attenuation_coefficient * ((wind_layer_heights / canopy_height) - 1)), nan=0, ).squeeze() - - -# def update_wind_profile( TODO this is outdated -# wind_heights_above_canopy: NDArray[np.float32], -# wind_heights_below_canopy: NDArray[np.float32], -# temperature: NDArray[np.float32], -# atmospheric_pressure: NDArray[np.float32], -# friction_velocity: NDArray[np.float32], -# canopy_height: NDArray[np.float32], -# leaf_area_index: NDArray[np.float32], -# sensible_heat_topofcanopy: NDArray[np.float32], -# abiotic_const: AbioticConsts, -# ) -> Tuple[DataArray, DataArray]: -# r"""Calculate wind profile above and below canopy. - -# Args: -# wind_heights_above_canopy: array of heights above canopy for which wind speed -# is calculated, [m] -# wind_heights_below_canopy: array of heights below canopy for which wind speed -# is calculated, [m], base on 'layer_heights' variable -# temperature: Air temperature, [C] -# atmospheric_pressure: Atmospheric pressure, [kPa] -# friction_velocity: friction velocity, TODO -# canopy_height: Canopy height, [m] -# leaf_area_index: leaf area index, [m m-1] -# sensible_heat_topofcanopy: Sensible heat flux at the top of the canopy, TODO -# abiotic_const: An AbioticConsts instance - -# Returns: -# vertical wind profile above canopy, [m s-1] -# vertical wind profile within canopy, [m s-1] -# """ - -# # preparatory calculations for wind profiles -# zero_plane_displacement = calculate_zero_plane_displacement( -# canopy_height=canopy_height, -# leaf_area_index=leaf_area_index, -# zero_plane_scaling_parameter=abiotic_const.zero_plane_scaling_parameter, -# ) -# roughness_length_momentum = calculate_roughness_length_momentum( -# canopy_height=canopy_height, -# leaf_area_index=leaf_area_index, -# zero_plane_displacement=zero_plane_displacement, -# substrate_surface_drag_coefficient=( -# abiotic_const.substrate_surface_drag_coefficient -# ), -# roughness_element_drag_coefficient=( -# abiotic_const.roughness_element_drag_coefficient -# ), -# roughness_sublayer_depth_parameter=( -# abiotic_const.roughness_sublayer_depth_parameter -# ), -# max_ratio_wind_to_friction_velocity=( -# abiotic_const.max_ratio_wind_to_friction_velocity -# ), -# ) - -# diabatic_correction_momentum_above = calculate_diabatic_correction_momentum_above( -# temperature=temperature, -# atmospheric_pressure=atmospheric_pressure, -# sensible_heat_flux=sensible_heat_topofcanopy, -# friction_velocity=friction_velocity, -# wind_heights=wind_heights_above_canopy, -# zero_plane_displacement=zero_plane_displacement, -# celsius_to_kelvin=abiotic_const.celsius_to_kelvin, -# gravity=abiotic_const.gravity, -# ) - -# mean_mixing_length = calculate_mean_mixing_length( -# canopy_height=canopy_height, -# zero_plane_displacement=zero_plane_displacement, -# roughness_length_momentum=roughness_length_momentum, -# mixing_length_factor=abiotic_const.mixing_length_factor, -# ) - -# wind_attenuation_coefficient = calculate_wind_attenuation_coefficient( -# canopy_height=canopy_height, -# leaf_area_index=leaf_area_index, -# mean_mixing_length=mean_mixing_length, -# drag_coefficient=abiotic_const.drag_coefficient, -# relative_turbulence_intensity=abiotic_const.relative_turbulence_intensity, -# diabatic_correction_factor_below=abiotic_const.diabatic_correction_factor_below, -# ) - -# # Calculate wind profiles and return as Tuple -# wind_profile_above = calculate_wind_above_canopy( -# wind_heights_above_canopy=wind_heights_above_canopy, -# zero_plane_displacement=zero_plane_displacement, -# roughness_length_momentum=roughness_length_momentum, -# diabatic_correction_momentum_above=diabatic_correction_momentum_above, -# friction_velocity=friction_velocity, -# ) - -# wind_profile_canopy = calculate_wind_below_canopy( -# canopy_node_heights=wind_heights_below_canopy, -# wind_profile_above=wind_profile_above, -# wind_attenuation_coefficient=wind_attenuation_coefficient, -# canopy_height=canopy_height, -# ) - -# return wind_profile_above, wind_profile_canopy From 4b01dec625a0a93b75a07a2b4401edc6dffe803d Mon Sep 17 00:00:00 2001 From: vgro Date: Fri, 24 Nov 2023 11:27:58 +0000 Subject: [PATCH 04/18] docstring typo fixed --- virtual_rainforest/models/hydrology/above_ground.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/virtual_rainforest/models/hydrology/above_ground.py b/virtual_rainforest/models/hydrology/above_ground.py index bb0dedc41..11741a423 100644 --- a/virtual_rainforest/models/hydrology/above_ground.py +++ b/virtual_rainforest/models/hydrology/above_ground.py @@ -30,7 +30,7 @@ def calculate_soil_evaporation( heat_transfer_coefficient: float, extinction_coefficient_global_radiation: float, ) -> NDArray[np.float32]: - r"""Calculate soil evaporation based classical bulk aerodynamic formulation. + r"""Calculate soil evaporation based on classical bulk aerodynamic formulation. This function uses the so-called 'alpha' method to estimate the evaporative flux :cite:p:`mahfouf_comparative_1991`. @@ -51,7 +51,7 @@ def calculate_soil_evaporation( :math:`E_{act} = E_{g} * exp(-\kappa_{gb}*LAI)` - where :math:`kappa_{gb}` is the extinction coefficient for global radiation, and + where :math:`\kappa_{gb}` is the extinction coefficient for global radiation, and :math:`LAI` is the total leaf area index. Args: From bea66558031ab06e98a458938c51d83b36039315 Mon Sep 17 00:00:00 2001 From: Taran Rallings Date: Fri, 24 Nov 2023 14:56:56 +0000 Subject: [PATCH 05/18] Constants rework into dataclass format. --- tests/models/animals/conftest.py | 39 +++++-- tests/models/animals/test_animal_cohorts.py | 54 ++++++---- .../models/animals/test_animal_communities.py | 33 +++--- tests/models/animals/test_animal_model.py | 16 ++- tests/models/animals/test_functional_group.py | 12 ++- .../models/animals/animal_cohorts.py | 17 +-- .../models/animals/animal_communities.py | 23 ++-- .../models/animals/animal_model.py | 10 +- .../models/animals/constants.py | 100 +++++++++++++++++- .../models/animals/functional_group.py | 42 ++++---- .../models/animals/plant_resources.py | 6 +- 11 files changed, 262 insertions(+), 90 deletions(-) diff --git a/tests/models/animals/conftest.py b/tests/models/animals/conftest.py index da188ca31..0c27262be 100644 --- a/tests/models/animals/conftest.py +++ b/tests/models/animals/conftest.py @@ -77,31 +77,49 @@ def plant_climate_data_instance(layer_roles_fixture): @pytest.fixture -def functional_group_list_instance(shared_datadir): +def constants_instance(): + """Fixture for an instance of animal constants.""" + from virtual_rainforest.models.animals.constants import AnimalConsts + + return AnimalConsts() + + +@pytest.fixture +def functional_group_list_instance(shared_datadir, constants_instance): """Fixture for an animal functional group used in tests.""" from virtual_rainforest.models.animals.functional_group import ( import_functional_groups, ) file = shared_datadir / "example_functional_group_import.csv" - fg_list = import_functional_groups(file) + fg_list = import_functional_groups(file, constants_instance) return fg_list @pytest.fixture -def animal_model_instance(data_instance, functional_group_list_instance): +def animal_model_instance( + data_instance, functional_group_list_instance, constants_instance +): """Fixture for an animal model object used in tests.""" from pint import Quantity from virtual_rainforest.models.animals.animal_model import AnimalModel - return AnimalModel(data_instance, Quantity("1 day"), functional_group_list_instance) + return AnimalModel( + data_instance, + Quantity("1 day"), + functional_group_list_instance, + constants_instance, + ) @pytest.fixture def animal_community_instance( - functional_group_list_instance, animal_model_instance, plant_data_instance + functional_group_list_instance, + animal_model_instance, + plant_data_instance, + constants_instance, ): """Fixture for an animal community used in tests.""" from virtual_rainforest.models.animals.animal_communities import AnimalCommunity @@ -112,28 +130,31 @@ def animal_community_instance( community_key=4, neighbouring_keys=[1, 3, 5, 7], get_destination=animal_model_instance.get_community_by_key, + constants=constants_instance, ) @pytest.fixture -def herbivore_functional_group_instance(shared_datadir): +def herbivore_functional_group_instance(shared_datadir, constants_instance): """Fixture for an animal functional group used in tests.""" from virtual_rainforest.models.animals.functional_group import ( import_functional_groups, ) file = shared_datadir / "example_functional_group_import.csv" - fg_list = import_functional_groups(file) + fg_list = import_functional_groups(file, constants_instance) return fg_list[3] @pytest.fixture -def herbivore_cohort_instance(herbivore_functional_group_instance): +def herbivore_cohort_instance(herbivore_functional_group_instance, constants_instance): """Fixture for an animal cohort used in tests.""" from virtual_rainforest.models.animals.animal_cohorts import AnimalCohort - return AnimalCohort(herbivore_functional_group_instance, 10000.0, 1, 10) + return AnimalCohort( + herbivore_functional_group_instance, 10000.0, 1, 10, constants_instance + ) @pytest.fixture diff --git a/tests/models/animals/test_animal_cohorts.py b/tests/models/animals/test_animal_cohorts.py index 343c1091b..6fbae3c54 100644 --- a/tests/models/animals/test_animal_cohorts.py +++ b/tests/models/animals/test_animal_cohorts.py @@ -5,53 +5,59 @@ @pytest.fixture -def predator_functional_group_instance(shared_datadir): +def predator_functional_group_instance(shared_datadir, constants_instance): """Fixture for an animal functional group used in tests.""" from virtual_rainforest.models.animals.functional_group import ( import_functional_groups, ) file = shared_datadir / "example_functional_group_import.csv" - fg_list = import_functional_groups(file) + fg_list = import_functional_groups(file, constants_instance) return fg_list[2] @pytest.fixture -def predator_cohort_instance(predator_functional_group_instance): +def predator_cohort_instance(predator_functional_group_instance, constants_instance): """Fixture for an animal cohort used in tests.""" from virtual_rainforest.models.animals.animal_cohorts import AnimalCohort - return AnimalCohort(predator_functional_group_instance, 10000.0, 1, 10) + return AnimalCohort( + predator_functional_group_instance, 10000.0, 1, 10, constants_instance + ) @pytest.fixture -def ectotherm_functional_group_instance(shared_datadir): +def ectotherm_functional_group_instance(shared_datadir, constants_instance): """Fixture for an animal functional group used in tests.""" from virtual_rainforest.models.animals.functional_group import ( import_functional_groups, ) file = shared_datadir / "example_functional_group_import.csv" - fg_list = import_functional_groups(file) + fg_list = import_functional_groups(file, constants_instance) return fg_list[5] @pytest.fixture -def ectotherm_cohort_instance(ectotherm_functional_group_instance): +def ectotherm_cohort_instance(ectotherm_functional_group_instance, constants_instance): """Fixture for an animal cohort used in tests.""" from virtual_rainforest.models.animals.animal_cohorts import AnimalCohort - return AnimalCohort(ectotherm_functional_group_instance, 100.0, 1, 10) + return AnimalCohort( + ectotherm_functional_group_instance, 100.0, 1, 10, constants_instance + ) @pytest.fixture -def prey_cohort_instance(herbivore_functional_group_instance): +def prey_cohort_instance(herbivore_functional_group_instance, constants_instance): """Fixture for an animal cohort used in tests.""" from virtual_rainforest.models.animals.animal_cohorts import AnimalCohort - return AnimalCohort(herbivore_functional_group_instance, 100.0, 1, 10) + return AnimalCohort( + herbivore_functional_group_instance, 100.0, 1, 10, constants_instance + ) @pytest.fixture @@ -85,6 +91,7 @@ def test_invalid_animal_cohort_initialization( age, individuals, error_type, + constants_instance, ): """Test for invalid inputs during AnimalCohort initialization.""" from virtual_rainforest.models.animals.animal_cohorts import AnimalCohort @@ -95,6 +102,7 @@ def test_invalid_animal_cohort_initialization( mass, age, individuals, + constants_instance, ) @pytest.mark.parametrize( @@ -346,38 +354,46 @@ def test_eat(self, herbivore_cohort_instance, mocker): ) assert herbivore_cohort_instance.mass_current == 0 - def test_is_below_mass_threshold(self, herbivore_cohort_instance): + def test_is_below_mass_threshold( + self, herbivore_cohort_instance, constants_instance + ): """Test the can_reproduce method of AnimalCohort.""" - from virtual_rainforest.models.animals.constants import BIRTH_MASS_THRESHOLD # TODO: test other mass thresholds # 1. Test when stored_energy is exactly equal to the threshold herbivore_cohort_instance.mass_current = ( - herbivore_cohort_instance.functional_group.adult_mass * BIRTH_MASS_THRESHOLD + herbivore_cohort_instance.functional_group.adult_mass + * constants_instance.birth_mass_threshold ) assert not herbivore_cohort_instance.is_below_mass_threshold( - BIRTH_MASS_THRESHOLD + constants_instance.birth_mass_threshold ) # 2. Test when stored_energy is just below the threshold herbivore_cohort_instance.mass_current = ( - herbivore_cohort_instance.functional_group.adult_mass * BIRTH_MASS_THRESHOLD + herbivore_cohort_instance.functional_group.adult_mass + * constants_instance.birth_mass_threshold - 0.01 ) - assert herbivore_cohort_instance.is_below_mass_threshold(BIRTH_MASS_THRESHOLD) + assert herbivore_cohort_instance.is_below_mass_threshold( + constants_instance.birth_mass_threshold + ) # 3. Test when stored_energy is above the threshold herbivore_cohort_instance.mass_current = ( - herbivore_cohort_instance.functional_group.adult_mass * BIRTH_MASS_THRESHOLD + herbivore_cohort_instance.functional_group.adult_mass + * constants_instance.birth_mass_threshold + 0.01 ) assert not herbivore_cohort_instance.is_below_mass_threshold( - BIRTH_MASS_THRESHOLD + constants_instance.birth_mass_threshold ) # 4. Test with stored_energy set to 0 herbivore_cohort_instance.mass_current = 0.0 - assert herbivore_cohort_instance.is_below_mass_threshold(BIRTH_MASS_THRESHOLD) + assert herbivore_cohort_instance.is_below_mass_threshold( + constants_instance.birth_mass_threshold + ) @pytest.mark.parametrize( "initial_individuals, number_days, mortality_prob", diff --git a/tests/models/animals/test_animal_communities.py b/tests/models/animals/test_animal_communities.py index 209f9e28f..7a033afd4 100644 --- a/tests/models/animals/test_animal_communities.py +++ b/tests/models/animals/test_animal_communities.py @@ -5,7 +5,10 @@ @pytest.fixture def animal_community_destination_instance( - functional_group_list_instance, animal_model_instance, plant_data_instance + functional_group_list_instance, + animal_model_instance, + plant_data_instance, + constants_instance, ): """Fixture for an animal community used in tests.""" from virtual_rainforest.models.animals.animal_communities import AnimalCommunity @@ -16,29 +19,34 @@ def animal_community_destination_instance( community_key=4, neighbouring_keys=[1, 3, 5, 7], get_destination=animal_model_instance.get_community_by_key, + constants=constants_instance, ) @pytest.fixture -def functional_group_instance(shared_datadir): +def functional_group_instance(shared_datadir, constants_instance): """Fixture for an animal functional group used in tests.""" from virtual_rainforest.models.animals.functional_group import ( import_functional_groups, ) file = shared_datadir / "example_functional_group_import.csv" - fg_list = import_functional_groups(file) + fg_list = import_functional_groups(file, constants_instance) return fg_list[3] @pytest.fixture -def animal_cohort_instance(functional_group_instance): +def animal_cohort_instance(functional_group_instance, constants_instance): """Fixture for an animal cohort used in tests.""" from virtual_rainforest.models.animals.animal_cohorts import AnimalCohort return AnimalCohort( - functional_group_instance, functional_group_instance.adult_mass, 1.0, 10 + functional_group_instance, + functional_group_instance.adult_mass, + 1.0, + 10, + constants_instance, ) @@ -164,9 +172,10 @@ def test_die_cohort(self, animal_cohort_instance, animal_community_instance): not in animal_community_instance.animal_cohorts["herbivorous_mammal"] ) - def test_birth(self, animal_community_instance, animal_cohort_instance): + def test_birth( + self, animal_community_instance, animal_cohort_instance, constants_instance + ): """Test the birth method in AnimalCommunity.""" - from virtual_rainforest.models.animals.constants import BIRTH_MASS_THRESHOLD # Setup initial conditions parent_cohort_name = animal_cohort_instance.name @@ -179,7 +188,8 @@ def test_birth(self, animal_community_instance, animal_cohort_instance): # Set the reproductive mass of the parent cohort to ensure it can reproduce required_mass_for_birth = ( - animal_cohort_instance.functional_group.adult_mass * BIRTH_MASS_THRESHOLD + animal_cohort_instance.functional_group.adult_mass + * constants_instance.birth_mass_threshold - animal_cohort_instance.functional_group.adult_mass ) @@ -198,13 +208,11 @@ def test_birth(self, animal_community_instance, animal_cohort_instance): # 2. Check that the reproductive mass of the parent cohort is reduced to 0 assert animal_cohort_instance.reproductive_mass == 0 - def test_birth_community(self, animal_community_instance): + def test_birth_community(self, animal_community_instance, constants_instance): """Test the thresholding behavior of birth_community.""" from itertools import chain - from virtual_rainforest.models.animals.constants import BIRTH_MASS_THRESHOLD - # Preparation: populate the community animal_community_instance.populate_community() @@ -216,7 +224,8 @@ def test_birth_community(self, animal_community_instance): # Set mass to just below the threshold threshold_mass = ( - initial_cohort.functional_group.adult_mass * BIRTH_MASS_THRESHOLD + initial_cohort.functional_group.adult_mass + * constants_instance.birth_mass_threshold - initial_cohort.functional_group.adult_mass ) diff --git a/tests/models/animals/test_animal_model.py b/tests/models/animals/test_animal_model.py index 9b5bb2469..d0a18ae63 100644 --- a/tests/models/animals/test_animal_model.py +++ b/tests/models/animals/test_animal_model.py @@ -11,7 +11,10 @@ def test_animal_model_initialization( - caplog, plant_climate_data_instance, functional_group_list_instance + caplog, + plant_climate_data_instance, + functional_group_list_instance, + constants_instance, ): """Test `AnimalModel` initialization.""" from virtual_rainforest.core.base_model import BaseModel @@ -22,6 +25,7 @@ def test_animal_model_initialization( plant_climate_data_instance, pint.Quantity("1 week"), functional_group_list_instance, + constants=constants_instance, ) # In cases where it passes then checks that the object has the right properties @@ -108,6 +112,7 @@ def test_generate_animal_model( time_interval, raises, expected_log_entries, + constants_instance, ): """Test that the function to initialise the animal model behaves as expected.""" from virtual_rainforest.models.animals.animal_model import AnimalModel @@ -115,9 +120,9 @@ def test_generate_animal_model( # Check whether model is initialised (or not) as expected with raises: model = AnimalModel.from_config( - plant_climate_data_instance, - config, - pint.Quantity(config["core"]["timing"]["update_interval"]), + data=plant_climate_data_instance, + config=config, + update_interval=pint.Quantity(config["core"]["timing"]["update_interval"]), ) assert model.update_interval == time_interval # Run the update step (once this does something should check output) @@ -152,7 +157,7 @@ def test_get_community_by_key(animal_model_instance): def test_update_method_sequence( - plant_climate_data_instance, functional_group_list_instance + plant_climate_data_instance, functional_group_list_instance, constants_instance ): """Test update to ensure it runs the community methods in order. @@ -166,6 +171,7 @@ def test_update_method_sequence( plant_climate_data_instance, pint.Quantity("1 week"), functional_group_list_instance, + constants=constants_instance, ) # Mock all the methods that are supposed to be called by update diff --git a/tests/models/animals/test_functional_group.py b/tests/models/animals/test_functional_group.py index c89bb1794..79cf752ad 100644 --- a/tests/models/animals/test_functional_group.py +++ b/tests/models/animals/test_functional_group.py @@ -99,10 +99,17 @@ def test_initialization( MetabolicType, TaxaType, ) + from virtual_rainforest.models.animals.constants import AnimalConsts from virtual_rainforest.models.animals.functional_group import FunctionalGroup func_group = FunctionalGroup( - name, taxa, diet, metabolic_type, birth_mass, adult_mass + name, + taxa, + diet, + metabolic_type, + birth_mass, + adult_mass, + constants=AnimalConsts(), ) assert func_group.name == name assert func_group.taxa == TaxaType(taxa) @@ -133,13 +140,14 @@ def test_import_functional_groups( MetabolicType, TaxaType, ) + from virtual_rainforest.models.animals.constants import AnimalConsts from virtual_rainforest.models.animals.functional_group import ( FunctionalGroup, import_functional_groups, ) file = shared_datadir / "example_functional_group_import.csv" - fg_list = import_functional_groups(file) + fg_list = import_functional_groups(file, constants=AnimalConsts()) assert len(fg_list) == 6 assert isinstance(fg_list[index], FunctionalGroup) assert fg_list[index].name == name diff --git a/virtual_rainforest/models/animals/animal_cohorts.py b/virtual_rainforest/models/animals/animal_cohorts.py index 375d5e770..33facb677 100644 --- a/virtual_rainforest/models/animals/animal_cohorts.py +++ b/virtual_rainforest/models/animals/animal_cohorts.py @@ -16,11 +16,7 @@ from virtual_rainforest.core.logger import LOGGER from virtual_rainforest.models.animals.animal_traits import DietType -from virtual_rainforest.models.animals.constants import ( - DECAY_FRACTION_CARCASSES, - DECAY_FRACTION_EXCREMENT, - FLOW_TO_REPRODUCTIVE_MASS_THRESHOLD, -) +from virtual_rainforest.models.animals.constants import AnimalConsts from virtual_rainforest.models.animals.decay import CarcassPool from virtual_rainforest.models.animals.functional_group import FunctionalGroup from virtual_rainforest.models.animals.protocols import Consumer, DecayPool, Resource @@ -42,6 +38,7 @@ def __init__( mass: float, age: float, individuals: int, + constants: AnimalConsts = AnimalConsts(), ) -> None: if age < 0: raise ValueError("Age must be a positive number.") @@ -63,6 +60,8 @@ def __init__( """The age of the animal cohort [days].""" self.individuals = individuals """The number of individuals in this cohort.""" + self.constants = constants + """Animal constants.""" self.damuth_density: int = damuths_law( self.functional_group.adult_mass, self.functional_group.damuths_law_terms ) @@ -92,9 +91,9 @@ def __init__( # TODO - In future this should be parameterised using a constants dataclass, but # this hasn't yet been implemented for the animal model - self.decay_fraction_excrement: float = DECAY_FRACTION_EXCREMENT + self.decay_fraction_excrement: float = self.constants.decay_fraction_excrement """The fraction of excrement which decays before it gets consumed.""" - self.decay_fraction_carcasses: float = DECAY_FRACTION_CARCASSES + self.decay_fraction_carcasses: float = self.constants.decay_fraction_carcasses """The fraction of carcass biomass which decays before it gets consumed.""" def metabolize(self, temperature: float, dt: timedelta64) -> None: @@ -326,7 +325,9 @@ def eat(self, food: Resource, pool: DecayPool) -> float: # get the per-individual energetic gain from the bulk value mass_consumed = food.get_eaten(self, pool) / self.individuals - if self.is_below_mass_threshold(FLOW_TO_REPRODUCTIVE_MASS_THRESHOLD): + if self.is_below_mass_threshold( + self.constants.flow_to_reproductive_mass_threshold + ): # if current mass equals or exceeds standard adult mass, gains to repro mass self.mass_current += mass_consumed else: diff --git a/virtual_rainforest/models/animals/animal_communities.py b/virtual_rainforest/models/animals/animal_communities.py index 94f2e9579..8c6b78abc 100644 --- a/virtual_rainforest/models/animals/animal_communities.py +++ b/virtual_rainforest/models/animals/animal_communities.py @@ -18,10 +18,7 @@ from virtual_rainforest.core.data import Data from virtual_rainforest.core.logger import LOGGER from virtual_rainforest.models.animals.animal_cohorts import AnimalCohort -from virtual_rainforest.models.animals.constants import ( - BIRTH_MASS_THRESHOLD, - DISPERSAL_MASS_THRESHOLD, -) +from virtual_rainforest.models.animals.constants import AnimalConsts from virtual_rainforest.models.animals.decay import CarcassPool, ExcrementPool from virtual_rainforest.models.animals.functional_group import FunctionalGroup from virtual_rainforest.models.animals.plant_resources import PlantResources @@ -47,6 +44,7 @@ def __init__( community_key: int, neighbouring_keys: list[int], get_destination: Callable[[int], "AnimalCommunity"], + constants: AnimalConsts = AnimalConsts(), ) -> None: # The constructor of the AnimalCommunity class. self.data = data @@ -59,6 +57,8 @@ def __init__( """List of integer keys of neighbouring communities.""" self.get_destination = get_destination """Callable get_destination from AnimalModel.""" + self.constants = constants + """Animal constants.""" self.animal_cohorts: dict[str, list[AnimalCohort]] = { k.name: [] for k in self.functional_groups @@ -96,7 +96,11 @@ def populate_community(self) -> None: ) cohort = AnimalCohort( - functional_group, functional_group.adult_mass, 0.0, individuals + functional_group, + functional_group.adult_mass, + 0.0, + individuals, + self.constants, ) self.animal_cohorts[functional_group.name].append(cohort) @@ -121,7 +125,7 @@ def migrate(self, migrant: AnimalCohort, destination: AnimalCommunity) -> None: def migrate_community(self) -> None: """This handles migrating all cohorts in a community.""" for cohort in self.all_animal_cohorts: - if cohort.is_below_mass_threshold(DISPERSAL_MASS_THRESHOLD): + if cohort.is_below_mass_threshold(self.constants.dispersal_mass_threshold): # Random walk destination from the neighbouring keys destination_key = choice(self.neighbouring_keys) destination = self.get_destination(destination_key) @@ -177,6 +181,7 @@ def birth(self, parent_cohort: AnimalCohort) -> None: parent_cohort.functional_group.birth_mass, 0.0, number_offspring, + self.constants, ) ) @@ -185,7 +190,7 @@ def birth_community(self) -> None: # reproduction occurs for cohorts with sufficient reproductive mass for cohort in self.all_animal_cohorts: - if not cohort.is_below_mass_threshold(BIRTH_MASS_THRESHOLD): + if not cohort.is_below_mass_threshold(self.constants.birth_mass_threshold): self.birth(cohort) def forage_community(self) -> None: @@ -203,7 +208,9 @@ def forage_community(self) -> None: """ # Generate the plant resources for foraging. plant_community: PlantResources = PlantResources( - data=self.data, cell_id=self.community_key + data=self.data, + cell_id=self.community_key, + constants=self.constants, ) plant_list = [plant_community] diff --git a/virtual_rainforest/models/animals/animal_model.py b/virtual_rainforest/models/animals/animal_model.py index e6de497fa..9aef7bf0e 100644 --- a/virtual_rainforest/models/animals/animal_model.py +++ b/virtual_rainforest/models/animals/animal_model.py @@ -31,6 +31,7 @@ from virtual_rainforest.core.data import Data from virtual_rainforest.core.logger import LOGGER from virtual_rainforest.models.animals.animal_communities import AnimalCommunity +from virtual_rainforest.models.animals.constants import AnimalConsts from virtual_rainforest.models.animals.functional_group import FunctionalGroup @@ -69,6 +70,7 @@ def __init__( data: Data, update_interval: Quantity, functional_groups: list[FunctionalGroup], + constants: AnimalConsts = AnimalConsts(), **kwargs: Any, ): super().__init__(data, update_interval, **kwargs) @@ -82,6 +84,8 @@ def __init__( self.communities: dict[int, AnimalCommunity] = {} """Set empty dict for populating with communities.""" + self.constants = constants + """Animal constants.""" self._initialize_communities(functional_groups) """Create the dictionary of animal communities and populate each community with animal cohorts.""" @@ -126,6 +130,7 @@ def _initialize_communities(self, functional_groups: list[FunctionalGroup]) -> N community_key=k, neighbouring_keys=list(self.data.grid.neighbours[k]), get_destination=self.get_community_by_key, + constants=self.constants, ) for k in self.data.grid.cell_id } @@ -155,7 +160,10 @@ def from_config( functional_groups_raw = config["animals"]["functional_groups"] - animal_functional_groups = [FunctionalGroup(**k) for k in functional_groups_raw] + animal_functional_groups = [ + FunctionalGroup(**k, constants=AnimalConsts()) + for k in functional_groups_raw + ] LOGGER.info( "Information required to initialise the animal model successfully " diff --git a/virtual_rainforest/models/animals/constants.py b/virtual_rainforest/models/animals/constants.py index dbd156c4c..1b599ea90 100644 --- a/virtual_rainforest/models/animals/constants.py +++ b/virtual_rainforest/models/animals/constants.py @@ -11,6 +11,7 @@ from dataclasses import dataclass, field +from virtual_rainforest.core.constants_class import ConstantsDataclass from virtual_rainforest.models.animals.animal_traits import ( DietType, MetabolicType, @@ -19,7 +20,7 @@ @dataclass(frozen=True) -class AnimalConsts: +class AnimalConsts(ConstantsDataclass): """Dataclass to store all constants related to metabolic rates. TODO: The entire constants fille will be reworked in this style after the energy to @@ -41,7 +42,99 @@ class AnimalConsts: } ) + damuths_law_terms: dict[TaxaType, dict[DietType, tuple[float, float]]] = field( + default_factory=lambda: { + TaxaType.MAMMAL: { + DietType.HERBIVORE: (-0.75, 4.23), + DietType.CARNIVORE: (-0.75, 1.00), + }, + TaxaType.BIRD: { + DietType.HERBIVORE: (-0.75, 5.00), + DietType.CARNIVORE: (-0.75, 2.00), + }, + TaxaType.INSECT: { + DietType.HERBIVORE: (-0.75, 5.00), + DietType.CARNIVORE: (-0.75, 2.00), + }, + } + ) + fat_mass_terms: dict[TaxaType, tuple[float, float]] = field( + default_factory=lambda: { + TaxaType.MAMMAL: (1.19, 0.02), # Scaling of mammalian herbivore fat mass + TaxaType.BIRD: (1.19, 0.05), # Toy Values + TaxaType.INSECT: (1.19, 0.05), # Toy Values + } + ) + + muscle_mass_terms: dict[TaxaType, tuple[float, float]] = field( + default_factory=lambda: { + TaxaType.MAMMAL: (1.0, 0.38), # Scaling of mammalian herbivore muscle mass + TaxaType.BIRD: (1.0, 0.40), # Toy Values + TaxaType.INSECT: (1.0, 0.40), # Toy Values + } + ) + + intake_rate_terms: dict[TaxaType, tuple[float, float]] = field( + default_factory=lambda: { + TaxaType.MAMMAL: (0.71, 0.63), # Mammalian maximum intake rate + TaxaType.BIRD: (0.7, 0.50), # Toy Values + TaxaType.INSECT: (0.7, 0.50), # Toy Values + } + ) + + energy_density: dict[str, float] = field( + default_factory=lambda: { + "meat": 7000.0, # Energy of mammal meat [J/g] + "plant": 18200000.0, # Energy of plant food [J/g] + } + ) + + conversion_efficiency: dict[DietType, float] = field( + default_factory=lambda: { + DietType.HERBIVORE: 0.1, # Toy value + DietType.CARNIVORE: 0.25, # Toy value + } + ) + + mechanical_efficiency: dict[DietType, float] = field( + default_factory=lambda: { + DietType.HERBIVORE: 0.9, # Toy value + DietType.CARNIVORE: 0.8, # Toy value + } + ) + + prey_mass_scaling_terms: dict[ + MetabolicType, dict[TaxaType, tuple[float, float]] + ] = field( + default_factory=lambda: { + MetabolicType.ENDOTHERMIC: { + TaxaType.MAMMAL: (1.0, 1.0), # Toy values + TaxaType.BIRD: (1.0, 1.0), # Toy values + }, + MetabolicType.ECTOTHERMIC: {TaxaType.INSECT: (1.0, 1.0)}, # Toy values + } + ) + + longevity_scaling_terms: dict[TaxaType, tuple[float, float]] = field( + default_factory=lambda: { + TaxaType.MAMMAL: (0.25, 0.02), # Toy values + TaxaType.BIRD: (0.25, 0.05), # Toy values + TaxaType.INSECT: (0.25, 0.05), # Toy values + } + ) + + birth_mass_threshold: float = 1.5 # Threshold for reproduction + flow_to_reproductive_mass_threshold: float = ( + 1.0 # Threshold of trophic flow to reproductive mass + ) + dispersal_mass_threshold: float = 0.75 # Threshold for dispersal + energy_percentile_threshold: float = 0.5 # Threshold for initiating migration + decay_fraction_excrement: float = 0.5 # Decay fraction for excrement + decay_fraction_carcasses: float = 0.2 # Decay fraction for carcasses + + +""" METABOLIC_RATE_TERMS: dict[MetabolicType, dict[str, tuple[float, float]]] = { # Parameters from Madingley, mass based metabolic rates MetabolicType.ENDOTHERMIC: { @@ -161,7 +254,7 @@ class AnimalConsts: DISPERSAL_MASS_THRESHOLD: float = 0.75 # Toy value for thesholding dispersal. ENERGY_PERCENTILE_THRESHOLD: float = 0.5 # Toy value for initiating migration - +""" DECAY_FRACTION_EXCREMENT: float = 0.5 """Fraction of excrement that is assumed to decay rather than be consumed [unitless]. @@ -176,3 +269,6 @@ class AnimalConsts: [unitless]. TODO - The number given here is very much made up, see :attr:`DECAY_FRACTION_EXCREMENT` for details of how this should be changed in future. """ +BOLTZMANN_CONSTANT: float = 8.617333262145e-5 # Boltzmann constant [eV/K] + +TEMPERATURE: float = 37.0 # Toy temperature for setting up metabolism [C]. diff --git a/virtual_rainforest/models/animals/functional_group.py b/virtual_rainforest/models/animals/functional_group.py index dc92adb62..1d6cd5a09 100644 --- a/virtual_rainforest/models/animals/functional_group.py +++ b/virtual_rainforest/models/animals/functional_group.py @@ -10,17 +10,7 @@ MetabolicType, TaxaType, ) -from virtual_rainforest.models.animals.constants import ( - CONVERSION_EFFICIENCY, - DAMUTHS_LAW_TERMS, - FAT_MASS_TERMS, - INTAKE_RATE_TERMS, - LONGEVITY_SCALING_TERMS, - MECHANICAL_EFFICIENCY, - METABOLIC_RATE_TERMS, - MUSCLE_MASS_TERMS, - PREY_MASS_SCALING_TERMS, -) +from virtual_rainforest.models.animals.constants import AnimalConsts class FunctionalGroup: @@ -44,6 +34,7 @@ def __init__( metabolic_type: str, birth_mass: float, adult_mass: float, + constants: AnimalConsts = AnimalConsts(), ) -> None: """The constructor for the FunctionalGroup class.""" @@ -59,27 +50,35 @@ def __init__( """The mass of the functional group at birth.""" self.adult_mass = adult_mass """The mass of the functional group at adulthood.""" - self.metabolic_rate_terms = METABOLIC_RATE_TERMS[self.metabolic_type] + self.constants = constants + """Animal constants.""" + self.metabolic_rate_terms = self.constants.metabolic_rate_terms[ + self.metabolic_type + ] """The coefficient and exponent of metabolic rate.""" - self.damuths_law_terms = DAMUTHS_LAW_TERMS[self.taxa][self.diet] + self.damuths_law_terms = self.constants.damuths_law_terms[self.taxa][self.diet] """The coefficient and exponent of damuth's law for population density.""" - self.muscle_mass_terms = MUSCLE_MASS_TERMS[self.taxa] + self.muscle_mass_terms = self.constants.muscle_mass_terms[self.taxa] """The coefficient and exponent of muscle mass allometry.""" - self.fat_mass_terms = FAT_MASS_TERMS[self.taxa] + self.fat_mass_terms = self.constants.fat_mass_terms[self.taxa] """The coefficient and exponent of fat mass allometry.""" - self.intake_rate_terms = INTAKE_RATE_TERMS[self.taxa] + self.intake_rate_terms = self.constants.intake_rate_terms[self.taxa] """The coefficient and exponent of intake allometry.""" - self.conversion_efficiency = CONVERSION_EFFICIENCY[self.diet] + self.conversion_efficiency = self.constants.conversion_efficiency[self.diet] """The conversion efficiency of the functional group based on diet.""" - self.mechanical_efficiency = MECHANICAL_EFFICIENCY[self.diet] + self.mechanical_efficiency = self.constants.mechanical_efficiency[self.diet] """The mechanical transfer efficiency of a functional group based on diet.""" - self.prey_scaling = PREY_MASS_SCALING_TERMS[self.metabolic_type][self.taxa] + self.prey_scaling = self.constants.prey_mass_scaling_terms[self.metabolic_type][ + self.taxa + ] """The predator-prey mass ratio scaling relationship.""" - self.longevity_scaling = LONGEVITY_SCALING_TERMS[self.taxa] + self.longevity_scaling = self.constants.longevity_scaling_terms[self.taxa] """The coefficient and exponent for lifespan allometry.""" -def import_functional_groups(fg_csv_file: str) -> list[FunctionalGroup]: +def import_functional_groups( + fg_csv_file: str, constants: AnimalConsts +) -> list[FunctionalGroup]: """The function to import pre-defined functional groups. This function is a first-pass of how we might import pre-defined functional groups. @@ -117,6 +116,7 @@ def import_functional_groups(fg_csv_file: str) -> list[FunctionalGroup]: row.metabolic_type, row.birth_mass, row.adult_mass, + constants=constants, ) for row in fg.itertuples() ] diff --git a/virtual_rainforest/models/animals/plant_resources.py b/virtual_rainforest/models/animals/plant_resources.py index 4d38af6d7..f6348c26b 100644 --- a/virtual_rainforest/models/animals/plant_resources.py +++ b/virtual_rainforest/models/animals/plant_resources.py @@ -5,7 +5,7 @@ from __future__ import annotations from virtual_rainforest.core.data import Data -from virtual_rainforest.models.animals.constants import ENERGY_DENSITY +from virtual_rainforest.models.animals.constants import AnimalConsts from virtual_rainforest.models.animals.protocols import Consumer, DecayPool @@ -26,7 +26,7 @@ class PlantResources: cell_id: The cell id for the plant community to expose. """ - def __init__(self, data: Data, cell_id: int) -> None: + def __init__(self, data: Data, cell_id: int, constants: AnimalConsts) -> None: # Store the data and extract the appropriate plant data self.data = data """A reference to the core data object.""" @@ -38,7 +38,7 @@ def __init__(self, data: Data, cell_id: int) -> None: # Calculate energy availability # TODO - this needs to be handed back to the plants model, which will define PFT # specific conversions to different resources. - self.energy_density: float = ENERGY_DENSITY["plant"] + self.energy_density: float = constants.energy_density["plant"] """The energy (J) in a kg of plant [currently set to toy value of Alfalfa].""" self.energy_max: float = self.mass_current * self.energy_density """The maximum amount of energy that the cohort can have [J] [Alfalfa].""" From c6943d339c1b8f32c5d4798470ece6fab8e29fdd Mon Sep 17 00:00:00 2001 From: vgro Date: Mon, 27 Nov 2023 13:19:11 +0000 Subject: [PATCH 06/18] simplified, constants description updated --- docs/source/refs.bib | 10 + tests/models/abiotic/test_abiotic_tools.py | 17 +- tests/models/abiotic/test_wind.py | 10 +- .../models/abiotic/constants.py | 178 ++++++++++++++---- virtual_rainforest/models/abiotic/wind.py | 42 +++-- 5 files changed, 198 insertions(+), 59 deletions(-) diff --git a/docs/source/refs.bib b/docs/source/refs.bib index d6a469727..01640a5b1 100644 --- a/docs/source/refs.bib +++ b/docs/source/refs.bib @@ -1,4 +1,14 @@ +@article{yasuda_turbulent_1988, + title = {Turbulent diffusivity and diurnal variations in the atmospheric boundary layer.}, + volume = {43}, + doi = {https://doi.org/10.1007/BF00128403}, + journal = {Boundary-layer meteorology}, + author = {Yasuda, N.}, + year = {1988}, + pages = {209--221}, +} + @article{henderson-sellers_new_1984, title = {A new formula for latent heat of vaporization of water as a function of temperature}, volume = {110}, diff --git a/tests/models/abiotic/test_abiotic_tools.py b/tests/models/abiotic/test_abiotic_tools.py index 418414b41..4d9fee3e6 100644 --- a/tests/models/abiotic/test_abiotic_tools.py +++ b/tests/models/abiotic/test_abiotic_tools.py @@ -23,11 +23,13 @@ def test_calculate_molar_density_air(): result, np.array( [ - [38.722469, 38.722469, 38.722469], - [39.382924, 39.382924, 39.382924], - [39.653457, 39.653457, 39.653457], + [38.749371, 38.749371, 38.749371], + [39.410285, 39.410285, 39.410285], + [39.681006, 39.681006, 39.681006], ] ), + rtol=1e-5, + atol=1e-5, ) @@ -46,7 +48,7 @@ def test_calculate_specific_heat_air(): ) exp_result = np.array( - [[29.2075, 29.207, 29.207], [29.202, 29.202, 29.202], [29.2, 29.2, 29.2]] + [[29.2075, 29.207, 29.207], [29.202, 29.202, 29.202], [29.2, 29.2, 29.2]], ) np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) @@ -73,4 +75,9 @@ def test_calculate_latent_heat_vaporisation(): ] ) - np.testing.assert_allclose(result, exp_result) + np.testing.assert_allclose( + result, + exp_result, + rtol=1e-5, + atol=1e-5, + ) diff --git a/tests/models/abiotic/test_wind.py b/tests/models/abiotic/test_wind.py index b9a7fa1ce..f1d44e55b 100644 --- a/tests/models/abiotic/test_wind.py +++ b/tests/models/abiotic/test_wind.py @@ -109,6 +109,10 @@ def test_calculate_diabatic_correction_above(dummy_data): zero_plane_displacement=np.array([0, 8, 43]), celsius_to_kelvin=AbioticConsts.celsius_to_kelvin, von_karmans_constant=AbioticConsts.von_karmans_constant, + yasuda_stability_parameter1=AbioticConsts.yasuda_stability_parameter1, + yasuda_stability_parameter2=AbioticConsts.yasuda_stability_parameter2, + yasuda_stability_parameter3=AbioticConsts.yasuda_stability_parameter3, + diabatic_heat_momentum_ratio=AbioticConsts.diabatic_heat_momentum_ratio, ) exp_result_h = np.array([0.003044, 0.026923, 0.012881]) @@ -254,9 +258,9 @@ def test_calculate_wind_canopy(dummy_data): ) exp_result = np.array( [ - [0.0, 1.797693e308, 1.797693e308], - [2.0, 7.841561, 57.16319], - [2.0, 3.188141, 8.051917], + [0.000000e000, 1.797693e308, 1.797693e308], + [2.000000e000, 6.718990e000, 4.247477e001], + [2.000000e000, 3.814989e000, 7.398743e000], ] ) np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) diff --git a/virtual_rainforest/models/abiotic/constants.py b/virtual_rainforest/models/abiotic/constants.py index 25f2d476f..7aa2faf09 100644 --- a/virtual_rainforest/models/abiotic/constants.py +++ b/virtual_rainforest/models/abiotic/constants.py @@ -6,6 +6,8 @@ from dataclasses import dataclass +from scipy.constants import Stefan_Boltzmann, atmosphere, gravitational_constant + from virtual_rainforest.core.constants_class import ConstantsDataclass @@ -14,59 +16,167 @@ class AbioticConsts(ConstantsDataclass): """Dataclass to store all constants for the `abiotic` model.""" celsius_to_kelvin: float = 273.15 - """Factor to convert temperature in Celsius to absolute temperature in Kelvin.""" - standard_pressure: float = 101.3 - """Standard atmospheric pressure, [kPa]""" - standard_mole: float = 44.6 - """Moles of ideal gas in 1 m^3 air at standard atmosphere.""" + """TODO CoreConst! Factor to convert temperature between Celsius and Kelvin.""" + + standard_pressure: float = atmosphere / 1000 + """TODO CoreConst! Standard atmospheric pressure, [kPa]""" + + standard_mole: float = 44.642 + """TODO CoreConst! Moles of ideal gas in 1 m^3 air at standard atmosphere.""" + molar_heat_capacity_air: float = 29.19 - """Molar heat capacity of air, [J mol-1 K-1].""" - gravity: float = 9.81 - """Gravitational acceleration constant, [m s-1].""" - stefan_boltzmann_constant: float = 5.67e-8 - """Stefan-Boltzmann constant, [W m-2 K-4]""" + """TODO CoreConst! Molar heat capacity of air, [J mol-1 K-1].""" + + gravity: float = gravitational_constant + """TODO CoreConst! Newtonian constant of gravitation, [m s-1].""" + + stefan_boltzmann_constant: float = Stefan_Boltzmann + """TODO CoreConst! Stefan-Boltzmann constant, [W m-2 K-4]. + + The Stefan-Boltzmann constant relates the energy radiated by a black body to its + temperature.""" + von_karmans_constant: float = 0.4 - """Von Karman's constant, unitless constant describing the logarithmic velocity - profile of a turbulent fluid near a no-slip boundary.""" + """TODO CoreConst? Von Karman's constant, [unitless]. + + The von Karman's constant describes the logarithmic velocity profile of a turbulent + fluid near a no-slip boundary.""" + specific_heat_equ_factor_1: float = 2e-05 - """Factor in calculation of molar specific heat of air.""" + """Factor in calculation of molar specific heat of air. + + Implementation after :cite:t:`maclean_microclimc_2021`.""" + specific_heat_equ_factor_2: float = 0.0002 - """Factor in calculation of molar specific heat of air.""" + """Factor in calculation of molar specific heat of air. + + Implementation after :cite:t:`maclean_microclimc_2021`.""" + latent_heat_vap_equ_factor_1: float = 1.91846e6 - """Factor in calculation of latent heat of vaporisation, (Henderson-Sellers, 1984). + """Factor in calculation of latent heat of vaporisation. + + Implementation after :cite:t:`maclean_microclimc_2021`, value is taken from + :cite:t:`henderson-sellers_new_1984`. """ latent_heat_vap_equ_factor_2: float = 33.91 - """Factor in calculation of latent heat of vaporisation, (Henderson-Sellers, 1984). + """Factor in calculation of latent heat of vaporisation. + + Implementation after :cite:t:`maclean_microclimc_2021`, value is taken from + :cite:t:`henderson-sellers_new_1984`. """ zero_plane_scaling_parameter: float = 7.5 - """Control parameter for scaling zero displacement/height - :cite:p:`raupach_simplified_1994`.""" + """Control parameter for scaling zero displacement to height, dimensionless. + + Implementation after :cite:t:`maclean_microclimc_2021`, value is taken from + :cite:p:`raupach_simplified_1994`.""" + substrate_surface_drag_coefficient: float = 0.003 - """Substrate-surface drag coefficient :cite:p:`maclean_microclimc_2021`.""" + """Substrate-surface drag coefficient, dimensionless. + + The substrate-surface drag coefficient represents the ratio of drag force exerted on + a surface to the dynamic pressure, crucial for modeling airflow over different + surface types, such as vegetation or rough terrain. Implementation and value from + :cite:t:`maclean_microclimc_2021`.""" + roughness_element_drag_coefficient: float = 0.3 - """Roughness-element drag coefficient :cite:p:`maclean_microclimc_2021`.""" + """Roughness-element drag coefficient, dimensionless. + + The roughness-element drag coefficient refers to the dimensionless coefficient used + to quantify the drag force exerted by individual roughness elements (such as + buildings, trees, or surface irregularities) on airflow, influencing the overall + aerodynamic characteristics of a surface within the atmospheric boundary layer. + Implementation and value from :cite:t:`maclean_microclimc_2021`.""" + roughness_sublayer_depth_parameter: float = 0.193 - """Parameter characterizes the roughness sublayer depth - :cite:p:`maclean_microclimc_2021`.""" + """Parameter characterizes the roughness sublayer depth. + + The roughness sublayer depth refers to the layer near the surface where the effects + of surface roughness significantly influence airflow, turbulence, and momentum + transfer, typically extending up to about 10% of the height of the roughness + elements or features on the surface. This layer is characterized by intense + turbulence and rapid velocity changes due to surface irregularities. + Implentation and value taken from :cite:p:`maclean_microclimc_2021`.""" + max_ratio_wind_to_friction_velocity: float = 0.3 - """Maximum ratio of wind velocity to friction velocity - :cite:p:`maclean_microclimc_2021`.""" + """Maximum ratio of wind velocity to friction velocity, dimensionless. + + Implementation and value from :cite:t:`maclean_microclimc_2021`.""" + drag_coefficient: float = 0.2 - """Drag coefficient :cite:p:`maclean_microclimc_2021`.""" + """Drag coefficient, dimensionless. + + The drag coefficient is a dimensionless quantity that characterizes the drag or + resistance experienced by an object moving through a fluid (here the atmosphere) and + is defined as the ratio of the drag force on the object to the dynamic pressure of + the fluid flow and the reference area of the object. + Implementation and value from :cite:t:`maclean_microclimc_2021`.""" + relative_turbulence_intensity: float = 0.5 - """Relative turbulence intensity :cite:p:`maclean_microclimc_2021`.""" + """Relative turbulence intensity. + + The relative turbulence intensity is a proportionality factor that relates the mean + eddy velocity is assumed to the local wind speed below the canopy. Implementation + and value from :cite:t:`maclean_microclimc_2021`.""" + diabatic_correction_factor_below: float = 1 - "Diabatic correction factor below canopy." + """Diabatic correction factor below canopy. + + The diabatic correction factor is a scaling adjustment used to compensate for the + effects of vertical heat transfer or thermal non-adiabaticity on atmospheric + variables or processes, particularly when estimating or interpreting measurements + across different heights or conditions. This factor is used to adjust wind profiles + below the canopy. Implementation and value from :cite:t:`maclean_microclimc_2021`. + """ + mixing_length_factor: float = 0.32 - """Factor in calculation of mixing length :cite:p:`maclean_microclimc_2021`.""" + """Factor in calculation of mixing length. + + Implementation and value from :cite:t:`maclean_microclimc_2021`.""" + min_relative_turbulence_intensity: float = 0.36 - """minimum relative turbulence intensity, default value from Shaw et al (1974) - Agricultural Meteorology, 13: 419-425. """ + """Minimum relative turbulence intensity. + + See :attr:`relative_turbulence_intensity`. + The default value is taken from Shaw et al (1974) Agricultural Meteorology, 13: + 419-425. TODO this is not representative of a rainforest environment and needs to be + adjusted. + """ + max_relative_turbulence_intensity: float = 0.9 - """maximum relative turbulence intensity, default value from Shaw et al (1974) - Agricultural Meteorology, 13: 419-425.""" + """Maximum relative turbulence intensity. + + See :attr:`relative_turbulence_intensity`. + The default value from Shaw et al (1974) Agricultural Meteorology, 13: 419-425. + TODO this is not representative of a rainforest environment and needs to be + adjusted.""" + min_wind_speed_above_canopy: float = 0.55 - """Minimum wind speed above the canopy, [m/s].""" + """Minimum wind speed above the canopy, [m/s]. + + Implementation and value from :cite:t:`maclean_microclimc_2021`.""" + min_roughness_length: float = 0.05 - """Minimum roughness length, [m].""" + """Minimum roughness length, [m]. + + Implementation and value from :cite:t:`maclean_microclimc_2021`.""" + + yasuda_stability_parameter1: float = 6 + """Parameter to approximate diabatic correction factors for heat and momentum. + + after :cite:t:`yasuda_turbulent_1988`.""" + + yasuda_stability_parameter2: float = 2 + """Parameter to approximate diabatic correction factors for heat and momentum. + + Value is taken from :cite:t:`yasuda_turbulent_1988.""" + + yasuda_stability_parameter3: float = 16 + """Parameter to approximate diabatic correction factors for heat and momentum. + + Value is taken :cite:t:`yasuda_turbulent_1988`.""" + + diabatic_heat_momentum_ratio: float = 0.6 + """Factor that relates diabatic correction factors for heat and momentum. + + Value is taken :cite:t:`yasuda_turbulent_1988`.""" diff --git a/virtual_rainforest/models/abiotic/wind.py b/virtual_rainforest/models/abiotic/wind.py index 875e11fea..12fde81b0 100644 --- a/virtual_rainforest/models/abiotic/wind.py +++ b/virtual_rainforest/models/abiotic/wind.py @@ -59,20 +59,16 @@ def calculate_zero_plane_displacement( zero place displacement height, [m] """ - # Calculate in presence of vegetation - leaf_area_index_displacement = np.where( - leaf_area_index > 0, leaf_area_index, np.nan - ) - - scale_displacement = np.sqrt( - zero_plane_scaling_parameter * leaf_area_index_displacement - ) + # Select grid cells where vegetation is present + displacement = np.where(leaf_area_index > 0, leaf_area_index, np.nan) + # Calculate zero displacement height + scale_displacement = np.sqrt(zero_plane_scaling_parameter * displacement) zero_place_displacement = ( (1 - (1 - np.exp(-scale_displacement)) / scale_displacement) * canopy_height, ) - # no displacement in absence of vegetation + # No displacement in absence of vegetation return np.nan_to_num(zero_place_displacement, nan=0.0).squeeze() @@ -154,6 +150,10 @@ def calculate_diabatic_correction_above( zero_plane_displacement: NDArray[np.float32], celsius_to_kelvin: float, von_karmans_constant: float, + yasuda_stability_parameter1: float, + yasuda_stability_parameter2: float, + yasuda_stability_parameter3: float, + diabatic_heat_momentum_ratio: float, ) -> dict[str, NDArray[np.float32]]: r"""Calculates the diabatic correction factors for momentum and heat above canopy. @@ -176,6 +176,14 @@ def calculate_diabatic_correction_above( celsius_to_kelvin: Factor to convert temperature in Celsius to absolute temperature in Kelvin von_karmans_constant + yasuda_stability_parameter1: Parameter to approximate diabatic correction + factors for heat and momentum after :cite:t:`yasuda_turbulent_1988` + yasuda_stability_parameter2: Parameter to approximate diabatic correction + factors for heat and momentum after :cite:t:`yasuda_turbulent_1988` + yasuda_stability_parameter3: Parameter to approximate diabatic correction + factors for heat and momentum after :cite:t:`yasuda_turbulent_1988` + diabatic_heat_momentum_ratio: Factor that relates diabatic correction + factors for heat and momentum after :cite:t:`yasuda_turbulent_1988` Returns: diabatic correction factors for heat (psi_h) and momentum (psi_m) transfer @@ -193,21 +201,21 @@ def calculate_diabatic_correction_above( * friction_velocity ) - stable_condition = 6 * np.log(1 - stability) - unstable_condition = -2 * np.log((1 + np.sqrt(1 - 16 * stability)) / 2) + stable_condition = yasuda_stability_parameter1 * np.log(1 - stability) + unstable_condition = -yasuda_stability_parameter2 * np.log( + (1 + np.sqrt(1 - yasuda_stability_parameter3 * stability)) / 2 + ) diabatic_correction_heat = np.where( - stability < 0, stable_condition, unstable_condition + sensible_heat_flux < 0, stable_condition, unstable_condition ) diabatic_correction_momentum = np.where( - stability < 0, diabatic_correction_heat, 0.6 * diabatic_correction_heat + sensible_heat_flux < 0, + diabatic_correction_heat, + diabatic_heat_momentum_ratio * diabatic_correction_heat, ) - # Apply upper threshold - diabatic_correction_momentum = np.minimum(diabatic_correction_momentum, 5) - diabatic_correction_heat = np.minimum(diabatic_correction_heat, 5) - return {"psi_m": diabatic_correction_momentum, "psi_h": diabatic_correction_heat} From 6ccb75077d9935cc09693513909e08b4751d3837 Mon Sep 17 00:00:00 2001 From: vgro Date: Mon, 27 Nov 2023 13:37:26 +0000 Subject: [PATCH 07/18] constants moved to CoreConsts --- tests/models/abiotic/test_abiotic_tools.py | 11 +++---- tests/models/abiotic/test_wind.py | 11 +++---- virtual_rainforest/core/constants.py | 28 ++++++++++++++++-- .../models/abiotic/constants.py | 29 ------------------- 4 files changed, 38 insertions(+), 41 deletions(-) diff --git a/tests/models/abiotic/test_abiotic_tools.py b/tests/models/abiotic/test_abiotic_tools.py index 4d9fee3e6..f45d03bab 100644 --- a/tests/models/abiotic/test_abiotic_tools.py +++ b/tests/models/abiotic/test_abiotic_tools.py @@ -2,6 +2,7 @@ import numpy as np +from virtual_rainforest.core.constants import CoreConsts from virtual_rainforest.models.abiotic.constants import AbioticConsts @@ -15,9 +16,9 @@ def test_calculate_molar_density_air(): result = calculate_molar_density_air( temperature=np.array([[25, 25, 25], [20, 20, 20], [18, 18, 18]]), atmospheric_pressure=np.full((3, 3), 96), - standard_mole=AbioticConsts.standard_mole, - standard_pressure=AbioticConsts.standard_pressure, - celsius_to_kelvin=AbioticConsts.celsius_to_kelvin, + standard_mole=CoreConsts.standard_mole, + standard_pressure=CoreConsts.standard_pressure, + celsius_to_kelvin=CoreConsts.zero_Celsius, ) np.testing.assert_allclose( result, @@ -42,7 +43,7 @@ def test_calculate_specific_heat_air(): result = calculate_specific_heat_air( temperature=np.array([[25, 25, 25], [20, 20, 20], [18, 18, 18]]), - molar_heat_capacity_air=AbioticConsts.molar_heat_capacity_air, + molar_heat_capacity_air=CoreConsts.molar_heat_capacity_air, specific_heat_equ_factor_1=AbioticConsts.specific_heat_equ_factor_1, specific_heat_equ_factor_2=AbioticConsts.specific_heat_equ_factor_2, ) @@ -63,7 +64,7 @@ def test_calculate_latent_heat_vaporisation(): result = calculate_latent_heat_vaporisation( temperature=np.array([[25, 25, 25], [20, 20, 20], [18, 18, 18]]), - celsius_to_kelvin=AbioticConsts.celsius_to_kelvin, + celsius_to_kelvin=CoreConsts.zero_Celsius, latent_heat_vap_equ_factor_1=AbioticConsts.latent_heat_vap_equ_factor_1, latent_heat_vap_equ_factor_2=AbioticConsts.latent_heat_vap_equ_factor_2, ) diff --git a/tests/models/abiotic/test_wind.py b/tests/models/abiotic/test_wind.py index f1d44e55b..79c3b5a0c 100644 --- a/tests/models/abiotic/test_wind.py +++ b/tests/models/abiotic/test_wind.py @@ -5,6 +5,7 @@ import pytest from xarray import DataArray +from virtual_rainforest.core.constants import CoreConsts from virtual_rainforest.models.abiotic.constants import AbioticConsts @@ -84,7 +85,7 @@ def test_calculate_roughness_length_momentum(dummy_data): AbioticConsts.max_ratio_wind_to_friction_velocity ), min_roughness_length=AbioticConsts.min_roughness_length, - von_karman_constant=AbioticConsts.von_karmans_constant, + von_karman_constant=CoreConsts.von_karmans_constant, ) np.testing.assert_allclose( @@ -107,8 +108,8 @@ def test_calculate_diabatic_correction_above(dummy_data): friction_velocity=dummy_data["friction_velocity"].to_numpy(), wind_heights=np.array([1, 15, 50]), zero_plane_displacement=np.array([0, 8, 43]), - celsius_to_kelvin=AbioticConsts.celsius_to_kelvin, - von_karmans_constant=AbioticConsts.von_karmans_constant, + celsius_to_kelvin=CoreConsts.zero_Celsius, + von_karmans_constant=CoreConsts.von_karmans_constant, yasuda_stability_parameter1=AbioticConsts.yasuda_stability_parameter1, yasuda_stability_parameter2=AbioticConsts.yasuda_stability_parameter2, yasuda_stability_parameter3=AbioticConsts.yasuda_stability_parameter3, @@ -217,7 +218,7 @@ def test_calculate_friction_velocity(dummy_data): zeroplane_displacement=np.array([0, 8, 43]), roughness_length_momentum=np.array([0.003, 0.435, 1.521]), diabatic_correction_momentum=np.array([-0.1, 0.0, 0.1]), - von_karmans_constant=AbioticConsts.von_karmans_constant, + von_karmans_constant=CoreConsts.von_karmans_constant, ) exp_result = np.array([0.0, 1.310997, 4.0]) np.testing.assert_allclose(result, exp_result) @@ -234,7 +235,7 @@ def test_calculate_wind_above_canopy(): zeroplane_displacement=np.array([0, 8, 43]), roughness_length_momentum=np.array([0.003, 0.435, 1.521]), diabatic_correction_momentum=np.array([-0.1, 0.0, 0.1]), - von_karmans_constant=AbioticConsts.von_karmans_constant, + von_karmans_constant=CoreConsts.von_karmans_constant, min_wind_speed_above_canopy=AbioticConsts.min_wind_speed_above_canopy, ) diff --git a/virtual_rainforest/core/constants.py b/virtual_rainforest/core/constants.py index eb8b550c9..ce992056c 100644 --- a/virtual_rainforest/core/constants.py +++ b/virtual_rainforest/core/constants.py @@ -18,8 +18,32 @@ class CoreConsts(ConstantsDataclass): """Core constants for use across the Virtual Rainforest modules.""" + placeholder: float = 123.4 + """A placeholder configurable constant.""" + zero_Celsius: ClassVar[float] = constants.zero_Celsius """Conversion constant from Kelvin to Celsius (°).""" - placeholder: float = 123.4 - """A placeholder configurable constant.""" + standard_pressure: float = constants.atmosphere / 1000 + """Standard atmospheric pressure, [kPa]""" + + standard_mole: float = 44.642 + """Moles of ideal gas in 1 m^3 air at standard atmosphere.""" + + molar_heat_capacity_air: float = 29.19 + """Molar heat capacity of air, [J mol-1 K-1].""" + + gravity: float = constants.gravitational_constant + """Newtonian constant of gravitation, [m s-1].""" + + stefan_boltzmann_constant: float = constants.Stefan_Boltzmann + """Stefan-Boltzmann constant, [W m-2 K-4]. + + The Stefan-Boltzmann constant relates the energy radiated by a black body to its + temperature.""" + + von_karmans_constant: float = 0.4 + """Von Karman's constant, [unitless]. + + The von Karman's constant describes the logarithmic velocity profile of a turbulent + fluid near a no-slip boundary.""" diff --git a/virtual_rainforest/models/abiotic/constants.py b/virtual_rainforest/models/abiotic/constants.py index 7aa2faf09..6609c37ac 100644 --- a/virtual_rainforest/models/abiotic/constants.py +++ b/virtual_rainforest/models/abiotic/constants.py @@ -6,8 +6,6 @@ from dataclasses import dataclass -from scipy.constants import Stefan_Boltzmann, atmosphere, gravitational_constant - from virtual_rainforest.core.constants_class import ConstantsDataclass @@ -15,33 +13,6 @@ class AbioticConsts(ConstantsDataclass): """Dataclass to store all constants for the `abiotic` model.""" - celsius_to_kelvin: float = 273.15 - """TODO CoreConst! Factor to convert temperature between Celsius and Kelvin.""" - - standard_pressure: float = atmosphere / 1000 - """TODO CoreConst! Standard atmospheric pressure, [kPa]""" - - standard_mole: float = 44.642 - """TODO CoreConst! Moles of ideal gas in 1 m^3 air at standard atmosphere.""" - - molar_heat_capacity_air: float = 29.19 - """TODO CoreConst! Molar heat capacity of air, [J mol-1 K-1].""" - - gravity: float = gravitational_constant - """TODO CoreConst! Newtonian constant of gravitation, [m s-1].""" - - stefan_boltzmann_constant: float = Stefan_Boltzmann - """TODO CoreConst! Stefan-Boltzmann constant, [W m-2 K-4]. - - The Stefan-Boltzmann constant relates the energy radiated by a black body to its - temperature.""" - - von_karmans_constant: float = 0.4 - """TODO CoreConst? Von Karman's constant, [unitless]. - - The von Karman's constant describes the logarithmic velocity profile of a turbulent - fluid near a no-slip boundary.""" - specific_heat_equ_factor_1: float = 2e-05 """Factor in calculation of molar specific heat of air. From fef089fc359d04a79faca24a53012cc6085cbddc Mon Sep 17 00:00:00 2001 From: vgro Date: Mon, 27 Nov 2023 14:37:16 +0000 Subject: [PATCH 08/18] bugfix documentation --- virtual_rainforest/models/abiotic/constants.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/virtual_rainforest/models/abiotic/constants.py b/virtual_rainforest/models/abiotic/constants.py index 6609c37ac..064068542 100644 --- a/virtual_rainforest/models/abiotic/constants.py +++ b/virtual_rainforest/models/abiotic/constants.py @@ -40,7 +40,7 @@ class AbioticConsts(ConstantsDataclass): """Control parameter for scaling zero displacement to height, dimensionless. Implementation after :cite:t:`maclean_microclimc_2021`, value is taken from - :cite:p:`raupach_simplified_1994`.""" + :cite:t:`raupach_simplified_1994`.""" substrate_surface_drag_coefficient: float = 0.003 """Substrate-surface drag coefficient, dimensionless. @@ -140,7 +140,7 @@ class AbioticConsts(ConstantsDataclass): yasuda_stability_parameter2: float = 2 """Parameter to approximate diabatic correction factors for heat and momentum. - Value is taken from :cite:t:`yasuda_turbulent_1988.""" + Value is taken from :cite:t:`yasuda_turbulent_1988`.""" yasuda_stability_parameter3: float = 16 """Parameter to approximate diabatic correction factors for heat and momentum. From 95e9f60535028389282137aec94bf90c866bf88f Mon Sep 17 00:00:00 2001 From: vgro Date: Tue, 28 Nov 2023 11:48:22 +0000 Subject: [PATCH 09/18] merged dummy data sets, fix dimensions attencoeff --- tests/conftest.py | 5 +- tests/models/abiotic/test_wind.py | 163 ++++++++++-------- .../models/abiotic/constants.py | 35 ++-- virtual_rainforest/models/abiotic/wind.py | 8 +- 4 files changed, 122 insertions(+), 89 deletions(-) diff --git a/tests/conftest.py b/tests/conftest.py index f2b432ccf..ef89a6311 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -472,5 +472,8 @@ def dummy_climate_data(layer_roles_fixture): np.full((2, 3), 450), dims=("groundwater_layers", "cell_id"), ) - + data["canopy_height"] = DataArray([0, 10, 50], dims=["cell_id"]) + data["sensible_heat_flux_toc"] = DataArray([100, 50, 10], dims=["cell_id"]) + data["friction_velocity"] = DataArray([12, 5, 2], dims=["cell_id"]) + data["wind_speed_ref"] = DataArray([0, 5, 10], dims=["cell_id"]) return data diff --git a/tests/models/abiotic/test_wind.py b/tests/models/abiotic/test_wind.py index 79c3b5a0c..cc5607009 100644 --- a/tests/models/abiotic/test_wind.py +++ b/tests/models/abiotic/test_wind.py @@ -2,74 +2,36 @@ import numpy as np -import pytest -from xarray import DataArray from virtual_rainforest.core.constants import CoreConsts from virtual_rainforest.models.abiotic.constants import AbioticConsts -@pytest.fixture -def dummy_data(): - """Creates a dummy data object for use in wind tests. - - One grid cell has no vegetation, two grid cells represent a range of values. - """ - - from virtual_rainforest.core.data import Data - from virtual_rainforest.core.grid import Grid - - # Setup the data object with two cells. - grid = Grid(cell_nx=3, cell_ny=1) - data = Data(grid) - - # Add the required data. - data["canopy_height"] = DataArray([0, 10, 50], dims=["cell_id"]) - data["wind_layer_heights"] = DataArray( - [[0.1, 5, 15], [1, 15, 25], [5, 25, 55]], dims=["cell_id", "layers"] - ) - data["leaf_area_index"] = DataArray( - [ - [0, 1, 5], - [0, 2, 5], - [0, 3, 5], - ], - dims=["cell_id", "layers"], - ) - data["air_temperature"] = DataArray([30, 20, 30], dims=["cell_id"]) - data["atmospheric_pressure"] = DataArray([101, 102, 103], dims=["cell_id"]) - data["sensible_heat_flux_topofcanopy"] = DataArray([100, 50, 10], dims=["cell_id"]) - data["friction_velocity"] = DataArray([12, 5, 2], dims=["cell_id"]) - data["wind_speed_ref"] = DataArray([0, 5, 10], dims=["cell_id"]) - - return data - - -def test_calculate_zero_plane_displacement(dummy_data): +def test_calculate_zero_plane_displacement(dummy_climate_data): """Test if calculated correctly and set to zero without vegetation.""" from virtual_rainforest.models.abiotic.wind import calculate_zero_plane_displacement - leaf_area_index = dummy_data.data["leaf_area_index"].sum(dim="layers") + leaf_area_index = dummy_climate_data.data["leaf_area_index"].sum(dim="layers") result = calculate_zero_plane_displacement( - canopy_height=dummy_data.data["canopy_height"].to_numpy(), + canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(), leaf_area_index=leaf_area_index.to_numpy(), zero_plane_scaling_parameter=AbioticConsts.zero_plane_scaling_parameter, ) - np.testing.assert_allclose(result, np.array([0.0, 8.620853, 43.547819])) + np.testing.assert_allclose(result, np.array([0.0, 7.910175, 39.550874])) -def test_calculate_roughness_length_momentum(dummy_data): +def test_calculate_roughness_length_momentum(dummy_climate_data): """Test roughness length governing momentum transfer.""" from virtual_rainforest.models.abiotic.wind import ( calculate_roughness_length_momentum, ) - leaf_area_index = dummy_data.data["leaf_area_index"].sum(dim="layers") + leaf_area_index = dummy_climate_data.data["leaf_area_index"].sum(dim="layers") result = calculate_roughness_length_momentum( - canopy_height=dummy_data.data["canopy_height"].to_numpy(), + canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(), leaf_area_index=leaf_area_index.to_numpy(), zero_plane_displacement=np.array([0, 8, 43]), substrate_surface_drag_coefficient=( @@ -93,19 +55,24 @@ def test_calculate_roughness_length_momentum(dummy_data): ) -def test_calculate_diabatic_correction_above(dummy_data): +def test_calculate_diabatic_correction_above(dummy_climate_data): """Test diabatic correction factors for heat and momentum.""" from virtual_rainforest.models.abiotic.wind import ( calculate_diabatic_correction_above, ) + air_temperature = ( + dummy_climate_data.data["air_temperature"] + .where(dummy_climate_data.data["air_temperature"].layer_roles != "soil") + .dropna(dim="layers") + ) result = calculate_diabatic_correction_above( molar_density_air=np.repeat(28.96, 3), specific_heat_air=np.repeat(1, 3), - temperature=dummy_data["air_temperature"].to_numpy(), - sensible_heat_flux=dummy_data["sensible_heat_flux_topofcanopy"].to_numpy(), - friction_velocity=dummy_data["friction_velocity"].to_numpy(), + temperature=air_temperature.to_numpy(), + sensible_heat_flux=dummy_climate_data.data["sensible_heat_flux_toc"].to_numpy(), + friction_velocity=dummy_climate_data.data["friction_velocity"].to_numpy(), wind_heights=np.array([1, 15, 50]), zero_plane_displacement=np.array([0, 8, 43]), celsius_to_kelvin=CoreConsts.zero_Celsius, @@ -116,19 +83,37 @@ def test_calculate_diabatic_correction_above(dummy_data): diabatic_heat_momentum_ratio=AbioticConsts.diabatic_heat_momentum_ratio, ) - exp_result_h = np.array([0.003044, 0.026923, 0.012881]) - exp_result_m = np.array([0.001827, 0.016154, 0.007729]) + exp_result_h = np.array( + [ + [0.003044, 0.026017, 0.012881], + [0.003046, 0.026031, 0.012888], + [0.003056, 0.026117, 0.01293], + [0.003073, 0.026264, 0.013002], + [0.00312, 0.026677, 0.013204], + [0.003191, 0.027289, 0.013504], + ] + ) + exp_result_m = np.array( + [ + [0.001827, 0.01561, 0.007729], + [0.001828, 0.015619, 0.007733], + [0.001833, 0.01567, 0.007758], + [0.001833, 0.01567, 0.007758], + [0.001872, 0.016006, 0.007923], + [0.001914, 0.016374, 0.008103], + ] + ) np.testing.assert_allclose(result["psi_h"], exp_result_h, rtol=1e-4, atol=1e-4) np.testing.assert_allclose(result["psi_m"], exp_result_m, rtol=1e-4, atol=1e-4) -def test_calculate_mean_mixing_length(dummy_data): +def test_calculate_mean_mixing_length(dummy_climate_data): """Test mixing length with and without vegetation.""" from virtual_rainforest.models.abiotic.wind import calculate_mean_mixing_length result = calculate_mean_mixing_length( - canopy_height=dummy_data.data["canopy_height"].to_numpy(), + canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(), zero_plane_displacement=np.array([0, 8, 43]), roughness_length_momentum=np.array([0.003, 0.435, 1.521]), mixing_length_factor=AbioticConsts.mixing_length_factor, @@ -139,15 +124,20 @@ def test_calculate_mean_mixing_length(dummy_data): ) -def test_generate_relative_turbulence_intensity(dummy_data): +def test_generate_relative_turbulence_intensity(dummy_climate_data): """Test relative turbulence intensity.""" from virtual_rainforest.models.abiotic.wind import ( generate_relative_turbulence_intensity, ) + layer_heights = ( + dummy_climate_data.data["layer_heights"] + .where(dummy_climate_data.data["layer_heights"].layer_roles != "soil") + .dropna(dim="layers") + ) result = generate_relative_turbulence_intensity( - layer_heights=dummy_data["wind_layer_heights"].to_numpy(), + layer_heights=layer_heights.to_numpy(), min_relative_turbulence_intensity=( AbioticConsts.min_relative_turbulence_intensity ), @@ -158,24 +148,39 @@ def test_generate_relative_turbulence_intensity(dummy_data): ) exp_result = np.array( - [[0.414, 3.06, 8.46], [0.9, 8.46, 13.86], [3.06, 13.86, 30.06]] + [ + [17.64, 17.64, 17.64], + [16.56, 16.56, 16.56], + [11.16, 11.16, 11.166], + [5.76, 5.76, 5.76], + [1.17, 1.17, 1.17], + [0.414, 0.414, 0.414], + ] ) - np.testing.assert_allclose(result, exp_result) + np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) -def test_calculate_wind_attenuation_coefficient(dummy_data): +def test_calculate_wind_attenuation_coefficient(dummy_climate_data): """Test wind attenuation coefficient with and without vegetation.""" from virtual_rainforest.models.abiotic.wind import ( calculate_wind_attenuation_coefficient, ) - leaf_area_index = dummy_data.data["leaf_area_index"].to_numpy() + leaf_area_index = ( + dummy_climate_data.data["leaf_area_index"] + .where(dummy_climate_data.data["leaf_area_index"].layer_roles != "soil") + .dropna(dim="layers") + ) relative_turbulence_intensity = np.array( - [[0.414, 3.06, 8.46], [0.9, 8.46, 13.86], [3.06, 13.86, 30.06]] + [ + [17.64, 17.64, 17.64], + [16.56, 16.56, 16.56], + [11.16, 11.16, 11.166], + ] ) result = calculate_wind_attenuation_coefficient( - canopy_height=dummy_data.data["canopy_height"].to_numpy(), + canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(), leaf_area_index=leaf_area_index, mean_mixing_length=np.array([0, 0.4, 1.5]), drag_coefficient=AbioticConsts.drag_coefficient, @@ -183,7 +188,11 @@ def test_calculate_wind_attenuation_coefficient(dummy_data): ) exp_result = np.array( - [[0.0, 0.816, 1.970], [0.0, 0.591, 1.202], [0.0, 0.541, 0.554]] + [ + [0.0, 0.141723, 0.188964], + [0.0, 0.150966, 0.201288], + [0.0, 0.224014, 0.298525], + ] ) np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) @@ -207,13 +216,13 @@ def test_wind_log_profile(): ) -def test_calculate_friction_velocity(dummy_data): +def test_calculate_friction_velocity(dummy_climate_data): """Calculate friction velocity.""" from virtual_rainforest.models.abiotic.wind import calculate_fricition_velocity result = calculate_fricition_velocity( - wind_speed_ref=dummy_data["wind_speed_ref"], + wind_speed_ref=dummy_climate_data.data["wind_speed_ref"], reference_height=10.0, zeroplane_displacement=np.array([0, 8, 43]), roughness_length_momentum=np.array([0.003, 0.435, 1.521]), @@ -243,25 +252,33 @@ def test_calculate_wind_above_canopy(): np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) -def test_calculate_wind_canopy(dummy_data): +def test_calculate_wind_canopy(dummy_climate_data): """Test below canopy wind profile.""" from virtual_rainforest.models.abiotic.wind import calculate_wind_canopy - atten_coeff = np.array( - [[0.0, 0.816, 1.970], [0.0, 0.591, 1.202], [0.0, 0.541, 0.554]] + layer_heights = ( + dummy_climate_data.data["layer_heights"] + .where(dummy_climate_data.data["layer_heights"].layer_roles != "soil") + .dropna(dim="layers") ) + result = calculate_wind_canopy( top_of_canopy_wind_speed=np.array([2, 5, 7]), - wind_layer_heights=dummy_data["wind_layer_heights"], - canopy_height=dummy_data["canopy_height"], - attenuation_coefficient=atten_coeff, + wind_layer_heights=layer_heights.to_numpy(), + canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(), + attenuation_coefficient=np.array([0.0, 0.64, 1.24]), + min_windspeed_below_canopy=0.001, ) + exp_result = np.array( [ - [0.000000e000, 1.797693e308, 1.797693e308], - [2.000000e000, 6.718990e000, 4.247477e001], - [2.000000e000, 3.814989e000, 7.398743e000], + [0.001, 20.438858, 4.479494], + [0.001, 17.983199, 4.262731], + [0.001, 9.482404, 3.326465], + [0.001, 5.0, 2.59584], + [0.001, 2.90211, 2.102464], + [0.001, 2.65339, 2.030719], ] ) np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) diff --git a/virtual_rainforest/models/abiotic/constants.py b/virtual_rainforest/models/abiotic/constants.py index 064068542..9c15529d0 100644 --- a/virtual_rainforest/models/abiotic/constants.py +++ b/virtual_rainforest/models/abiotic/constants.py @@ -45,10 +45,10 @@ class AbioticConsts(ConstantsDataclass): substrate_surface_drag_coefficient: float = 0.003 """Substrate-surface drag coefficient, dimensionless. - The substrate-surface drag coefficient represents the ratio of drag force exerted on - a surface to the dynamic pressure, crucial for modeling airflow over different - surface types, such as vegetation or rough terrain. Implementation and value from - :cite:t:`maclean_microclimc_2021`.""" + The substrate-surface drag coefficient represents the resistance encountered by an + object moving on or through a surface and varies based on the nature of the surface + and the object's properties. Here, it affects how wind speed is altered by a surface + . Implementation and value from :cite:t:`maclean_microclimc_2021`.""" roughness_element_drag_coefficient: float = 0.3 """Roughness-element drag coefficient, dimensionless. @@ -84,14 +84,14 @@ class AbioticConsts(ConstantsDataclass): Implementation and value from :cite:t:`maclean_microclimc_2021`.""" relative_turbulence_intensity: float = 0.5 - """Relative turbulence intensity. + """Relative turbulence intensity, dimensionless. The relative turbulence intensity is a proportionality factor that relates the mean eddy velocity is assumed to the local wind speed below the canopy. Implementation and value from :cite:t:`maclean_microclimc_2021`.""" diabatic_correction_factor_below: float = 1 - """Diabatic correction factor below canopy. + """Diabatic correction factor below canopy, dimensionless. The diabatic correction factor is a scaling adjustment used to compensate for the effects of vertical heat transfer or thermal non-adiabaticity on atmospheric @@ -101,12 +101,12 @@ class AbioticConsts(ConstantsDataclass): """ mixing_length_factor: float = 0.32 - """Factor in calculation of mixing length. + """Factor in calculation of mixing length, dimensionless. Implementation and value from :cite:t:`maclean_microclimc_2021`.""" min_relative_turbulence_intensity: float = 0.36 - """Minimum relative turbulence intensity. + """Minimum relative turbulence intensity, dimensionless. See :attr:`relative_turbulence_intensity`. The default value is taken from Shaw et al (1974) Agricultural Meteorology, 13: @@ -115,7 +115,7 @@ class AbioticConsts(ConstantsDataclass): """ max_relative_turbulence_intensity: float = 0.9 - """Maximum relative turbulence intensity. + """Maximum relative turbulence intensity, dimensionless. See :attr:`relative_turbulence_intensity`. The default value from Shaw et al (1974) Agricultural Meteorology, 13: 419-425. @@ -127,27 +127,36 @@ class AbioticConsts(ConstantsDataclass): Implementation and value from :cite:t:`maclean_microclimc_2021`.""" + min_windspeed_below_canopy: float = 0.001 + """Minimum wind speed below the canopy or in absence of vegetation, [m/s].""" + min_roughness_length: float = 0.05 """Minimum roughness length, [m]. + The minimum roughness length represents the lowest height at which the surface + roughness significantly affects the wind flow over a particular terrain or surface. Implementation and value from :cite:t:`maclean_microclimc_2021`.""" yasuda_stability_parameter1: float = 6 """Parameter to approximate diabatic correction factors for heat and momentum. - after :cite:t:`yasuda_turbulent_1988`.""" + Dimenionless parameter, implementation after :cite:t:`maclean_microclimc_2021` and + values taken from :cite:t:`yasuda_turbulent_1988`.""" yasuda_stability_parameter2: float = 2 """Parameter to approximate diabatic correction factors for heat and momentum. - Value is taken from :cite:t:`yasuda_turbulent_1988`.""" + Dimenionless parameter, implementation after :cite:t:`maclean_microclimc_2021` and + values taken from :cite:t:`yasuda_turbulent_1988`.""" yasuda_stability_parameter3: float = 16 """Parameter to approximate diabatic correction factors for heat and momentum. - Value is taken :cite:t:`yasuda_turbulent_1988`.""" + Dimenionless parameter, implementation after :cite:t:`maclean_microclimc_2021` and + values taken from :cite:t:`yasuda_turbulent_1988`.""" diabatic_heat_momentum_ratio: float = 0.6 """Factor that relates diabatic correction factors for heat and momentum. - Value is taken :cite:t:`yasuda_turbulent_1988`.""" + Dimenionless parameter, implementation after :cite:t:`maclean_microclimc_2021` and + values taken from :cite:t:`yasuda_turbulent_1988`.""" diff --git a/virtual_rainforest/models/abiotic/wind.py b/virtual_rainforest/models/abiotic/wind.py index 12fde81b0..83ddf5d48 100644 --- a/virtual_rainforest/models/abiotic/wind.py +++ b/virtual_rainforest/models/abiotic/wind.py @@ -439,6 +439,7 @@ def calculate_wind_canopy( wind_layer_heights: NDArray[np.float32], canopy_height: NDArray[np.float32], attenuation_coefficient: NDArray[np.float32], + min_windspeed_below_canopy: float, ) -> NDArray[np.float32]: """Calculate wind profile for individual canopy layers with variable leaf area. @@ -448,8 +449,11 @@ def calculate_wind_canopy( top_of_canopy_wind_speed: Wind speed at top of canopy layer, [m /s] wind_layer_heights: Height of canopy layer node, [m] canopy_height: Height to top of canopy layer, [m] - attenuation_coefficient: Attenuation coefficient as returned by + attenuation_coefficient: Mean attenuation coefficient based on the profile + calculated by :func:`~virtual_rainforest.models.abiotic.wind.calculate_wind_attenuation_coefficient` + min_windspeed_below_canopy: Minimum wind speed below the canopy or in absence of + vegetation, [m/s]. This value is set to avoid dividion by zero. Returns: wind speed at height of canopy node, [m/s] @@ -458,5 +462,5 @@ def calculate_wind_canopy( return np.nan_to_num( top_of_canopy_wind_speed * np.exp(attenuation_coefficient * ((wind_layer_heights / canopy_height) - 1)), - nan=0, + nan=min_windspeed_below_canopy, ).squeeze() From 8b788746640e13232aca96a2e909317728234f09 Mon Sep 17 00:00:00 2001 From: Taran Rallings Date: Tue, 28 Nov 2023 12:07:46 +0000 Subject: [PATCH 10/18] Fixed a small bug in plant_resource testing. --- tests/models/animals/conftest.py | 6 ++++-- virtual_rainforest/models/animals/plant_resources.py | 5 +++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/tests/models/animals/conftest.py b/tests/models/animals/conftest.py index 0c27262be..0b2039bf4 100644 --- a/tests/models/animals/conftest.py +++ b/tests/models/animals/conftest.py @@ -166,8 +166,10 @@ def excrement_instance(): @pytest.fixture -def plant_instance(plant_data_instance): +def plant_instance(plant_data_instance, constants_instance): """Fixture for a plant community used in tests.""" from virtual_rainforest.models.animals.plant_resources import PlantResources - return PlantResources(data=plant_data_instance, cell_id=4) + return PlantResources( + data=plant_data_instance, cell_id=4, constants=constants_instance + ) diff --git a/virtual_rainforest/models/animals/plant_resources.py b/virtual_rainforest/models/animals/plant_resources.py index f6348c26b..b355137e6 100644 --- a/virtual_rainforest/models/animals/plant_resources.py +++ b/virtual_rainforest/models/animals/plant_resources.py @@ -34,11 +34,12 @@ def __init__(self, data: Data, cell_id: int, constants: AnimalConsts) -> None: data["layer_leaf_mass"].sel(cell_id=cell_id).sum(dim="layers").item() ) """The mass of the plant leaf mass [kg].""" - + self.constants = constants + """The animals constants.""" # Calculate energy availability # TODO - this needs to be handed back to the plants model, which will define PFT # specific conversions to different resources. - self.energy_density: float = constants.energy_density["plant"] + self.energy_density: float = self.constants.energy_density["plant"] """The energy (J) in a kg of plant [currently set to toy value of Alfalfa].""" self.energy_max: float = self.mass_current * self.energy_density """The maximum amount of energy that the cohort can have [J] [Alfalfa].""" From 15f8648c4891af97de8617c83cc1387b2e2b138f Mon Sep 17 00:00:00 2001 From: vgro Date: Thu, 30 Nov 2023 16:36:46 +0000 Subject: [PATCH 11/18] wind tests revisited and wind in model.update --- tests/conftest.py | 2 +- tests/models/abiotic/test_wind.py | 191 +++++++++----- .../models/abiotic/abiotic_model.py | 62 +++++ .../models/abiotic/abiotic_tools.py | 1 + virtual_rainforest/models/abiotic/wind.py | 240 +++++++++++++++--- .../abiotic_simple/abiotic_simple_model.py | 2 + .../models/abiotic_simple/microclimate.py | 2 + .../models/hydrology/above_ground.py | 2 + .../models/hydrology/hydrology_model.py | 2 + 9 files changed, 410 insertions(+), 94 deletions(-) diff --git a/tests/conftest.py b/tests/conftest.py index ef89a6311..889b1a5d4 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -472,7 +472,7 @@ def dummy_climate_data(layer_roles_fixture): np.full((2, 3), 450), dims=("groundwater_layers", "cell_id"), ) - data["canopy_height"] = DataArray([0, 10, 50], dims=["cell_id"]) + data["canopy_height"] = DataArray([32, 32, 32], dims=["cell_id"]) data["sensible_heat_flux_toc"] = DataArray([100, 50, 10], dims=["cell_id"]) data["friction_velocity"] = DataArray([12, 5, 2], dims=["cell_id"]) data["wind_speed_ref"] = DataArray([0, 5, 10], dims=["cell_id"]) diff --git a/tests/models/abiotic/test_wind.py b/tests/models/abiotic/test_wind.py index cc5607009..e1b527cf4 100644 --- a/tests/models/abiotic/test_wind.py +++ b/tests/models/abiotic/test_wind.py @@ -12,14 +12,13 @@ def test_calculate_zero_plane_displacement(dummy_climate_data): from virtual_rainforest.models.abiotic.wind import calculate_zero_plane_displacement - leaf_area_index = dummy_climate_data.data["leaf_area_index"].sum(dim="layers") result = calculate_zero_plane_displacement( canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(), - leaf_area_index=leaf_area_index.to_numpy(), + leaf_area_index=np.array([0, 3, 7]), zero_plane_scaling_parameter=AbioticConsts.zero_plane_scaling_parameter, ) - np.testing.assert_allclose(result, np.array([0.0, 7.910175, 39.550874])) + np.testing.assert_allclose(result, np.array([0.0, 25.312559, 27.58673])) def test_calculate_roughness_length_momentum(dummy_climate_data): @@ -29,11 +28,10 @@ def test_calculate_roughness_length_momentum(dummy_climate_data): calculate_roughness_length_momentum, ) - leaf_area_index = dummy_climate_data.data["leaf_area_index"].sum(dim="layers") result = calculate_roughness_length_momentum( canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(), - leaf_area_index=leaf_area_index.to_numpy(), - zero_plane_displacement=np.array([0, 8, 43]), + leaf_area_index=np.array([0, 3, 7]), + zero_plane_displacement=np.array([0.0, 25.312559, 27.58673]), substrate_surface_drag_coefficient=( AbioticConsts.substrate_surface_drag_coefficient ), @@ -51,7 +49,7 @@ def test_calculate_roughness_length_momentum(dummy_climate_data): ) np.testing.assert_allclose( - result, np.array([0.003, 0.434662, 1.521318]), rtol=1e-3, atol=1e-3 + result, np.array([0.017, 1.4533, 0.9591]), rtol=1e-3, atol=1e-3 ) @@ -74,7 +72,7 @@ def test_calculate_diabatic_correction_above(dummy_climate_data): sensible_heat_flux=dummy_climate_data.data["sensible_heat_flux_toc"].to_numpy(), friction_velocity=dummy_climate_data.data["friction_velocity"].to_numpy(), wind_heights=np.array([1, 15, 50]), - zero_plane_displacement=np.array([0, 8, 43]), + zero_plane_displacement=np.array([0.0, 25.312559, 27.58673]), celsius_to_kelvin=CoreConsts.zero_Celsius, von_karmans_constant=CoreConsts.von_karmans_constant, yasuda_stability_parameter1=AbioticConsts.yasuda_stability_parameter1, @@ -85,22 +83,22 @@ def test_calculate_diabatic_correction_above(dummy_climate_data): exp_result_h = np.array( [ - [0.003044, 0.026017, 0.012881], - [0.003046, 0.026031, 0.012888], - [0.003056, 0.026117, 0.01293], - [0.003073, 0.026264, 0.013002], - [0.00312, 0.026677, 0.013204], - [0.003191, 0.027289, 0.013504], + [0.003044, -0.036571, 0.042159], + [0.003046, -0.03659, 0.042182], + [0.003056, -0.036704, 0.042322], + [0.003073, -0.036902, 0.042564], + [0.00312, -0.037455, 0.043242], + [0.003191, -0.038274, 0.044247], ] ) exp_result_m = np.array( [ - [0.001827, 0.01561, 0.007729], - [0.001828, 0.015619, 0.007733], - [0.001833, 0.01567, 0.007758], - [0.001833, 0.01567, 0.007758], - [0.001872, 0.016006, 0.007923], - [0.001914, 0.016374, 0.008103], + [0.001827, -0.021943, 0.025296], + [0.001828, -0.021954, 0.025309], + [0.001833, -0.022023, 0.025393], + [0.001844, -0.022141, 0.025539], + [0.001872, -0.022473, 0.025945], + [0.001914, -0.022964, 0.026548], ] ) np.testing.assert_allclose(result["psi_h"], exp_result_h, rtol=1e-4, atol=1e-4) @@ -114,13 +112,13 @@ def test_calculate_mean_mixing_length(dummy_climate_data): result = calculate_mean_mixing_length( canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(), - zero_plane_displacement=np.array([0, 8, 43]), - roughness_length_momentum=np.array([0.003, 0.435, 1.521]), + zero_plane_displacement=np.array([0.0, 25.312559, 27.58673]), + roughness_length_momentum=np.array([0.017, 1.4533, 0.9591]), mixing_length_factor=AbioticConsts.mixing_length_factor, ) np.testing.assert_allclose( - result, np.array([0, 0.419, 1.467]), rtol=1e-3, atol=1e-3 + result, np.array([1.35804, 1.401984, 0.925228]), rtol=1e-4, atol=1e-4 ) @@ -177,44 +175,60 @@ def test_calculate_wind_attenuation_coefficient(dummy_climate_data): [17.64, 17.64, 17.64], [16.56, 16.56, 16.56], [11.16, 11.16, 11.166], + [5.76, 5.76, 5.76], + [1.17, 1.17, 1.17], + [0.414, 0.414, 0.414], ] ) result = calculate_wind_attenuation_coefficient( canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(), leaf_area_index=leaf_area_index, - mean_mixing_length=np.array([0, 0.4, 1.5]), + mean_mixing_length=np.array([1.35804, 1.401984, 0.925228]), drag_coefficient=AbioticConsts.drag_coefficient, relative_turbulence_intensity=relative_turbulence_intensity, ) exp_result = np.array( [ - [0.0, 0.141723, 0.188964], - [0.0, 0.150966, 0.201288], - [0.0, 0.224014, 0.298525], + [0.133579, 0.129392, 0.196066], + [0.142291, 0.137831, 0.208853], + [0.211141, 0.204523, 0.309744], ] ) np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) -def test_wind_log_profile(): +def test_wind_log_profile(dummy_climate_data): """Test log wind profile.""" from virtual_rainforest.models.abiotic.wind import wind_log_profile + layer_heights = ( + dummy_climate_data.data["layer_heights"] + .where(dummy_climate_data.data["layer_heights"].layer_roles != "soil") + .dropna(dim="layers") + ) + diab_correction = np.array([0.003044, -0.036571, 0.042159]) result = wind_log_profile( - height=np.array([[1, 15, 50], [5, 20, 60]]), - zeroplane_displacement=np.array([0, 8, 43]), - roughness_length_momentum=np.array([0.003, 0.435, 1.521]), - diabatic_correction_momentum=np.array([-0.1, 0.0, 0.1]), + height=layer_heights.to_numpy(), + zeroplane_displacement=np.array([0.0, 25.312559, 27.58673]), + roughness_length_momentum=np.array([0.017, 1.4533, 0.9591]), + diabatic_correction_momentum=diab_correction, ) - np.testing.assert_allclose( - result, - np.array([[5.709, 2.778, 1.627], [7.319, 3.317, 2.514]]), - rtol=1e-3, - atol=1e-3, + + exp_result = np.array( + [ + [7.543322, 1.489823, 1.568535], + [7.478785, 1.13446, 0.964925], + [7.07333, 1.0, 1.0], + [6.3802, 1.0, 1.0], + [4.483127, 1.0, 1.0], + [1.775148, 1.0, 1.0], + ] ) + np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) + def test_calculate_friction_velocity(dummy_climate_data): """Calculate friction velocity.""" @@ -223,32 +237,32 @@ def test_calculate_friction_velocity(dummy_climate_data): result = calculate_fricition_velocity( wind_speed_ref=dummy_climate_data.data["wind_speed_ref"], - reference_height=10.0, - zeroplane_displacement=np.array([0, 8, 43]), - roughness_length_momentum=np.array([0.003, 0.435, 1.521]), + reference_height=(dummy_climate_data.data["canopy_height"] + 10).to_numpy(), + zeroplane_displacement=np.array([0.0, 25.312559, 27.58673]), + roughness_length_momentum=np.array([0.017, 1.4533, 0.9591]), diabatic_correction_momentum=np.array([-0.1, 0.0, 0.1]), von_karmans_constant=CoreConsts.von_karmans_constant, ) - exp_result = np.array([0.0, 1.310997, 4.0]) - np.testing.assert_allclose(result, exp_result) + exp_result = np.array([0.0, 0.819397, 1.423534]) + np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) -def test_calculate_wind_above_canopy(): +def test_calculate_wind_above_canopy(dummy_climate_data): """Wind speed above canopy.""" from virtual_rainforest.models.abiotic.wind import calculate_wind_above_canopy result = calculate_wind_above_canopy( - friction_velocity=np.array([0, 1.3, 4]), - wind_layer_heights=np.array([2, 12, 52]), - zeroplane_displacement=np.array([0, 8, 43]), - roughness_length_momentum=np.array([0.003, 0.435, 1.521]), - diabatic_correction_momentum=np.array([-0.1, 0.0, 0.1]), + friction_velocity=np.array([0.0, 0.819397, 1.423534]), + wind_height_above=(dummy_climate_data.data["canopy_height"] + 15).to_numpy(), + zeroplane_displacement=np.array([0.0, 25.312559, 27.58673]), + roughness_length_momentum=np.array([0.017, 1.4533, 0.9591]), + diabatic_correction_momentum=np.array([0.003, 0.026, 0.013]), von_karmans_constant=CoreConsts.von_karmans_constant, min_wind_speed_above_canopy=AbioticConsts.min_wind_speed_above_canopy, ) - exp_result = np.array([0.0, 7.211, 18.779]) + exp_result = np.array([0.55, 5.590124, 10.750233]) np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) @@ -264,21 +278,84 @@ def test_calculate_wind_canopy(dummy_climate_data): ) result = calculate_wind_canopy( - top_of_canopy_wind_speed=np.array([2, 5, 7]), + top_of_canopy_wind_speed=np.array([0.55, 5.590124, 10.750233]), wind_layer_heights=layer_heights.to_numpy(), canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(), - attenuation_coefficient=np.array([0.0, 0.64, 1.24]), + attenuation_coefficient=np.array([0.133579, 0.129392, 0.196066]), min_windspeed_below_canopy=0.001, ) exp_result = np.array( [ - [0.001, 20.438858, 4.479494], - [0.001, 17.983199, 4.262731], - [0.001, 9.482404, 3.326465], - [0.001, 5.0, 2.59584], - [0.001, 2.90211, 2.102464], - [0.001, 2.65339, 2.030719], + [0.55, 5.590124, 10.750233], + [0.545427, 5.545099, 10.619302], + [0.523128, 5.325355, 9.988183], + [0.50174, 5.11432, 9.394572], + [0.48425, 4.941529, 8.917824], + [0.481428, 4.913634, 8.841656], ] ) np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) + + +def test_calculate_wind_profile(dummy_climate_data): + """Test full update of wind profile.""" + + from virtual_rainforest.models.abiotic.wind import calculate_wind_profile + + wind_layer_heights = ( + dummy_climate_data.data["layer_heights"] + .where(dummy_climate_data.data["layer_heights"].layer_roles != "soil") + .dropna(dim="layers") + ) + leaf_area_index = ( + dummy_climate_data.data["leaf_area_index"] + .where(dummy_climate_data.data["leaf_area_index"].layer_roles != "soil") + .dropna(dim="layers") + ) + air_temperature = ( + dummy_climate_data.data["air_temperature"] + .where(dummy_climate_data.data["air_temperature"].layer_roles != "soil") + .dropna(dim="layers") + ) + + wind_update = calculate_wind_profile( + canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(), + wind_height_above=(dummy_climate_data.data["canopy_height"] + 15).to_numpy(), + wind_layer_heights=wind_layer_heights.to_numpy(), + leaf_area_index=leaf_area_index.to_numpy(), + air_temperature=air_temperature.to_numpy(), + atmospheric_pressure=np.array([96, 96, 96]), + sensible_heat_flux_topofcanopy=np.array([100, 50, 10]), + wind_speed_ref=dummy_climate_data.data["wind_speed_ref"].to_numpy(), + wind_reference_height=( + dummy_climate_data.data["canopy_height"] + 10 + ).to_numpy(), + turbulence_sign=True, + abiotic_constants=AbioticConsts, + core_constants=CoreConsts, + ) + + friction_velocity_exp = np.array([0.0, 0.818637, 1.638679]) + wind_speed_exp = np.array( + [ + [0.55, 5.536364, 11.07365], + [0.54557, 5.491774, 10.984462], + [0.523951, 5.274152, 10.549181], + [0.503188, 5.065153, 10.13115], + [0.486188, 4.89403, 9.788873], + [0.483444, 4.866404, 9.733618], + ] + ) + + wind_above_exp = np.array([0.55, 5.536364, 11.07365]) + + np.testing.assert_allclose( + wind_update["wind_speed_above_canopy"], wind_above_exp, rtol=1e-3, atol=1e-3 + ) + np.testing.assert_allclose( + wind_update["friction_velocity"], friction_velocity_exp, rtol=1e-3, atol=1e-3 + ) + np.testing.assert_allclose( + wind_update["wind_speed_canopy"], wind_speed_exp, rtol=1e-3, atol=1e-3 + ) diff --git a/virtual_rainforest/models/abiotic/abiotic_model.py b/virtual_rainforest/models/abiotic/abiotic_model.py index 347ab46e4..d5dc9c2c7 100644 --- a/virtual_rainforest/models/abiotic/abiotic_model.py +++ b/virtual_rainforest/models/abiotic/abiotic_model.py @@ -12,20 +12,26 @@ class as a child of the :class:`~virtual_rainforest.core.base_model.BaseModel` c config dictionary to the required types they are caught and then logged, and at the end of the unpacking an error is thrown. This error should be caught and handled by downstream functions so that all model configuration failures can be reported as one. + +TODO temperatures in Kelvin """ # noqa: D205, D415 from __future__ import annotations from typing import Any +import numpy as np from pint import Quantity +from xarray import DataArray from virtual_rainforest.core.base_model import BaseModel from virtual_rainforest.core.config import Config +from virtual_rainforest.core.constants import CoreConsts from virtual_rainforest.core.constants_loader import load_constants from virtual_rainforest.core.data import Data from virtual_rainforest.core.logger import LOGGER from virtual_rainforest.core.utils import set_layer_roles +from virtual_rainforest.models.abiotic import wind from virtual_rainforest.models.abiotic.constants import AbioticConsts @@ -122,5 +128,61 @@ def update(self, time_index: int, **kwargs: Any) -> None: time_index: The index of the current time step in the data object. """ + wind_update_inputs: dict[str, DataArray] = {} + + for var in ["layer_heights", "leaf_area_index", "air_temperature"]: + selection = ( + self.data[var] + .where(self.data[var].layer_roles != "soil") + .dropna(dim="layers") + ) + wind_update_inputs[var].append(selection[var]) + + wind_update = wind.calculate_wind_profile( + canopy_height=self.data["canopy_height"].to_numpy(), + wind_height_above=(self.data["canopy_height"] + 15).to_numpy(), # TODO + wind_layer_heights=wind_update_inputs["layer_heights"].to_numpy(), + leaf_area_index=wind_update_inputs["leaf_area_index"].to_numpy(), + air_temperature=wind_update_inputs["air_temperature"].to_numpy(), + atmospheric_pressure=self.data["atmospheric_pressure_ref"].to_numpy(), + sensible_heat_flux_topofcanopy=( + self.data["sensible_heat_flux_topofcanopy"].to_numpy() + ), + wind_speed_ref=self.data["wind_speed_ref"].to_numpy(), + wind_reference_height=(self.data["canopy_height"] + 10).to_numpy(), + turbulence_sign=True, + abiotic_constants=AbioticConsts(), + core_constants=CoreConsts(), + ) + + wind_output = {} + wind_speed_above_canopy = DataArray( + wind_update["wind_speed_above_canopy"], + dims="cell_id", + coords={"cell_id": self.data.grid.cell_id}, + ) + wind_output["wind_speed_above_canopy"] = wind_speed_above_canopy + + for var in ["wind_speed_canopy", "friction_velocity"]: + var_out = DataArray( + np.concatenate( + ( + wind_update[var], + np.full( + ( + self.layer_roles.count("soil"), + self.data.grid.n_cells, + ), + np.nan, + ), + ), + ), + dims=self.data["layer_heights"].dims, + coords=self.data["layer_heights"].coords, + ) + wind_output[var] = var_out + + self.data.add_from_dict(output_dict=wind_output) + def cleanup(self) -> None: """Placeholder function for abiotic model cleanup.""" diff --git a/virtual_rainforest/models/abiotic/abiotic_tools.py b/virtual_rainforest/models/abiotic/abiotic_tools.py index 758fac76d..62aae0a2a 100644 --- a/virtual_rainforest/models/abiotic/abiotic_tools.py +++ b/virtual_rainforest/models/abiotic/abiotic_tools.py @@ -2,6 +2,7 @@ are shared across submodules in the :mod:`~virtual_rainforest.models.abiotic` model. TODO cross-check with pyrealm for duplication/ different implementation +TODO change temperatures to Kelvin """ # noqa: D205, D415 import numpy as np diff --git a/virtual_rainforest/models/abiotic/wind.py b/virtual_rainforest/models/abiotic/wind.py index 83ddf5d48..d0ed15e2b 100644 --- a/virtual_rainforest/models/abiotic/wind.py +++ b/virtual_rainforest/models/abiotic/wind.py @@ -3,31 +3,10 @@ :math:`CO_{2}` between soil and atmosphere below the canopy as well as the exchange with the atmosphere above the canopy. -The wind profile above the canopy is described as follows (based on -:cite:p:`campbell_introduction_1998` as implemented in -:cite:t:`maclean_microclimc_2021`): - -:math:`u_z = \frac{u^{*}}{0.4} ln \frac{z-d}{z_M} + \phi_M` - -where :math:`u_z` is wind speed at height :math:`z` above the canopy, :math:`d` is -the height above ground within the canopy where the wind profile extrapolates to -zero, :math:`z_m` the roughness length for momentum, :math:`\phi_M` is a diabatic -correction for momentum and :math:`u^{*}` is the friction velocity, which gives the -wind speed at height :math:`d + z_m`. - -The wind profile below canopy is derived as follows: - -:math:`u_z = u_h exp(a(\frac{z}{h} - 1))` - -where :math:`u_z` is wind speed at height :math:`z` within the canopy, :math:`u_h` -is wind speed at the top of the canopy at height :math:`h`, and :math:`a` is a wind -attenuation coefficient given by :math:`a = 2 l_m i_w`, where :math:`c_d` is a drag -coefficient that varies with leaf inclination and shape, :math:`i_w` is a -coefficient describing relative turbulence intensity and :math:`l_m` is the mean -mixing length, equivalent to the free space between the leaves and stems. For -details, see :cite:t:`maclean_microclimc_2021`. - TODO: add sanity checks, errors and logging +TODO replace leaf area index by plant area index when we have more info about vertical +distribution of leaf and woody parts +TODO change temperatures to Kelvin """ # noqa: D205, D415 from typing import Union @@ -35,6 +14,13 @@ import numpy as np from numpy.typing import NDArray +from virtual_rainforest.core.constants import CoreConsts +from virtual_rainforest.models.abiotic.abiotic_tools import ( + calculate_molar_density_air, + calculate_specific_heat_air, +) +from virtual_rainforest.models.abiotic.constants import AbioticConsts + def calculate_zero_plane_displacement( canopy_height: NDArray[np.float32], @@ -313,7 +299,11 @@ def calculate_wind_attenuation_coefficient( intermediate_coefficient = ( (drag_coefficient * leaf_area_index * canopy_height) - / (2 * mean_mixing_length * relative_turbulence_intensity), + / ( + 2 + * mean_mixing_length + * relative_turbulence_intensity[0 : len(leaf_area_index)] + ), ) return np.nan_to_num(intermediate_coefficient, nan=0).squeeze() @@ -348,10 +338,10 @@ def wind_log_profile( def calculate_fricition_velocity( wind_speed_ref: NDArray[np.float32], - reference_height: float, + reference_height: Union[float, NDArray[np.float32]], zeroplane_displacement: NDArray[np.float32], roughness_length_momentum: NDArray[np.float32], - diabatic_correction_momentum: NDArray[np.float32], + diabatic_correction_momentum: Union[float, NDArray[np.float32]], von_karmans_constant: float, ) -> NDArray[np.float32]: """Calculate friction velocity from wind speed at reference height. @@ -384,7 +374,7 @@ def calculate_fricition_velocity( def calculate_wind_above_canopy( friction_velocity: NDArray[np.float32], - wind_layer_heights: NDArray[np.float32], + wind_height_above: NDArray[np.float32], zeroplane_displacement: NDArray[np.float32], roughness_length_momentum: NDArray[np.float32], diabatic_correction_momentum: NDArray[np.float32], @@ -404,8 +394,8 @@ def calculate_wind_above_canopy( :cite:t:`maclean_microclimc_2021`. Args: - friction_velocity: friction velocity - wind_layer_heights: Heights for which wind speed is required, [m] + friction_velocity: friction velocity, [m s-1] + wind_height_above: Height above canopy for which wind speed is required, [m] zero_plane_displacement: Height above ground within the canopy where the wind profile extrapolates to zero, [m] roughness_length_momentum: Momentum roughness length, [m] @@ -420,19 +410,19 @@ def calculate_wind_above_canopy( """ wind_profile_above = wind_log_profile( - height=wind_layer_heights, + height=wind_height_above, zeroplane_displacement=zeroplane_displacement, roughness_length_momentum=roughness_length_momentum, diabatic_correction_momentum=diabatic_correction_momentum, ) - np.where( - wind_profile_above < min_wind_speed_above_canopy, + wind_profile = (friction_velocity / von_karmans_constant) * wind_profile_above + + return np.where( + wind_profile < min_wind_speed_above_canopy, min_wind_speed_above_canopy, - wind_profile_above, + wind_profile, ) - return (friction_velocity / von_karmans_constant) * wind_profile_above - def calculate_wind_canopy( top_of_canopy_wind_speed: NDArray[np.float32], @@ -464,3 +454,181 @@ def calculate_wind_canopy( * np.exp(attenuation_coefficient * ((wind_layer_heights / canopy_height) - 1)), nan=min_windspeed_below_canopy, ).squeeze() + + +def calculate_wind_profile( + canopy_height: NDArray[np.float32], + wind_height_above: NDArray[np.float32], + wind_layer_heights: NDArray[np.float32], + leaf_area_index: NDArray[np.float32], + air_temperature: NDArray[np.float32], + atmospheric_pressure: NDArray[np.float32], + sensible_heat_flux_topofcanopy: NDArray[np.float32], + wind_speed_ref: NDArray[np.float32], + wind_reference_height: Union[float, NDArray[np.float32]], + turbulence_sign: bool, + abiotic_constants: AbioticConsts, + core_constants: CoreConsts, +) -> dict[str, NDArray[np.float32]]: + r"""Calculate wind speed above and below the canoopy, [m s-1]. + + The wind profile above the canopy is described as follows (based on + :cite:p:`campbell_introduction_1998` as implemented in + :cite:t:`maclean_microclimc_2021`): + + :math:`u_z = \frac{u^{*}}{0.4} ln \frac{z-d}{z_M} + \phi_M` + + where :math:`u_z` is wind speed at height :math:`z` above the canopy, :math:`d` is + the height above ground within the canopy where the wind profile extrapolates to + zero, :math:`z_m` the roughness length for momentum, :math:`\phi_M` is a diabatic + correction for momentum and :math:`u^{*}` is the friction velocity, which gives the + wind speed at height :math:`d + z_m`. + + The wind profile below canopy is derived as follows: + + :math:`u_z = u_h exp(a(\frac{z}{h} - 1))` + + where :math:`u_z` is wind speed at height :math:`z` within the canopy, :math:`u_h` + is wind speed at the top of the canopy at height :math:`h`, and :math:`a` is a wind + attenuation coefficient given by :math:`a = 2 l_m i_w`, where :math:`c_d` is a drag + coefficient that varies with leaf inclination and shape, :math:`i_w` is a + coefficient describing relative turbulence intensity and :math:`l_m` is the mean + mixing length, equivalent to the free space between the leaves and stems. For + details, see :cite:t:`maclean_microclimc_2021`. + + Args: + + Returns: + dictionnary that contains wind speed above the canopy, [m s-1], wind speed + within and below the canopy, [m s-1], and friction velocity, [m s-1] + """ + + output = {} + + # TODO adjust wind to 2m above canopy? + + molar_density_air = calculate_molar_density_air( + temperature=air_temperature, + atmospheric_pressure=atmospheric_pressure, + standard_mole=core_constants.standard_mole, + standard_pressure=core_constants.standard_pressure, + celsius_to_kelvin=core_constants.zero_Celsius, + ) + + specific_heat_air = calculate_specific_heat_air( + temperature=air_temperature, + molar_heat_capacity_air=core_constants.molar_heat_capacity_air, + specific_heat_equ_factor_1=abiotic_constants.specific_heat_equ_factor_1, + specific_heat_equ_factor_2=abiotic_constants.specific_heat_equ_factor_2, + ) + + leaf_area_index_sum = leaf_area_index.sum(axis=1) + + zero_plane_displacement = calculate_zero_plane_displacement( + canopy_height=canopy_height, + leaf_area_index=leaf_area_index_sum, + zero_plane_scaling_parameter=abiotic_constants.zero_plane_scaling_parameter, + ) + + roughness_length_momentum = calculate_roughness_length_momentum( + canopy_height=canopy_height, + leaf_area_index=leaf_area_index_sum, + zero_plane_displacement=zero_plane_displacement, + substrate_surface_drag_coefficient=( + AbioticConsts.substrate_surface_drag_coefficient + ), + roughness_element_drag_coefficient=( + AbioticConsts.roughness_element_drag_coefficient + ), + roughness_sublayer_depth_parameter=( + AbioticConsts.roughness_sublayer_depth_parameter + ), + max_ratio_wind_to_friction_velocity=( + AbioticConsts.max_ratio_wind_to_friction_velocity + ), + min_roughness_length=AbioticConsts.min_roughness_length, + von_karman_constant=CoreConsts.von_karmans_constant, + ) + + friction_velocity_uncorrected = calculate_fricition_velocity( + wind_speed_ref=wind_speed_ref, + reference_height=wind_reference_height, + zeroplane_displacement=zero_plane_displacement, + roughness_length_momentum=roughness_length_momentum, + diabatic_correction_momentum=0.0, + von_karmans_constant=core_constants.von_karmans_constant, + ) + + diabatic_correction_above = calculate_diabatic_correction_above( + molar_density_air=molar_density_air, + specific_heat_air=specific_heat_air, + temperature=air_temperature, + sensible_heat_flux=sensible_heat_flux_topofcanopy, + friction_velocity=friction_velocity_uncorrected, + wind_heights=wind_layer_heights, + zero_plane_displacement=zero_plane_displacement, + celsius_to_kelvin=CoreConsts.zero_Celsius, + von_karmans_constant=CoreConsts.von_karmans_constant, + yasuda_stability_parameter1=AbioticConsts.yasuda_stability_parameter1, + yasuda_stability_parameter2=AbioticConsts.yasuda_stability_parameter2, + yasuda_stability_parameter3=AbioticConsts.yasuda_stability_parameter3, + diabatic_heat_momentum_ratio=AbioticConsts.diabatic_heat_momentum_ratio, + ) + + friction_velocity = calculate_fricition_velocity( + wind_speed_ref=wind_speed_ref, + reference_height=wind_reference_height, + zeroplane_displacement=zero_plane_displacement, + roughness_length_momentum=roughness_length_momentum, + diabatic_correction_momentum=diabatic_correction_above["psi_m"], + von_karmans_constant=core_constants.von_karmans_constant, + ) + output["friction_velocity"] = friction_velocity[0] + + mean_mixing_length = calculate_mean_mixing_length( + canopy_height=canopy_height, + zero_plane_displacement=zero_plane_displacement, + roughness_length_momentum=roughness_length_momentum, + mixing_length_factor=abiotic_constants.mixing_length_factor, + ) + + relative_turbulence_intensity = generate_relative_turbulence_intensity( + layer_heights=wind_layer_heights, + min_relative_turbulence_intensity=( + abiotic_constants.min_relative_turbulence_intensity + ), + max_relative_turbulence_intensity=( + abiotic_constants.max_relative_turbulence_intensity + ), + increasing_with_height=turbulence_sign, + ) + + attennuation_coefficient = calculate_wind_attenuation_coefficient( + canopy_height=canopy_height, + leaf_area_index=leaf_area_index, + mean_mixing_length=mean_mixing_length, + drag_coefficient=abiotic_constants.drag_coefficient, + relative_turbulence_intensity=relative_turbulence_intensity, + ) + + wind_speed_above_canopy = calculate_wind_above_canopy( + friction_velocity=friction_velocity[0], + wind_height_above=wind_height_above, + zeroplane_displacement=zero_plane_displacement, + roughness_length_momentum=roughness_length_momentum, + diabatic_correction_momentum=diabatic_correction_above["psi_m"][0], + von_karmans_constant=core_constants.von_karmans_constant, + min_wind_speed_above_canopy=abiotic_constants.min_wind_speed_above_canopy, + ) + output["wind_speed_above_canopy"] = wind_speed_above_canopy + + wind_speed_canopy = calculate_wind_canopy( + top_of_canopy_wind_speed=wind_speed_above_canopy, + wind_layer_heights=wind_layer_heights, + canopy_height=canopy_height, + attenuation_coefficient=attennuation_coefficient[0], + min_windspeed_below_canopy=abiotic_constants.min_windspeed_below_canopy, + ) + output["wind_speed_canopy"] = wind_speed_canopy + + return output diff --git a/virtual_rainforest/models/abiotic_simple/abiotic_simple_model.py b/virtual_rainforest/models/abiotic_simple/abiotic_simple_model.py index 1bd888263..9a944aa73 100644 --- a/virtual_rainforest/models/abiotic_simple/abiotic_simple_model.py +++ b/virtual_rainforest/models/abiotic_simple/abiotic_simple_model.py @@ -13,6 +13,8 @@ class as a child of the :class:`~virtual_rainforest.core.base_model.BaseModel` c config dictionary to the required types they are caught and then logged, and at the end of the unpacking an error is thrown. This error should be caught and handled by downstream functions so that all model configuration failures can be reported as one. + +TODO update temperatures to Kelvin """ # noqa: D205, D415 from __future__ import annotations diff --git a/virtual_rainforest/models/abiotic_simple/microclimate.py b/virtual_rainforest/models/abiotic_simple/microclimate.py index b045d3d44..765344034 100644 --- a/virtual_rainforest/models/abiotic_simple/microclimate.py +++ b/virtual_rainforest/models/abiotic_simple/microclimate.py @@ -8,6 +8,8 @@ 1 m depth which equals the mean annual temperature. The module also provides a constant vertical profile of atmospheric pressure and :math:`\ce{CO2}`. + +TODO change tenperatures to Kelvin """ # noqa: D205, D415 import numpy as np diff --git a/virtual_rainforest/models/hydrology/above_ground.py b/virtual_rainforest/models/hydrology/above_ground.py index 11741a423..37ec45c85 100644 --- a/virtual_rainforest/models/hydrology/above_ground.py +++ b/virtual_rainforest/models/hydrology/above_ground.py @@ -2,6 +2,8 @@ processes for the Virtual Rainforest. At the moment, this includes rain water interception by the canopy, soil evaporation, and functions related to surface runoff, bypass flow, and river discharge. + +TODO change temperatures to Kelvin """ # noqa: D205, D415 from math import sqrt diff --git a/virtual_rainforest/models/hydrology/hydrology_model.py b/virtual_rainforest/models/hydrology/hydrology_model.py index 14c1b4ca2..b443adb8d 100644 --- a/virtual_rainforest/models/hydrology/hydrology_model.py +++ b/virtual_rainforest/models/hydrology/hydrology_model.py @@ -13,6 +13,8 @@ class as a child of the :class:`~virtual_rainforest.core.base_model.BaseModel` c config dictionary to the required types they are caught and then logged, and at the end of the unpacking an error is thrown. This error should be caught and handled by downstream functions so that all model configuration failures can be reported as one. + +TODO change temperature to Kelvin """ # noqa: D205, D415 from __future__ import annotations From 8e109493592981ac264c603f3c09288800b82df8 Mon Sep 17 00:00:00 2001 From: vgro Date: Mon, 4 Dec 2023 11:57:01 +0000 Subject: [PATCH 12/18] required vars in abiotic model and tests updated --- tests/conftest.py | 2 +- tests/models/abiotic/test_abiotic_model.py | 224 +++++++++++++++--- tests/models/abiotic/test_wind.py | 4 +- .../models/abiotic/abiotic_model.py | 10 +- 4 files changed, 200 insertions(+), 40 deletions(-) diff --git a/tests/conftest.py b/tests/conftest.py index 889b1a5d4..35130a930 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -473,7 +473,7 @@ def dummy_climate_data(layer_roles_fixture): dims=("groundwater_layers", "cell_id"), ) data["canopy_height"] = DataArray([32, 32, 32], dims=["cell_id"]) - data["sensible_heat_flux_toc"] = DataArray([100, 50, 10], dims=["cell_id"]) + data["sensible_heat_flux_topofcanopy"] = DataArray([100, 50, 10], dims=["cell_id"]) data["friction_velocity"] = DataArray([12, 5, 2], dims=["cell_id"]) data["wind_speed_ref"] = DataArray([0, 5, 10], dims=["cell_id"]) return data diff --git a/tests/models/abiotic/test_abiotic_model.py b/tests/models/abiotic/test_abiotic_model.py index bea3b7b2d..dcee82858 100644 --- a/tests/models/abiotic/test_abiotic_model.py +++ b/tests/models/abiotic/test_abiotic_model.py @@ -1,7 +1,7 @@ """Test module for abiotic.abiotic_model.py.""" from contextlib import nullcontext as does_not_raise -from logging import DEBUG, ERROR, INFO +from logging import CRITICAL, DEBUG, ERROR, INFO import pint import pytest @@ -10,25 +10,9 @@ from virtual_rainforest.core.exceptions import ConfigurationError -@pytest.mark.parametrize( - "raises,expected_log_entries", - [ - ( - does_not_raise(), - ( - ( - DEBUG, - "abiotic model: required var 'air_temperature_ref' checked", - ), - ), - ), - ], -) def test_abiotic_model_initialization( caplog, dummy_climate_data, - raises, - expected_log_entries, layer_roles_fixture, ): """Test `AbioticModel` initialization.""" @@ -36,47 +20,213 @@ def test_abiotic_model_initialization( from virtual_rainforest.models.abiotic.abiotic_model import AbioticModel from virtual_rainforest.models.abiotic.constants import AbioticConsts - with raises: - # Initialize model - model = AbioticModel( - dummy_climate_data, + # Initialize model + model = AbioticModel( + dummy_climate_data, + pint.Quantity("1 day"), + soil_layers=[-0.5, -1.0], + canopy_layers=10, + constants=AbioticConsts(), + ) + + # In cases where it passes then checks that the object has the right properties + assert isinstance(model, BaseModel) + assert model.model_name == "abiotic" + assert str(model) == "A abiotic model instance" + assert repr(model) == "AbioticModel(update_interval = 1 day)" + assert model.layer_roles == layer_roles_fixture + + # Final check that expected logging entries are produced + log_check( + caplog, + expected_log=( + (DEBUG, "abiotic model: required var 'air_temperature' checked"), + (DEBUG, "abiotic model: required var 'canopy_height' checked"), + (DEBUG, "abiotic model: required var 'layer_heights' checked"), + (DEBUG, "abiotic model: required var 'leaf_area_index' checked"), + (DEBUG, "abiotic model: required var 'atmospheric_pressure_ref' checked"), + ( + DEBUG, + ( + "abiotic model: required var 'sensible_heat_flux_topofcanopy'" + " checked" + ), + ), + (DEBUG, "abiotic model: required var 'wind_speed_ref' checked"), + ), + ) + + +def test_abiotic_model_initialization_no_data(caplog, dummy_climate_data): + """Test `AbioticModel` initialization with no data.""" + + from virtual_rainforest.core.data import Data + from virtual_rainforest.core.grid import Grid + from virtual_rainforest.models.abiotic.abiotic_model import AbioticModel + from virtual_rainforest.models.abiotic.constants import AbioticConsts + + with pytest.raises(ValueError): + # Make four cell grid + grid = Grid(cell_nx=4, cell_ny=1) + empty_data = Data(grid) + + # Try and initialise model with empty data object + _ = AbioticModel( + empty_data, pint.Quantity("1 day"), soil_layers=[-0.5, -1.0], canopy_layers=10, - constants=AbioticConsts(), + constants=AbioticConsts, ) - # In cases where it passes then checks that the object has the right properties - assert isinstance(model, BaseModel) - assert model.model_name == "abiotic" - assert repr(model) == "AbioticModel(update_interval = 1 day)" - assert model.layer_roles == layer_roles_fixture - # Final check that expected logging entries are produced - log_check(caplog, expected_log_entries) + log_check( + caplog, + expected_log=( + (ERROR, "abiotic model: init data missing required var 'air_temperature'"), + (ERROR, "abiotic model: init data missing required var 'canopy_height'"), + (ERROR, "abiotic model: init data missing required var 'layer_heights'"), + (ERROR, "abiotic model: init data missing required var 'leaf_area_index'"), + ( + ERROR, + ( + "abiotic model: init data missing required var" + " 'atmospheric_pressure_ref'" + ), + ), + ( + ERROR, + ( + "abiotic model: init data missing required var" + " 'sensible_heat_flux_topofcanopy'" + ), + ), + (ERROR, "abiotic model: init data missing required var 'wind_speed_ref'"), + (ERROR, "abiotic model: error checking required_init_vars, see log."), + ), + ) @pytest.mark.parametrize( - "cfg_string,time_interval,raises,expected_log_entries", + "cfg_string, time_interval, drag_coeff, raises, expected_log_entries", [ pytest.param( - "[core]\n[core.timing]\nupdate_interval = '1 day'\n[abiotic]\n", - pint.Quantity("1 day"), + "[core]\n[core.timing]\nupdate_interval = '12 hours'\n[abiotic]\n", + pint.Quantity("12 hours"), + 0.2, does_not_raise(), ( (INFO, "Initialised abiotic.AbioticConsts from config"), ( INFO, - "Information required to initialise the abiotic model " - "successfully extracted.", + "Information required to initialise the abiotic model successfully " + "extracted.", + ), + (DEBUG, "abiotic model: required var 'air_temperature' checked"), + (DEBUG, "abiotic model: required var 'canopy_height' checked"), + (DEBUG, "abiotic model: required var 'layer_heights' checked"), + (DEBUG, "abiotic model: required var 'leaf_area_index' checked"), + ( + DEBUG, + "abiotic model: required var 'atmospheric_pressure_ref' checked", + ), + ( + DEBUG, + ( + "abiotic model: required var 'sensible_heat_flux_topofcanopy'" + " checked" + ), + ), + (DEBUG, "abiotic model: required var 'wind_speed_ref' checked"), + ), + id="default_config", + ), + pytest.param( + "[core]\n[core.timing]\nupdate_interval = '12 hours'\n" + "[abiotic.constants.AbioticConsts]\ndrag_coefficient = 0.05\n", + pint.Quantity("12 hours"), + 0.05, + does_not_raise(), + ( + (INFO, "Initialised abiotic.AbioticConsts from config"), + ( + INFO, + "Information required to initialise the abiotic model successfully " + "extracted.", ), + (DEBUG, "abiotic model: required var 'air_temperature' checked"), + (DEBUG, "abiotic model: required var 'canopy_height' checked"), + (DEBUG, "abiotic model: required var 'layer_heights' checked"), + (DEBUG, "abiotic model: required var 'leaf_area_index' checked"), ( DEBUG, - "abiotic model: required var 'air_temperature_ref' checked", + "abiotic model: required var 'atmospheric_pressure_ref' checked", ), + ( + DEBUG, + ( + "abiotic model: required var 'sensible_heat_flux_topofcanopy'" + " checked" + ), + ), + (DEBUG, "abiotic model: required var 'wind_speed_ref' checked"), ), - id="initialises correctly", + id="modified_config_correct", ), + pytest.param( + "[core]\n[core.timing]\nupdate_interval = '12 hours'\n" + "[abiotic.constants.AbioticConsts]\ndrag_coefficients = 0.05\n", + None, + None, + pytest.raises(ConfigurationError), + ( + (ERROR, "Unknown names supplied for AbioticConsts: drag_coefficients"), + (INFO, "Valid names are: "), + (CRITICAL, "Could not initialise abiotic.AbioticConsts from config"), + ), + id="modified_config_incorrect", + ), + ], +) +def test_generate_abiotic_model( + caplog, + dummy_climate_data, + cfg_string, + time_interval, + drag_coeff, + raises, + expected_log_entries, +): + """Test that the function to initialise the abiotic model behaves as expected.""" + + from virtual_rainforest.core.config import Config + from virtual_rainforest.core.registry import register_module + from virtual_rainforest.models.abiotic.abiotic_model import AbioticModel + + # Register the module components to access constants classes + register_module("virtual_rainforest.models.abiotic") + + # Build the config object + config = Config(cfg_strings=cfg_string) + caplog.clear() + + # Check whether model is initialised (or not) as expected + with raises: + model = AbioticModel.from_config( + dummy_climate_data, + config, + pint.Quantity(config["core"]["timing"]["update_interval"]), + ) + assert model.update_interval == time_interval + assert model.constants.drag_coefficient == drag_coeff + + # Final check that expected logging entries are produced + log_check(caplog, expected_log_entries) + + +@pytest.mark.parametrize( + "cfg_string, time_interval, raises, expected_log_entries", + [ pytest.param( "[core]\n[core.timing]\nupdate_interval = '1 month'\n[abiotic]\n", pint.Quantity("1 month"), @@ -94,7 +244,7 @@ def test_abiotic_model_initialization( ), ], ) -def test_generate_abiotic_model( +def test_generate_abiotic_model_bounds_error( caplog, dummy_climate_data, cfg_string, diff --git a/tests/models/abiotic/test_wind.py b/tests/models/abiotic/test_wind.py index e1b527cf4..eddf6f98c 100644 --- a/tests/models/abiotic/test_wind.py +++ b/tests/models/abiotic/test_wind.py @@ -69,7 +69,9 @@ def test_calculate_diabatic_correction_above(dummy_climate_data): molar_density_air=np.repeat(28.96, 3), specific_heat_air=np.repeat(1, 3), temperature=air_temperature.to_numpy(), - sensible_heat_flux=dummy_climate_data.data["sensible_heat_flux_toc"].to_numpy(), + sensible_heat_flux=( + dummy_climate_data.data["sensible_heat_flux_topofcanopy"].to_numpy() + ), friction_velocity=dummy_climate_data.data["friction_velocity"].to_numpy(), wind_heights=np.array([1, 15, 50]), zero_plane_displacement=np.array([0.0, 25.312559, 27.58673]), diff --git a/virtual_rainforest/models/abiotic/abiotic_model.py b/virtual_rainforest/models/abiotic/abiotic_model.py index d5dc9c2c7..8aa1afb30 100644 --- a/virtual_rainforest/models/abiotic/abiotic_model.py +++ b/virtual_rainforest/models/abiotic/abiotic_model.py @@ -52,7 +52,15 @@ class AbioticModel(BaseModel): """Shortest time scale that abiotic model can sensibly capture.""" upper_bound_on_time_scale = "1 day" """Longest time scale that abiotic model can sensibly capture.""" - required_init_vars = (("air_temperature_ref", ("spatial",)),) + required_init_vars = ( + ("air_temperature", ("spatial",)), + ("canopy_height", ("spatial",)), + ("layer_heights", ("spatial",)), + ("leaf_area_index", ("spatial",)), + ("atmospheric_pressure_ref", ("spatial",)), + ("sensible_heat_flux_topofcanopy", ("spatial",)), + ("wind_speed_ref", ("spatial",)), + ) """The required variables and axes for the abiotic model""" vars_updated = () """Variables updated by the abiotic model""" From 2c1a01ba97114b81aa75c9dc6eabb058ebde75c7 Mon Sep 17 00:00:00 2001 From: vgro Date: Mon, 4 Dec 2023 12:12:53 +0000 Subject: [PATCH 13/18] updated docstrings --- virtual_rainforest/models/abiotic/wind.py | 78 +++++++++++++---------- 1 file changed, 43 insertions(+), 35 deletions(-) diff --git a/virtual_rainforest/models/abiotic/wind.py b/virtual_rainforest/models/abiotic/wind.py index d0ed15e2b..2e5dc9395 100644 --- a/virtual_rainforest/models/abiotic/wind.py +++ b/virtual_rainforest/models/abiotic/wind.py @@ -27,7 +27,7 @@ def calculate_zero_plane_displacement( leaf_area_index: NDArray[np.float32], zero_plane_scaling_parameter: float, ) -> NDArray[np.float32]: - """Calculate zero plane displacement. + """Calculate zero plane displacement height, [m]. The zero plane displacement height of a vegetated surface is the height at which the wind speed would go to zero if the logarithmic wind profile was maintained from the @@ -38,7 +38,7 @@ def calculate_zero_plane_displacement( Args: canopy_height: Canopy height, [m] leaf_area_index: Total leaf area index, [m m-1] - zero_plane_scaling_parameter: Control parameter for scaling d/h + zero_plane_scaling_parameter: Control parameter for scaling d/h, dimensionless :cite:p:`raupach_simplified_1994` Returns: @@ -69,7 +69,7 @@ def calculate_roughness_length_momentum( min_roughness_length: float, von_karman_constant: float, ) -> NDArray[np.float32]: - """Calculate roughness length governing momentum transfer. + """Calculate roughness length governing momentum transfer, [m]. Roughness length is defined as the height at which the mean velocity is zero due to substrate roughness. Real surfaces such as the ground or vegetation are not smooth @@ -81,15 +81,17 @@ def calculate_roughness_length_momentum( leaf_area_index: Total leaf area index, [m m-1] zero_plane_displacement: Height above ground within the canopy where the wind profile extrapolates to zero, [m] - substrate_surface_drag_coefficient: Substrate-surface drag coefficient + substrate_surface_drag_coefficient: Substrate-surface drag coefficient, + dimensionless roughness_element_drag_coefficient: Roughness-element drag coefficient roughness_sublayer_depth_parameter: Parameter that characterizes the roughness - sublayer depth + sublayer depth, dimensionless max_ratio_wind_to_friction_velocity: Maximum ratio of wind velocity to friction - velocity + velocity, dimensionless min_roughness_length: Minimum roughness length, [m] - von_karman_constant: Von Karman's constant, unitless constant describing the - logarithmic velocity profile of a turbulent fluid near a no-slip boundary. + von_karman_constant: Von Karman's constant, dimensionless constant describing + the logarithmic velocity profile of a turbulent fluid near a no-slip + boundary. Returns: momentum roughness length, [m] @@ -141,7 +143,7 @@ def calculate_diabatic_correction_above( yasuda_stability_parameter3: float, diabatic_heat_momentum_ratio: float, ) -> dict[str, NDArray[np.float32]]: - r"""Calculates the diabatic correction factors for momentum and heat above canopy. + r"""Calculate the diabatic correction factors for momentum and heat above canopy. Diabatic correction factor for heat and momentum are used to adjust wind profiles for surface heating and cooling :cite:p:`maclean_microclimc_2021`. When the surface @@ -152,16 +154,18 @@ def calculate_diabatic_correction_above( Args: molar_density_air: molar density of air, [mol m-3] specific_heat_air: specific heat of air, [J mol-1 K-1] - temperature: 2m temperature + temperature: 2m temperature, [C] sensible_heat_flux: Sensible heat flux from canopy to atmosphere above, [W m-2], # TODO: could be the top entry of the general sensible heat flux - friction_velocity: Friction velocity + friction_velocity: Friction velocity, [m s-1] wind_heights: Vector of heights for which wind speed is calculated, [m] zero_plane_displacement: Height above ground within the canopy where the wind profile extrapolates to zero, [m] celsius_to_kelvin: Factor to convert temperature in Celsius to absolute temperature in Kelvin - von_karmans_constant + von_karmans_constant: Von Karman's constant, dimensionless constant describing + the logarithmic velocity profile of a turbulent fluid near a no-slip + boundary. yasuda_stability_parameter1: Parameter to approximate diabatic correction factors for heat and momentum after :cite:t:`yasuda_turbulent_1988` yasuda_stability_parameter2: Parameter to approximate diabatic correction @@ -211,7 +215,7 @@ def calculate_mean_mixing_length( roughness_length_momentum: NDArray[np.float32], mixing_length_factor: float, ) -> NDArray[np.float32]: - """Calculate mixing length for canopy air transport. + """Calculate mixing length for canopy air transport, [m]. The mean mixing length is used to calculate turbulent air transport inside vegetated canopies. It is made equivalent to the above canopy value at the canopy surface. In @@ -223,7 +227,7 @@ def calculate_mean_mixing_length( zero_plane_displacement: Height above ground within the canopy where the wind profile extrapolates to zero, [m] roughness_length_momentum: Momentum roughness length, [m] - mixing_length_factor: Factor in calculation of mean mixing length + mixing_length_factor: Factor in calculation of mean mixing length, dimensionless Returns: mixing length for canopy air transport, [m] @@ -242,20 +246,22 @@ def generate_relative_turbulence_intensity( max_relative_turbulence_intensity: float, increasing_with_height: bool, ) -> NDArray[np.float32]: - """Generate relative turbulence intensity profile. + """Generate relative turbulence intensity profile, dimensionless. At the moment, default values are for a maize crop Shaw et al (1974) - Agricultural Meteorology, 13: 419-425. TODO adjust to environment + Agricultural Meteorology, 13: 419-425. TODO adjust default to environment Args: - layer_heights: heights of above ground layers, m - min_relative_turbulence_intensity: minimum relative turbulence intensity - max_relative_turbulence_intensity: maximum relative turbulence intensity + layer_heights: heights of above ground layers, [m] + min_relative_turbulence_intensity: minimum relative turbulence intensity, + dimensionless + max_relative_turbulence_intensity: maximum relative turbulence intensity, + dimensionless increasing_with_height: increasing logical indicating whether turbulence intensity increases (TRUE) or decreases (FALSE) with height Returns: - relative turbulence intensity for each node + relative turbulence intensity for each node, dimensionless """ if increasing_with_height: @@ -280,7 +286,7 @@ def calculate_wind_attenuation_coefficient( drag_coefficient: float, relative_turbulence_intensity: NDArray[np.float32], ) -> NDArray[np.float32]: - """Calculate wind attenuation coefficient. + """Calculate wind attenuation coefficient, dimensionless. The wind attenuation coefficient describes how wind is slowed down by the presence of vegetation. In absence of vegetation, the coefficient is set to zero. @@ -290,11 +296,11 @@ def calculate_wind_attenuation_coefficient( canopy_height: Canopy height, [m] leaf_area_index: Total leaf area index, [m m-1] mean_mixing_length: Mixing length for canopy air transport, [m] - drag_coefficient: Drag coefficient - relative_turbulence_intensity: Relative turbulence intensity + drag_coefficient: Drag coefficient, dimensionless + relative_turbulence_intensity: Relative turbulence intensity, dimensionless Returns: - wind attenuation coefficient + wind attenuation coefficient, dimensionless """ intermediate_coefficient = ( @@ -325,7 +331,7 @@ def wind_log_profile( diabatic_correction_momentum: Diabatic correction factor for momentum Returns: - logarithmic wind profile, [m s-1] + logarithmic wind profile """ wind_profile = ( @@ -344,7 +350,7 @@ def calculate_fricition_velocity( diabatic_correction_momentum: Union[float, NDArray[np.float32]], von_karmans_constant: float, ) -> NDArray[np.float32]: - """Calculate friction velocity from wind speed at reference height. + """Calculate friction velocity from wind speed at reference height, [m s-1]. Args: wind_speed_ref: Wind speed at reference height, [m s-1] @@ -355,8 +361,9 @@ def calculate_fricition_velocity( diabatic_correction_momentum: Diabatic correction factor for momentum as returned by :func:`~virtual_rainforest.models.abiotic.wind.calculate_diabatic_correction_above` - von_karmans_constant: Von Karman's constant, unitless, describes the logarithmic - velocity profile of a turbulent fluid near a no-slip boundary + von_karmans_constant: Von Karman's constant, dimensionless constant describing + the logarithmic velocity profile of a turbulent fluid near a no-slip + boundary. Returns: friction velocity @@ -381,7 +388,7 @@ def calculate_wind_above_canopy( von_karmans_constant: float, min_wind_speed_above_canopy: float, ) -> NDArray[np.float32]: - """Calculate wind speed above canopy from wind speed at reference height. + """Calculate wind speed above canopy from wind speed at reference height, [m s-1]. Wind speed above the canopy dictates heat and vapor exchange between the canopy and the air above it, and therefore ultimately determines temperature and vapor @@ -402,8 +409,9 @@ def calculate_wind_above_canopy( diabatic_correction_momentum: Diabatic correction factor for momentum as returned by :func:`~virtual_rainforest.models.abiotic.wind.calculate_diabatic_correction_above` - von_karmans_constant: Von Karman's constant, unitless, describes the logarithmic - velocity profile of a turbulent fluid near a no-slip boundary + von_karmans_constant: Von Karman's constant, dimensionless constant describing + the logarithmic velocity profile of a turbulent fluid near a no-slip + boundary. Returns: wind speed at required heights above canopy, [m s-1] @@ -431,12 +439,12 @@ def calculate_wind_canopy( attenuation_coefficient: NDArray[np.float32], min_windspeed_below_canopy: float, ) -> NDArray[np.float32]: - """Calculate wind profile for individual canopy layers with variable leaf area. + """Calculate wind speed in a multi-layer canopy, [m s-1]. This function can be extended to account for edge distance effects. Args: - top_of_canopy_wind_speed: Wind speed at top of canopy layer, [m /s] + top_of_canopy_wind_speed: Wind speed at top of canopy layer, [m s-1] wind_layer_heights: Height of canopy layer node, [m] canopy_height: Height to top of canopy layer, [m] attenuation_coefficient: Mean attenuation coefficient based on the profile @@ -446,7 +454,7 @@ def calculate_wind_canopy( vegetation, [m/s]. This value is set to avoid dividion by zero. Returns: - wind speed at height of canopy node, [m/s] + wind speed at height of canopy node, [m s-1] """ return np.nan_to_num( @@ -470,7 +478,7 @@ def calculate_wind_profile( abiotic_constants: AbioticConsts, core_constants: CoreConsts, ) -> dict[str, NDArray[np.float32]]: - r"""Calculate wind speed above and below the canoopy, [m s-1]. + r"""Calculate wind speed above and below the canopy, [m s-1]. The wind profile above the canopy is described as follows (based on :cite:p:`campbell_introduction_1998` as implemented in From 67cbfbdfca0812bbf3da4f11232f4c49594d4f91 Mon Sep 17 00:00:00 2001 From: vgro Date: Mon, 4 Dec 2023 14:39:16 +0000 Subject: [PATCH 14/18] updateds on time index and dimesnions --- tests/conftest.py | 17 +++- tests/models/abiotic/test_abiotic_model.py | 85 ++++++++++++++++++- tests/models/abiotic/test_wind.py | 19 ++++- .../models/abiotic/abiotic_model.py | 14 +-- virtual_rainforest/models/abiotic/wind.py | 2 +- 5 files changed, 122 insertions(+), 15 deletions(-) diff --git a/tests/conftest.py b/tests/conftest.py index 35130a930..65278ff65 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -319,6 +319,9 @@ def dummy_climate_data(layer_roles_fixture): np.full((3, 3), 30), dims=["cell_id", "time_index"], ) + data["wind_speed_ref"] = DataArray( + [[0, 5, 10], [0, 5, 10], [0, 5, 10]], dims=["time_index", "cell_id"] + ) data["mean_annual_temperature"] = DataArray( np.full((3), 20), dims=["cell_id"], @@ -339,6 +342,18 @@ def dummy_climate_data(layer_roles_fixture): np.full((3, 3), 400), dims=["cell_id", "time_index"], ) + pressure = np.repeat(a=[96.0, np.nan], repeats=[13, 2]) + data["atmospheric_pressure"] = DataArray( + np.broadcast_to(pressure, (3, 15)).T, + dims=["layers", "cell_id"], + coords={ + "layers": np.arange(15), + "layer_roles": ("layers", layer_roles_fixture), + "cell_id": data.grid.cell_id, + }, + name="atmospheric_pressure", + ) + evapotranspiration = np.repeat(a=[np.nan, 20.0, np.nan], repeats=[1, 3, 11]) data["evapotranspiration"] = DataArray( np.broadcast_to(evapotranspiration, (3, 15)).T, @@ -475,5 +490,5 @@ def dummy_climate_data(layer_roles_fixture): data["canopy_height"] = DataArray([32, 32, 32], dims=["cell_id"]) data["sensible_heat_flux_topofcanopy"] = DataArray([100, 50, 10], dims=["cell_id"]) data["friction_velocity"] = DataArray([12, 5, 2], dims=["cell_id"]) - data["wind_speed_ref"] = DataArray([0, 5, 10], dims=["cell_id"]) + return data diff --git a/tests/models/abiotic/test_abiotic_model.py b/tests/models/abiotic/test_abiotic_model.py index dcee82858..419059946 100644 --- a/tests/models/abiotic/test_abiotic_model.py +++ b/tests/models/abiotic/test_abiotic_model.py @@ -3,8 +3,10 @@ from contextlib import nullcontext as does_not_raise from logging import CRITICAL, DEBUG, ERROR, INFO +import numpy as np import pint import pytest +from xarray import DataArray from tests.conftest import log_check from virtual_rainforest.core.exceptions import ConfigurationError @@ -44,7 +46,7 @@ def test_abiotic_model_initialization( (DEBUG, "abiotic model: required var 'canopy_height' checked"), (DEBUG, "abiotic model: required var 'layer_heights' checked"), (DEBUG, "abiotic model: required var 'leaf_area_index' checked"), - (DEBUG, "abiotic model: required var 'atmospheric_pressure_ref' checked"), + (DEBUG, "abiotic model: required var 'atmospheric_pressure' checked"), ( DEBUG, ( @@ -91,7 +93,7 @@ def test_abiotic_model_initialization_no_data(caplog, dummy_climate_data): ERROR, ( "abiotic model: init data missing required var" - " 'atmospheric_pressure_ref'" + " 'atmospheric_pressure'" ), ), ( @@ -128,7 +130,7 @@ def test_abiotic_model_initialization_no_data(caplog, dummy_climate_data): (DEBUG, "abiotic model: required var 'leaf_area_index' checked"), ( DEBUG, - "abiotic model: required var 'atmospheric_pressure_ref' checked", + "abiotic model: required var 'atmospheric_pressure' checked", ), ( DEBUG, @@ -160,7 +162,7 @@ def test_abiotic_model_initialization_no_data(caplog, dummy_climate_data): (DEBUG, "abiotic model: required var 'leaf_area_index' checked"), ( DEBUG, - "abiotic model: required var 'atmospheric_pressure_ref' checked", + "abiotic model: required var 'atmospheric_pressure' checked", ), ( DEBUG, @@ -278,3 +280,78 @@ def test_generate_abiotic_model_bounds_error( # Final check that expected logging entries are produced log_check(caplog, expected_log_entries) + + +@pytest.mark.parametrize( + "cfg_string,time_interval", + [ + pytest.param( + "[core]\n[core.timing]\nupdate_interval = '1 day'\n[abiotic]\n", + pint.Quantity("1 day"), + id="updates correctly", + ), + ], +) +def test_update_abiotic_model( + dummy_climate_data, layer_roles_fixture, cfg_string, time_interval +): + """Test that update() returns expected output in data object.""" + + from virtual_rainforest.core.config import Config + from virtual_rainforest.core.registry import register_module + from virtual_rainforest.models.abiotic.abiotic_model import AbioticModel + + # Register the module components to access constants classes + register_module("virtual_rainforest.models.abiotic") + + # Build the config object + config = Config(cfg_strings=cfg_string) + + # initialise model + model = AbioticModel.from_config( + dummy_climate_data, + config, + pint.Quantity(config["core"]["timing"]["update_interval"]), + ) + + model.setup() + model.update(time_index=0) + + friction_velocity_exp = np.array( + [ + [0.0, 0.818637, 1.638679], + [0.0, 0.81887, 1.638726], + [0.0, 0.820036, 1.638959], + [0.0, 0.821194, 1.639192], + [0.0, 0.822174, 1.63939], + [0.0, 0.822336, 1.639422], + ] + ) + wind_speed_exp = np.array( + [ + [0.55, 5.536364, 11.07365], + [0.54557, 5.491774, 10.984462], + [0.523951, 5.274152, 10.549181], + [0.503188, 5.065153, 10.13115], + [0.486188, 4.89403, 9.788873], + [0.483444, 4.866404, 9.733618], + ] + ) + + wind_above_exp = np.array([0.55, 5.536364, 11.07365]) + + np.testing.assert_allclose( + model.data["wind_speed_above_canopy"], wind_above_exp, rtol=1e-3, atol=1e-3 + ) + np.testing.assert_allclose( + model.data["friction_velocity"], + DataArray(np.concatenate((friction_velocity_exp, np.full((9, 3), np.nan)))), + rtol=1e-3, + atol=1e-3, + ) + np.testing.assert_allclose( + model.data["wind_speed_canopy"], + DataArray(np.concatenate((wind_speed_exp, np.full((9, 3), np.nan)))), + rtol=1e-3, + atol=1e-3, + ) diff --git a/tests/models/abiotic/test_wind.py b/tests/models/abiotic/test_wind.py index eddf6f98c..d76b447f4 100644 --- a/tests/models/abiotic/test_wind.py +++ b/tests/models/abiotic/test_wind.py @@ -238,7 +238,9 @@ def test_calculate_friction_velocity(dummy_climate_data): from virtual_rainforest.models.abiotic.wind import calculate_fricition_velocity result = calculate_fricition_velocity( - wind_speed_ref=dummy_climate_data.data["wind_speed_ref"], + wind_speed_ref=( + dummy_climate_data.data["wind_speed_ref"].isel(time_index=0).to_numpy() + ), reference_height=(dummy_climate_data.data["canopy_height"] + 10).to_numpy(), zeroplane_displacement=np.array([0.0, 25.312559, 27.58673]), roughness_length_momentum=np.array([0.017, 1.4533, 0.9591]), @@ -329,7 +331,9 @@ def test_calculate_wind_profile(dummy_climate_data): air_temperature=air_temperature.to_numpy(), atmospheric_pressure=np.array([96, 96, 96]), sensible_heat_flux_topofcanopy=np.array([100, 50, 10]), - wind_speed_ref=dummy_climate_data.data["wind_speed_ref"].to_numpy(), + wind_speed_ref=( + dummy_climate_data.data["wind_speed_ref"].isel(time_index=0).to_numpy() + ), wind_reference_height=( dummy_climate_data.data["canopy_height"] + 10 ).to_numpy(), @@ -338,7 +342,16 @@ def test_calculate_wind_profile(dummy_climate_data): core_constants=CoreConsts, ) - friction_velocity_exp = np.array([0.0, 0.818637, 1.638679]) + friction_velocity_exp = np.array( + [ + [0.0, 0.818637, 1.638679], + [0.0, 0.81887, 1.638726], + [0.0, 0.820036, 1.638959], + [0.0, 0.821194, 1.639192], + [0.0, 0.822174, 1.63939], + [0.0, 0.822336, 1.639422], + ] + ) wind_speed_exp = np.array( [ [0.55, 5.536364, 11.07365], diff --git a/virtual_rainforest/models/abiotic/abiotic_model.py b/virtual_rainforest/models/abiotic/abiotic_model.py index 8aa1afb30..ab2d4e6fa 100644 --- a/virtual_rainforest/models/abiotic/abiotic_model.py +++ b/virtual_rainforest/models/abiotic/abiotic_model.py @@ -57,7 +57,7 @@ class AbioticModel(BaseModel): ("canopy_height", ("spatial",)), ("layer_heights", ("spatial",)), ("leaf_area_index", ("spatial",)), - ("atmospheric_pressure_ref", ("spatial",)), + ("atmospheric_pressure", ("spatial",)), ("sensible_heat_flux_topofcanopy", ("spatial",)), ("wind_speed_ref", ("spatial",)), ) @@ -144,7 +144,7 @@ def update(self, time_index: int, **kwargs: Any) -> None: .where(self.data[var].layer_roles != "soil") .dropna(dim="layers") ) - wind_update_inputs[var].append(selection[var]) + wind_update_inputs[var] = selection wind_update = wind.calculate_wind_profile( canopy_height=self.data["canopy_height"].to_numpy(), @@ -152,11 +152,13 @@ def update(self, time_index: int, **kwargs: Any) -> None: wind_layer_heights=wind_update_inputs["layer_heights"].to_numpy(), leaf_area_index=wind_update_inputs["leaf_area_index"].to_numpy(), air_temperature=wind_update_inputs["air_temperature"].to_numpy(), - atmospheric_pressure=self.data["atmospheric_pressure_ref"].to_numpy(), + atmospheric_pressure=self.data["atmospheric_pressure"].to_numpy()[0], sensible_heat_flux_topofcanopy=( self.data["sensible_heat_flux_topofcanopy"].to_numpy() ), - wind_speed_ref=self.data["wind_speed_ref"].to_numpy(), + wind_speed_ref=( + self.data["wind_speed_ref"].isel(time_index=time_index).to_numpy() + ), wind_reference_height=(self.data["canopy_height"] + 10).to_numpy(), turbulence_sign=True, abiotic_constants=AbioticConsts(), @@ -178,7 +180,7 @@ def update(self, time_index: int, **kwargs: Any) -> None: wind_update[var], np.full( ( - self.layer_roles.count("soil"), + len(self.layer_roles) - len(wind_update[var]), self.data.grid.n_cells, ), np.nan, @@ -188,7 +190,7 @@ def update(self, time_index: int, **kwargs: Any) -> None: dims=self.data["layer_heights"].dims, coords=self.data["layer_heights"].coords, ) - wind_output[var] = var_out + wind_output[var] = var_out self.data.add_from_dict(output_dict=wind_output) diff --git a/virtual_rainforest/models/abiotic/wind.py b/virtual_rainforest/models/abiotic/wind.py index 2e5dc9395..50380ef43 100644 --- a/virtual_rainforest/models/abiotic/wind.py +++ b/virtual_rainforest/models/abiotic/wind.py @@ -591,7 +591,7 @@ def calculate_wind_profile( diabatic_correction_momentum=diabatic_correction_above["psi_m"], von_karmans_constant=core_constants.von_karmans_constant, ) - output["friction_velocity"] = friction_velocity[0] + output["friction_velocity"] = friction_velocity mean_mixing_length = calculate_mean_mixing_length( canopy_height=canopy_height, From 62443ffb934940d2b484366554c204dff05f208f Mon Sep 17 00:00:00 2001 From: vgro Date: Mon, 4 Dec 2023 14:46:28 +0000 Subject: [PATCH 15/18] constants values in some tests to make it clearer --- tests/models/abiotic/test_wind.py | 30 +++++++++--------------------- 1 file changed, 9 insertions(+), 21 deletions(-) diff --git a/tests/models/abiotic/test_wind.py b/tests/models/abiotic/test_wind.py index d76b447f4..5a7e863f3 100644 --- a/tests/models/abiotic/test_wind.py +++ b/tests/models/abiotic/test_wind.py @@ -15,7 +15,7 @@ def test_calculate_zero_plane_displacement(dummy_climate_data): result = calculate_zero_plane_displacement( canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(), leaf_area_index=np.array([0, 3, 7]), - zero_plane_scaling_parameter=AbioticConsts.zero_plane_scaling_parameter, + zero_plane_scaling_parameter=7.5, ) np.testing.assert_allclose(result, np.array([0.0, 25.312559, 27.58673])) @@ -32,19 +32,11 @@ def test_calculate_roughness_length_momentum(dummy_climate_data): canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(), leaf_area_index=np.array([0, 3, 7]), zero_plane_displacement=np.array([0.0, 25.312559, 27.58673]), - substrate_surface_drag_coefficient=( - AbioticConsts.substrate_surface_drag_coefficient - ), - roughness_element_drag_coefficient=( - AbioticConsts.roughness_element_drag_coefficient - ), - roughness_sublayer_depth_parameter=( - AbioticConsts.roughness_sublayer_depth_parameter - ), - max_ratio_wind_to_friction_velocity=( - AbioticConsts.max_ratio_wind_to_friction_velocity - ), - min_roughness_length=AbioticConsts.min_roughness_length, + substrate_surface_drag_coefficient=0.003, + roughness_element_drag_coefficient=0.3, + roughness_sublayer_depth_parameter=0.193, + max_ratio_wind_to_friction_velocity=0.3, + min_roughness_length=0.05, von_karman_constant=CoreConsts.von_karmans_constant, ) @@ -138,12 +130,8 @@ def test_generate_relative_turbulence_intensity(dummy_climate_data): ) result = generate_relative_turbulence_intensity( layer_heights=layer_heights.to_numpy(), - min_relative_turbulence_intensity=( - AbioticConsts.min_relative_turbulence_intensity - ), - max_relative_turbulence_intensity=( - AbioticConsts.max_relative_turbulence_intensity - ), + min_relative_turbulence_intensity=0.36, + max_relative_turbulence_intensity=0.9, increasing_with_height=True, ) @@ -263,7 +251,7 @@ def test_calculate_wind_above_canopy(dummy_climate_data): roughness_length_momentum=np.array([0.017, 1.4533, 0.9591]), diabatic_correction_momentum=np.array([0.003, 0.026, 0.013]), von_karmans_constant=CoreConsts.von_karmans_constant, - min_wind_speed_above_canopy=AbioticConsts.min_wind_speed_above_canopy, + min_wind_speed_above_canopy=0.55, ) exp_result = np.array([0.55, 5.590124, 10.750233]) From f736c7d82e5b179d26740c525e8292c9bf639986 Mon Sep 17 00:00:00 2001 From: vgro Date: Tue, 5 Dec 2023 15:21:10 +0000 Subject: [PATCH 16/18] Jacob's suggestions implemented --- tests/models/abiotic/test_wind.py | 25 ++++++++++++++++--- .../models/abiotic/abiotic_model.py | 1 - .../models/abiotic/constants.py | 3 +++ virtual_rainforest/models/abiotic/wind.py | 17 ++++++++++--- 4 files changed, 38 insertions(+), 8 deletions(-) diff --git a/tests/models/abiotic/test_wind.py b/tests/models/abiotic/test_wind.py index 5a7e863f3..d39d87eac 100644 --- a/tests/models/abiotic/test_wind.py +++ b/tests/models/abiotic/test_wind.py @@ -128,14 +128,14 @@ def test_generate_relative_turbulence_intensity(dummy_climate_data): .where(dummy_climate_data.data["layer_heights"].layer_roles != "soil") .dropna(dim="layers") ) - result = generate_relative_turbulence_intensity( + result_T = generate_relative_turbulence_intensity( layer_heights=layer_heights.to_numpy(), min_relative_turbulence_intensity=0.36, max_relative_turbulence_intensity=0.9, increasing_with_height=True, ) - exp_result = np.array( + exp_result_T = np.array( [ [17.64, 17.64, 17.64], [16.56, 16.56, 16.56], @@ -145,7 +145,25 @@ def test_generate_relative_turbulence_intensity(dummy_climate_data): [0.414, 0.414, 0.414], ] ) - np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3) + result_F = generate_relative_turbulence_intensity( + layer_heights=layer_heights.to_numpy(), + min_relative_turbulence_intensity=0.36, + max_relative_turbulence_intensity=0.9, + increasing_with_height=False, + ) + + exp_result_F = np.array( + [ + [-16.38, -16.38, -16.38], + [-15.3, -15.3, -15.3], + [-9.9, -9.9, -9.9], + [-4.5 , -4.5, -4.5], + [0.09, 0.09,0.09], + [0.846, 0.846, 0.846], + ] + ) + np.testing.assert_allclose(result_T, exp_result_T, rtol=1e-3, atol=1e-3) + np.testing.assert_allclose(result_F, exp_result_F, rtol=1e-3, atol=1e-3) def test_calculate_wind_attenuation_coefficient(dummy_climate_data): @@ -325,7 +343,6 @@ def test_calculate_wind_profile(dummy_climate_data): wind_reference_height=( dummy_climate_data.data["canopy_height"] + 10 ).to_numpy(), - turbulence_sign=True, abiotic_constants=AbioticConsts, core_constants=CoreConsts, ) diff --git a/virtual_rainforest/models/abiotic/abiotic_model.py b/virtual_rainforest/models/abiotic/abiotic_model.py index ab2d4e6fa..2ed96eea7 100644 --- a/virtual_rainforest/models/abiotic/abiotic_model.py +++ b/virtual_rainforest/models/abiotic/abiotic_model.py @@ -160,7 +160,6 @@ def update(self, time_index: int, **kwargs: Any) -> None: self.data["wind_speed_ref"].isel(time_index=time_index).to_numpy() ), wind_reference_height=(self.data["canopy_height"] + 10).to_numpy(), - turbulence_sign=True, abiotic_constants=AbioticConsts(), core_constants=CoreConsts(), ) diff --git a/virtual_rainforest/models/abiotic/constants.py b/virtual_rainforest/models/abiotic/constants.py index 9c15529d0..a06bc3ce8 100644 --- a/virtual_rainforest/models/abiotic/constants.py +++ b/virtual_rainforest/models/abiotic/constants.py @@ -160,3 +160,6 @@ class AbioticConsts(ConstantsDataclass): Dimenionless parameter, implementation after :cite:t:`maclean_microclimc_2021` and values taken from :cite:t:`yasuda_turbulent_1988`.""" + + turbulence_sign: bool = True + """Flag indicating if turbulence increases or decreases with height.""" diff --git a/virtual_rainforest/models/abiotic/wind.py b/virtual_rainforest/models/abiotic/wind.py index 50380ef43..42106f8ff 100644 --- a/virtual_rainforest/models/abiotic/wind.py +++ b/virtual_rainforest/models/abiotic/wind.py @@ -445,7 +445,7 @@ def calculate_wind_canopy( Args: top_of_canopy_wind_speed: Wind speed at top of canopy layer, [m s-1] - wind_layer_heights: Height of canopy layer node, [m] + wind_layer_heights: Heights of canopy layer nodes, [m] canopy_height: Height to top of canopy layer, [m] attenuation_coefficient: Mean attenuation coefficient based on the profile calculated by @@ -474,7 +474,6 @@ def calculate_wind_profile( sensible_heat_flux_topofcanopy: NDArray[np.float32], wind_speed_ref: NDArray[np.float32], wind_reference_height: Union[float, NDArray[np.float32]], - turbulence_sign: bool, abiotic_constants: AbioticConsts, core_constants: CoreConsts, ) -> dict[str, NDArray[np.float32]]: @@ -505,6 +504,18 @@ def calculate_wind_profile( details, see :cite:t:`maclean_microclimc_2021`. Args: + canopy_height: Canopy height, [m] + wind_height_above: Height above canopy for which wind speed is required, [m] + wind_layer_heights: Heights of canopy layer nodes, [m] + leaf_area_index: Leaf area index, [m m-1] + air_temperature: Air temperature, [C] + atmospheric_pressure: Atmospheric pressure, [kPa] + sensible_heat_flux_topofcanopy: Sensible heat flux from the top of the canopy to + the atmosphere, [W m-2], + wind_speed_ref: Wind speed at reference height, [m s-1] + wind_reference_height: Reference height for wind measurement, [m] + abiotic_constants: Specific constants for the abiotic model + core_constants: Universal constants shared across all models Returns: dictionnary that contains wind speed above the canopy, [m s-1], wind speed @@ -608,7 +619,7 @@ def calculate_wind_profile( max_relative_turbulence_intensity=( abiotic_constants.max_relative_turbulence_intensity ), - increasing_with_height=turbulence_sign, + increasing_with_height=AbioticConsts.turbulence_sign, ) attennuation_coefficient = calculate_wind_attenuation_coefficient( From 02a70686f1ad0f466033344f0460fc2ccec98af7 Mon Sep 17 00:00:00 2001 From: vgro Date: Tue, 5 Dec 2023 15:21:41 +0000 Subject: [PATCH 17/18] flake8 fix --- tests/models/abiotic/test_wind.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/models/abiotic/test_wind.py b/tests/models/abiotic/test_wind.py index d39d87eac..e356f02f9 100644 --- a/tests/models/abiotic/test_wind.py +++ b/tests/models/abiotic/test_wind.py @@ -157,8 +157,8 @@ def test_generate_relative_turbulence_intensity(dummy_climate_data): [-16.38, -16.38, -16.38], [-15.3, -15.3, -15.3], [-9.9, -9.9, -9.9], - [-4.5 , -4.5, -4.5], - [0.09, 0.09,0.09], + [-4.5, -4.5, -4.5], + [0.09, 0.09, 0.09], [0.846, 0.846, 0.846], ] ) From bd9f8d77bfa62a69665c25499c4a03f88f034de8 Mon Sep 17 00:00:00 2001 From: vgro Date: Tue, 5 Dec 2023 15:52:13 +0000 Subject: [PATCH 18/18] use of constants in AbioticModel fixed --- tests/models/abiotic/test_abiotic_model.py | 7 +++++++ virtual_rainforest/models/abiotic/abiotic_model.py | 11 ++++++++--- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/tests/models/abiotic/test_abiotic_model.py b/tests/models/abiotic/test_abiotic_model.py index 419059946..3d3bf21e0 100644 --- a/tests/models/abiotic/test_abiotic_model.py +++ b/tests/models/abiotic/test_abiotic_model.py @@ -19,6 +19,7 @@ def test_abiotic_model_initialization( ): """Test `AbioticModel` initialization.""" from virtual_rainforest.core.base_model import BaseModel + from virtual_rainforest.core.constants import CoreConsts from virtual_rainforest.models.abiotic.abiotic_model import AbioticModel from virtual_rainforest.models.abiotic.constants import AbioticConsts @@ -29,6 +30,7 @@ def test_abiotic_model_initialization( soil_layers=[-0.5, -1.0], canopy_layers=10, constants=AbioticConsts(), + core_constants=CoreConsts(), ) # In cases where it passes then checks that the object has the right properties @@ -62,6 +64,7 @@ def test_abiotic_model_initialization( def test_abiotic_model_initialization_no_data(caplog, dummy_climate_data): """Test `AbioticModel` initialization with no data.""" + from virtual_rainforest.core.constants import CoreConsts from virtual_rainforest.core.data import Data from virtual_rainforest.core.grid import Grid from virtual_rainforest.models.abiotic.abiotic_model import AbioticModel @@ -79,6 +82,7 @@ def test_abiotic_model_initialization_no_data(caplog, dummy_climate_data): soil_layers=[-0.5, -1.0], canopy_layers=10, constants=AbioticConsts, + core_constants=CoreConsts, ) # Final check that expected logging entries are produced @@ -119,6 +123,7 @@ def test_abiotic_model_initialization_no_data(caplog, dummy_climate_data): does_not_raise(), ( (INFO, "Initialised abiotic.AbioticConsts from config"), + (INFO, "Initialised core.CoreConsts from config"), ( INFO, "Information required to initialise the abiotic model successfully " @@ -151,6 +156,7 @@ def test_abiotic_model_initialization_no_data(caplog, dummy_climate_data): does_not_raise(), ( (INFO, "Initialised abiotic.AbioticConsts from config"), + (INFO, "Initialised core.CoreConsts from config"), ( INFO, "Information required to initialise the abiotic model successfully " @@ -235,6 +241,7 @@ def test_generate_abiotic_model( pytest.raises(ConfigurationError), ( (INFO, "Initialised abiotic.AbioticConsts from config"), + (INFO, "Initialised core.CoreConsts from config"), ( INFO, "Information required to initialise the abiotic model " diff --git a/virtual_rainforest/models/abiotic/abiotic_model.py b/virtual_rainforest/models/abiotic/abiotic_model.py index 2ed96eea7..8e34dfbc7 100644 --- a/virtual_rainforest/models/abiotic/abiotic_model.py +++ b/virtual_rainforest/models/abiotic/abiotic_model.py @@ -72,6 +72,7 @@ def __init__( soil_layers: list[float], canopy_layers: int, constants: AbioticConsts, + core_constants: CoreConsts, **kwargs: Any, ): super().__init__(data, update_interval, **kwargs) @@ -86,7 +87,9 @@ def __init__( self.update_interval """The time interval between model updates.""" self.constants = constants - """Set of constants for the abiotic model""" + """Set of constants for the abiotic model.""" + self.core_constants = core_constants + """Set of universal constants that are used across all models.""" @classmethod def from_config( @@ -110,6 +113,7 @@ def from_config( # Load in the relevant constants constants = load_constants(config, "abiotic", "AbioticConsts") + core_constants = load_constants(config, "core", "CoreConsts") LOGGER.info( "Information required to initialise the abiotic model successfully " @@ -121,6 +125,7 @@ def from_config( soil_layers, canopy_layers, constants, + core_constants, ) def setup(self) -> None: @@ -160,8 +165,8 @@ def update(self, time_index: int, **kwargs: Any) -> None: self.data["wind_speed_ref"].isel(time_index=time_index).to_numpy() ), wind_reference_height=(self.data["canopy_height"] + 10).to_numpy(), - abiotic_constants=AbioticConsts(), - core_constants=CoreConsts(), + abiotic_constants=self.constants, + core_constants=self.core_constants, ) wind_output = {}