|
| 1 | +# pylint: disable=missing-module-docstring |
| 2 | +import numpy as np |
| 3 | +import pytest |
| 4 | +from matplotlib import pyplot |
| 5 | +from PySDM import Formulae |
| 6 | +from PySDM import physics |
| 7 | +from PySDM.physics.dimensional_analysis import DimensionalAnalysis |
| 8 | + |
| 9 | + |
| 10 | +class TestIsotopeKineticFractionationFactors: |
| 11 | + @staticmethod |
| 12 | + @pytest.mark.parametrize( |
| 13 | + "variant, kwargs", |
| 14 | + ( |
| 15 | + ( |
| 16 | + physics.isotope_kinetic_fractionation_factors.CraigGordon, |
| 17 | + { |
| 18 | + "turbulence_parameter_n": 1, |
| 19 | + "delta_diff": 1, |
| 20 | + "theta": 1, |
| 21 | + }, |
| 22 | + ), |
| 23 | + ( |
| 24 | + physics.isotope_kinetic_fractionation_factors.JouzelAndMerlivat1984, |
| 25 | + { |
| 26 | + "alpha_equilibrium": 1, |
| 27 | + "heavy_to_light_diffusivity_ratio": 1, |
| 28 | + }, |
| 29 | + ), |
| 30 | + ), |
| 31 | + ) |
| 32 | + def test_units(variant, kwargs): |
| 33 | + """checks that alphas are dimensionless""" |
| 34 | + with DimensionalAnalysis(): |
| 35 | + # arrange |
| 36 | + sut = variant.alpha_kinetic |
| 37 | + |
| 38 | + # act |
| 39 | + result = sut(relative_humidity=1 * physics.si.dimensionless, **kwargs) |
| 40 | + |
| 41 | + # assert |
| 42 | + assert result.check("[]") |
| 43 | + |
| 44 | + @staticmethod |
| 45 | + def test_fig_9_from_jouzel_and_merlivat_1984(plot=False): |
| 46 | + """[Jouzel & Merlivat 1984](https://doi.org/10.1029/JD089iD07p11749)""" |
| 47 | + # arrange |
| 48 | + formulae = Formulae( |
| 49 | + isotope_kinetic_fractionation_factors="JouzelAndMerlivat1984", |
| 50 | + isotope_equilibrium_fractionation_factors="Majoube1970", |
| 51 | + isotope_diffusivity_ratios="Stewart1975", |
| 52 | + ) |
| 53 | + temperatures = formulae.trivia.C2K(np.asarray([-30, -20, -10])) |
| 54 | + saturation = np.linspace(start=1, stop=1.35) |
| 55 | + alpha_s = formulae.isotope_equilibrium_fractionation_factors.alpha_i_18O |
| 56 | + heavy_to_light_diffusivity_ratio = formulae.isotope_diffusivity_ratios.ratio_18O |
| 57 | + sut = formulae.isotope_kinetic_fractionation_factors.alpha_kinetic |
| 58 | + |
| 59 | + # act |
| 60 | + alpha_s = {temperature: alpha_s(temperature) for temperature in temperatures} |
| 61 | + alpha_k = { |
| 62 | + temperature: sut( |
| 63 | + alpha_equilibrium=alpha_s[temperature], |
| 64 | + relative_humidity=saturation, |
| 65 | + heavy_to_light_diffusivity_ratio=heavy_to_light_diffusivity_ratio( |
| 66 | + temperature |
| 67 | + ), |
| 68 | + ) |
| 69 | + for temperature in temperatures |
| 70 | + } |
| 71 | + alpha_s_times_alpha_k = { |
| 72 | + f"{formulae.trivia.K2C(temperature):.3g}C": alpha_k[temperature] |
| 73 | + * alpha_s[temperature] |
| 74 | + for temperature in temperatures |
| 75 | + } |
| 76 | + |
| 77 | + # plot |
| 78 | + pyplot.xlim(saturation[0], saturation[-1]) |
| 79 | + pyplot.ylim(1.003, 1.022) |
| 80 | + pyplot.xlabel("S") |
| 81 | + pyplot.ylabel("alpha_k * alpha_s") |
| 82 | + for k, v in alpha_s_times_alpha_k.items(): |
| 83 | + pyplot.plot(saturation, v, label=k) |
| 84 | + pyplot.legend() |
| 85 | + if plot: |
| 86 | + pyplot.show() |
| 87 | + else: |
| 88 | + pyplot.clf() |
| 89 | + |
| 90 | + # assert |
| 91 | + assert (alpha_s_times_alpha_k["-30C"] > alpha_s_times_alpha_k["-20C"]).all() |
| 92 | + assert (alpha_s_times_alpha_k["-20C"] > alpha_s_times_alpha_k["-10C"]).all() |
| 93 | + for alpha_alpha in alpha_s_times_alpha_k.values(): |
| 94 | + assert (np.diff(alpha_alpha) < 0).all() |
| 95 | + np.testing.assert_approx_equal( |
| 96 | + actual=alpha_s_times_alpha_k["-30C"][0], |
| 97 | + desired=1.021, |
| 98 | + significant=4, |
| 99 | + ) |
| 100 | + np.testing.assert_approx_equal( |
| 101 | + actual=alpha_s_times_alpha_k["-30C"][-1], |
| 102 | + desired=1.0075, |
| 103 | + significant=4, |
| 104 | + ) |
| 105 | + np.testing.assert_approx_equal( |
| 106 | + actual=alpha_s_times_alpha_k["-10C"][0], |
| 107 | + desired=1.0174, |
| 108 | + significant=4, |
| 109 | + ) |
| 110 | + np.testing.assert_approx_equal( |
| 111 | + actual=alpha_s_times_alpha_k["-10C"][-1], |
| 112 | + desired=1.004, |
| 113 | + significant=4, |
| 114 | + ) |
0 commit comments