From 1217524dee1622d0338d6f9d5aeafa0b612699a4 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Tue, 26 Nov 2024 01:47:17 +0100 Subject: [PATCH] new examples: R.R. Rogers 1975 (Atmosphere 13(4)) fig 1 (parcel model using SciPy+Pint, without PySDM) --- .../PySDM_examples/Rogers_1975/__init__.py | 9 + .../PySDM_examples/Rogers_1975/fig_1.ipynb | 1566 +++++++++++++++++ 2 files changed, 1575 insertions(+) create mode 100644 examples/PySDM_examples/Rogers_1975/__init__.py create mode 100644 examples/PySDM_examples/Rogers_1975/fig_1.ipynb diff --git a/examples/PySDM_examples/Rogers_1975/__init__.py b/examples/PySDM_examples/Rogers_1975/__init__.py new file mode 100644 index 000000000..0bc8eb68a --- /dev/null +++ b/examples/PySDM_examples/Rogers_1975/__init__.py @@ -0,0 +1,9 @@ +""" +based on Rogers 1975 (Atmosphere) +https://doi.org/10.1080/00046973.1975.9648397 + +fig_1.ipynb: +.. include:: ./fig_1.ipynb.badges.md +""" + +# pylint: disable=invalid-name diff --git a/examples/PySDM_examples/Rogers_1975/fig_1.ipynb b/examples/PySDM_examples/Rogers_1975/fig_1.ipynb new file mode 100644 index 000000000..435580565 --- /dev/null +++ b/examples/PySDM_examples/Rogers_1975/fig_1.ipynb @@ -0,0 +1,1566 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "bfd1a1e7-29d9-443d-9eb4-ef75151ffcc7", + "metadata": {}, + "source": [ + "[![preview notebook](https://img.shields.io/static/v1?label=render%20on&logo=github&color=87ce3e&message=GitHub)](https://github.com/open-atmos/PySDM/blob/main/examples/PySDM_examples/Rogers_1975/fig_1.ipynb) \n", + "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/open-atmos/PySDM.git/main?urlpath=examples/PySDM_examples/Rogers_1975/fig_1.ipynb) \n", + "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/open-atmos/PySDM/blob/main/examples/PySDM_examples/Rogers_1975/fig_1.ipynb)" + ] + }, + { + "cell_type": "markdown", + "id": "ca764a8c-a315-4dd3-9814-bd59401f6ccd", + "metadata": {}, + "source": [ + "basic parcel simulation based on [Rogers 1975](https://doi.org/10.1080/00046973.1975.9648397) implemented using SciPy+Pint but without PySDM - for illustrative purposes." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "99b453fa", + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "if 'google.colab' in sys.modules:\n", + " !pip --quiet install open-atmos-jupyter-utils\n", + " from open_atmos_jupyter_utils import pip_install_on_colab\n", + " pip_install_on_colab('PySDM-examples')" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "7dadd922", + "metadata": { + "ExecuteTime": { + "end_time": "2023-07-07T15:38:26.856395076Z", + "start_time": "2023-07-07T15:38:22.157465793Z" + } + }, + "outputs": [], + "source": [ + "import collections\n", + "import functools\n", + "from open_atmos_jupyter_utils import show_plot\n", + "from matplotlib import pyplot\n", + "import numpy as np\n", + "import scipy\n", + "from scipy import constants\n", + "import chempy" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "0ecc07d1-88e4-4a4f-a469-fe4f14b31ac4", + "metadata": {}, + "outputs": [], + "source": [ + "import pint\n", + "si = pint.UnitRegistry()\n", + "si.setup_matplotlib()\n", + "si.define('fraction = [] = frac')\n", + "si.define('percent = 1e-2 frac = pct')" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "31cdfaaa-31be-46c6-9dcf-257024240cd6", + "metadata": {}, + "outputs": [], + "source": [ + "class Storage:\n", + " \"\"\" state vector representation with each element having its own Pint-compatible\n", + " physical dimension, thus allowing to seamlessly couple Pint and scipy.odeint\n", + " (all methods return objects that inherit from `numpy.ndarray` but are additionally\n", + " equipped with .VAR unit-aware setters and getters, allowing both unit-unaware\n", + " whole-array expressions (e.g., `state += dt * deriv`) as well as unit-aware \n", + " operations on state vars (e.g., `state.T = 300 * si.K` or `state.m[:] = ...`) \"\"\"\n", + "\n", + " var_units = {\n", + " 'p': si.Pa,\n", + " 'T': si.K,\n", + " 'S': si.dimensionless,\n", + " 'r': si.metre,\n", + " }\n", + " \n", + " der_unit = si.second\n", + "\n", + " @staticmethod\n", + " def __make_storage(shape, deriv=False):\n", + " def getter(self, idx, unit):\n", + " return self[idx] * unit\n", + "\n", + " def setter(self, value, idx, unit):\n", + " self[idx] = (value.to(unit) / unit).magnitude\n", + "\n", + " properties = {}\n", + " for i, key in enumerate(Storage.var_units.keys()):\n", + " kwargs = {\n", + " 'unit': Storage.var_units[key] / (Storage.der_unit if deriv else 1),\n", + " 'idx': i\n", + " }\n", + " properties[key] = property(\n", + " functools.partial(getter, **kwargs),\n", + " functools.partial(setter, **kwargs),\n", + " ) \n", + " \n", + " return type(\"StorageImpl\", (np.ndarray,), properties)(shape)\n", + "\n", + " @staticmethod\n", + " def make_state():\n", + " \"\"\" returns a newly allocated unit-aware storage \"\"\"\n", + " return Storage.__make_storage((len(Storage.var_units),))\n", + "\n", + " @staticmethod\n", + " def make_deriv(state):\n", + " \"\"\" returns a newly allocated unit-aware storage with size of `state` and derivative dimensions \"\"\"\n", + " return Storage.__make_storage(state.shape, deriv=True)\n", + "\n", + " @staticmethod\n", + " def view_state(array):\n", + " \"\"\" returns a newly allocated unit-aware storage with size and data from unit-unaware `array` \"\"\"\n", + " storage = Storage.__make_storage(array.shape)\n", + " storage[:] = array[:]\n", + " return storage" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "69a2d947-8581-4382-9939-07088584933b", + "metadata": {}, + "outputs": [], + "source": [ + "# constants (coefficients with modified 10 exponents to match SI units)\n", + "C = {\n", + " # gas constant for dry air\n", + " \"R\": scipy.constants.R * si.joule / si.kelvin / si.mole / (\n", + " 0.76 * chempy.Substance.from_formula(\"N2\").mass * si.gram / si.mole\n", + " + 0.23 * chempy.Substance.from_formula(\"O2\").mass * si.gram / si.mole\n", + " + 0.01 * chempy.Substance.from_formula(\"Ar\").mass * si.gram / si.mole\n", + " ),\n", + " # standard acceleration of gravity\n", + " \"g\": scipy.constants.g * si.metre / si.second**2,\n", + " # latent heat of condensation\n", + " \"L\": 2.5e3 * si.J / si.g,\n", + " # specific gravity of vapour vs. dry air\n", + " \"eps\": .622,\n", + " # specific heat at constant pressure of dry air\n", + " \"cp\": 1005 * si.joule / si.kilogram / si.kelvin,\n", + " # density of liquid water\n", + " \"rho_L\": 1 * si.kg / si.litre,\n", + " # eq A.1\n", + " \"K\": lambda T: 2.42e-2 * (393/(T.to(si.K).magnitude + 120)) * (T.to(si.K).magnitude/273)**(3/2) * si.J / si.m / si.s / si.K,\n", + " # eq A.2\n", + " \"D_over_K\": lambda p, T: 8.28 * T.to(si.K).magnitude / p.to(si.dyne / si.cm**2).magnitude * si.m**3 / si.J * si.K,\n", + " # eq A.3\n", + " \"e_s\": lambda T: 2.75e12 * np.exp(-5.44e3 / (T/si.K)) * si.ubar,\n", + "}\n", + "C = collections.namedtuple(\"C\", C.keys())(**C)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "7d9b98c0-971c-41a0-87c2-f7adcfc958a1", + "metadata": {}, + "outputs": [], + "source": [ + "np.testing.assert_approx_equal(actual=C.e_s(273.15 * si.K).to(si.mbar).magnitude, desired=6.1, significant=2)\n", + "np.testing.assert_approx_equal(actual=C.K(273.15 * si.K).to(si.J / si.m / si.s / si.K).magnitude, desired=2.4e-2, significant=2)\n", + "np.testing.assert_approx_equal(\n", + " actual=(C.K(169.75 * si.K) * C.D_over_K(1000 * si.hPa, 169.75 * si.K)).to(si.m**2 / si.s).magnitude,\n", + " desired=2.26e-5,\n", + " significant=2\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "d2011171-3291-4c0e-ad4e-23ba8bf25f52", + "metadata": {}, + "outputs": [], + "source": [ + "# ODE system\n", + "class System:\n", + " def __init__(self, *, U, nu_0):\n", + " self.U = U\n", + " self.nu_0 = nu_0\n", + "\n", + " def __call__(self, _, state):\n", + " state = Storage.view_state(state)\n", + " deriv = Storage.make_deriv(state)\n", + "\n", + " # eq (8)\n", + " rho = state.p / C.R / state.T\n", + " \n", + " # eq (5) multiplied by Δt\n", + " deriv.p = -rho * C.g * self.U \n", + "\n", + " # eq (2)\n", + " K = C.K(state.T)\n", + " Fk = C.L**2 * C.eps * C.rho_L / K / C.R / state.T**2\n", + " Fd = C.R * state.T * C.rho_L / C.eps / C.D_over_K(state.p, state.T) / K / C.e_s(state.T)\n", + " sigma = (state.S - 1) / (Fk + Fd)\n", + " \n", + " # eq (1)\n", + " deriv.r = sigma / state.r\n", + " \n", + " # time derivative of eq (4)\n", + " dksi_dt = 4 * np.pi * C.rho_L * self.nu_0 * state.r**2 * deriv.r\n", + " \n", + " # eq (6) divided by dt, multiplied by T\n", + " deriv.T = (\n", + " state.T * C.R / C.cp * deriv.p / state.p + \n", + " C.L / C.cp * dksi_dt\n", + " )\n", + "\n", + " # eq (10)\n", + " Q1 = C.L * C.g * C.eps / C.R / C.cp / state.T**2 - C.g / C.R / state.T\n", + " Q2 = C.R * state.T / C.eps / C.e_s(state.T) + C.eps * C.L**2 / C.cp / state.T / state.p\n", + " deriv.S = Q1 * self.U - rho * Q2 * dksi_dt\n", + " \n", + " return deriv" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "0dc34252-4ac9-4d23-89f5-793239183eb7", + "metadata": {}, + "outputs": [], + "source": [ + "def solve(system, state, t_max, max_step):\n", + " integ = scipy.integrate.solve_ivp(\n", + " system,\n", + " [0, t_max / Storage.der_unit],\n", + " state,\n", + " max_step=max_step.to(Storage.der_unit).magnitude,\n", + " method='LSODA'\n", + " )\n", + " assert integ.success, integ.message\n", + " return Storage.view_state(integ.y), integ.t * Storage.der_unit" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "ceff8ab7-8b6d-46a1-83e8-83b2fc893380", + "metadata": {}, + "outputs": [], + "source": [ + "state = Storage.make_state()\n", + "state.p = 800 * si.mbar\n", + "state.T = (273.15 + 7) * si.K\n", + "state.S = 1 * si.dimensionless\n", + "state.r = 8 * si.um\n", + "\n", + "rho0 = state.p / C.R / state.T\n", + "system = System(U=10 * si.m / si.s, nu_0=200 / si.cm**3 / rho0)\n", + "solution, tsteps = solve(system, state, t_max=20 * si.s, max_step=.5 * si.s)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "f00aa5f8-219a-41dd-98fd-c8a411a502b9", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-11-26T01:45:01.286777\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.8.1, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "bd3178a6a902418d8219e28273c4ca7a", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(HTML(value=\"./fig_1.pdf
\"), HTML(value=\"