diff --git a/setup.py b/setup.py index fdf7c4aabb..e7e1781ab5 100644 --- a/setup.py +++ b/setup.py @@ -30,7 +30,7 @@ def get_long_description(): "numpy" + ("==1.24.4" if CI else ""), "Pint" + ("==0.21.1" if CI else ""), "chempy" + ("==0.8.3" if CI else ""), - "scipy" + ("==1.10.1" if CI and not _32bit else ""), + "scipy" + ("==1.12.0" if CI and not _32bit else ">=1.12.0"), "pyevtk" + ("==1.2.0" if CI else ""), ], extras_require={ diff --git a/tutorials/parcel/parcel_ode.ipynb b/tutorials/parcel/parcel_ode.ipynb new file mode 100644 index 0000000000..07c6a39800 --- /dev/null +++ b/tutorials/parcel/parcel_ode.ipynb @@ -0,0 +1,453 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "initial_id", + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-04T20:06:25.738448Z", + "start_time": "2024-04-04T20:06:24.822258Z" + } + }, + "outputs": [], + "source": [ + "\"\"\" minimalistic adiabatic parcel model using Numpy, pint, SciPy/odeint and matplotlib \"\"\"\n", + "import math\n", + "import pint\n", + "import chempy\n", + "import scipy\n", + "import ipywidgets\n", + "import IPython\n", + "import numpy as np\n", + "from matplotlib import pyplot\n", + "si = pint.UnitRegistry()\n", + "si.setup_matplotlib()" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "80db779823de1332", + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-04T20:06:25.743429Z", + "start_time": "2024-04-04T20:06:25.739647Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "outputs": [], + "source": [ + "IDX_S = -1\n", + "IDX_T = -2" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "613863dd840e1a6e", + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-04T20:06:25.762740Z", + "start_time": "2024-04-04T20:06:25.749899Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "outputs": [], + "source": [ + "# pylint: disable=missing-module-docstring,missing-function-docstring,disable-msg=too-many-locals\n", + "def sample_spectrum(parameters):\n", + " spectrum = scipy.stats.lognorm(*(\n", + " math.log(parameters['lognormal_geometric_stdev']),\n", + " 0,\n", + " parameters['lognormal_median_radius'].to(si.m).magnitude\n", + " ))\n", + " n_sd = parameters['n_sd']\n", + " rd_magnitude = [spectrum.ppf((2 * i - 1) / (2 * n_sd)) for i in range(1, n_sd + 1)]\n", + " rd_unit = si.m\n", + " sampled_values = {\n", + " 'multiplicity': [parameters['total_aerosol_number'] / n_sd] * n_sd,\n", + " 'dry radius': [rd * rd_unit for rd in rd_magnitude]\n", + " }\n", + " pyplot.plot(rd_magnitude * si.m,spectrum.pdf(rd_magnitude) / rd_unit)\n", + " pyplot.step(\n", + " (rd_magnitude[:-1] + np.diff(rd_magnitude)/2) * rd_unit,\n", + " 1 / n_sd / np.diff(rd_magnitude) / rd_unit ,\n", + " where='mid',\n", + " label='sampled distribution'\n", + " )\n", + " pyplot.gca().xaxis.set_units(si.um)\n", + " pyplot.gca().yaxis.set_units(si.um**-1)\n", + " pyplot.legend()\n", + "\n", + " pyplot.grid(True)\n", + " pyplot.show()\n", + " return sampled_values\n", + "\n", + "CONST = {\n", + " 'M_dry': (\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", + " 'M_vap': chempy.Substance.from_formula(\"H2O\").mass * si.gram / si.mole,\n", + " 'g': scipy.constants.g * si.m / si.s**2,\n", + " 'water_density': 1000 * si.kg / si.m**3,\n", + " \"specific_heat_constant_pressure\": 1005 * si.J / si.kg / si.K,\n", + " \"latent_heat\": 2.27e6 * si.J / si.kg,\n", + " \"surface_tension_water\": 7.5e-2 * si.N / si.m,\n", + "}\n", + "CONST['R_dry'] = scipy.constants.R * si.J / si.mole / si.K / CONST['M_dry']\n", + "CONST['eps'] = CONST['M_vap'] / CONST['M_dry']\n", + "CONST['R_vap'] = CONST['R_dry'] / CONST['eps']\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "ce741be9fda22c8", + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-04T20:06:26.107804Z", + "start_time": "2024-04-04T20:06:26.105196Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "outputs": [], + "source": [ + "def saturated_vapor_pressure(temperature):\n", + " \"\"\" eq. 2.17 in Rogers & Yau \"\"\"\n", + " return 6.112 * si.hPa * math.exp(\n", + " 17.67 * temperature / (temperature + 243.5 * si.K)\n", + " )\n", + "\n", + "def equil_saturation(radius, kappa, dry_radius, temperature):\n", + " \"\"\" eq. 6.6 in R&Y \"\"\"\n", + " drop_surface_term = (2 * CONST['surface_tension_water']) / (\n", + " CONST['R_dry'] * CONST['water_density'] * temperature)\n", + " drop_solute_term = kappa * dry_radius**3\n", + " return 1 + drop_surface_term/radius - drop_solute_term/radius**3\n", + "\n", + "def funct(_, state, parameters):\n", + " \"\"\" returns time derivative of the state vector \"\"\"\n", + " saturation = state[IDX_S]\n", + " temperature = state[IDX_T] * parameters['units'][IDX_T]\n", + " condensation_rate_constant = (CONST['water_density'] *\n", + " 4 * math.pi /\n", + " parameters['total_aerosol_number'] /\n", + " parameters['mass_of_dry_air'])\n", + "\n", + " water_vapor_diffusivity = 1e-5 * si.m**2 / si.s * (\n", + " (0.015 / si.K * temperature) - 1.9)\n", + " thermal_conductivity = (\n", + " (1.5e-11 * si.W / si.m / si.K**4) * temperature**3\n", + " + (-4.8e-8 * si.W / si.m / si.K**3) * temperature**2\n", + " + (1e-4 * si.W / si.m / si.K**2) * temperature\n", + " + -3.9e-4 * si.W / si.m / si.K\n", + " )\n", + "\n", + " p_vs = saturated_vapor_pressure(temperature)\n", + " vapour_pressure = saturation * p_vs\n", + "\n", + " vapour_density = vapour_pressure / (CONST['R_vap'] * temperature)\n", + " a_coeff = (CONST['g'] / (CONST['R_dry'] * temperature) * (\n", + " (CONST['eps'] * CONST[\"latent_heat\"]) /\n", + " (CONST[\"specific_heat_constant_pressure\"] * temperature) - 1 )\n", + " ) # equation 7.23\n", + "\n", + " a_value = ((p_vs * water_vapor_diffusivity)\n", + " / (CONST['water_density'] * temperature * CONST['R_dry']))\n", + " b_value = (((thermal_conductivity * temperature)/\n", + " (CONST['water_density'] * CONST[\"latent_heat\"]))\n", + " * ((CONST[\"latent_heat\"] / (CONST['R_dry'] * temperature)) - 1))\n", + "\n", + " g_term = a_value + b_value\n", + "\n", + " deriv = [None] * len(state)\n", + "\n", + " condensation_rate = 0 / si.s\n", + " drop_mass = 0 * si.kg\n", + " for i in range(parameters['n_sd']):\n", + " radius = state[i]* parameters['units'][i]\n", + " multiplicity = parameters['sampled_spectrum']['multiplicity'][i]\n", + "\n", + " drop_mass += multiplicity * CONST['water_density'] * 4/3 * np.pi * radius**3\n", + "\n", + " saturation_eq = equil_saturation(\n", + " radius = radius,\n", + " kappa = parameters['kappa'],\n", + " dry_radius = parameters['sampled_spectrum']['dry radius'][i],\n", + " temperature = temperature\n", + " )\n", + " deriv[i] = (g_term/radius * (saturation - saturation_eq) /\n", + " parameters[\"updraft\"]) # equation 7.18\n", + "\n", + " condensation_rate += (condensation_rate_constant *\n", + " (multiplicity * radius**2 * deriv[i] * parameters[\"updraft\"])\n", + " )\n", + " mixing_ratio = ((parameters[\"mass of vapour at t=0\"] -\n", + " drop_mass + parameters[\"mass of droplets at t=0\"]) /\n", + " parameters['mass_of_dry_air'])\n", + "\n", + " dry_air_density = vapour_density / mixing_ratio\n", + " total_density = vapour_density + dry_air_density\n", + "\n", + " pressure = vapour_pressure + (dry_air_density * temperature * CONST['R_dry'])\n", + "\n", + " b_coeff = total_density * (\n", + " (CONST['R_dry'] * temperature / (CONST['eps'] * p_vs))\n", + " +\n", + " (CONST['eps'] * CONST[\"latent_heat\"] ** 2 /\n", + " (pressure * temperature * CONST[\"specific_heat_constant_pressure\"])\n", + " )\n", + " ) # equation 7.24\n", + "\n", + " deriv[IDX_S] = ((a_coeff * parameters['updraft'] +\n", + " b_coeff * condensation_rate) /\n", + " parameters[\"updraft\"]) # equation 7.22\n", + " deriv[IDX_T] = (- ((CONST[\"latent_heat\"] / CONST[\"specific_heat_constant_pressure\"])\n", + " * condensation_rate) - ((CONST['g'] / CONST[\"specific_heat_constant_pressure\"])\n", + " * parameters[\"updraft\"])) / parameters[\"updraft\"]\n", + "\n", + " return [\n", + " deriv[i].to(parameters['units'][i] / si.m).magnitude\n", + " for i in range(len(parameters['units']))\n", + " ]" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "73acd4bfd7c62714", + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-04T20:06:26.676432Z", + "start_time": "2024-04-04T20:06:26.669197Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "outputs": [], + "source": [ + "def compute_derived_prameters(parameters, initial_condition):\n", + " p_vs = saturated_vapor_pressure(initial_condition[IDX_T])\n", + " parameters[\"mass of vapour at t=0\"] = (\n", + " (CONST['eps'] * initial_condition[IDX_S]) /\n", + " ((parameters['initial total pressure'] / p_vs)- 1)\n", + " ) * parameters[\"mass_of_dry_air\"]\n", + " parameters[\"mass of droplets at t=0\"] = 0\n", + "\n", + " parameters['units'] = [None] * (parameters['n_sd'] + 2)\n", + " parameters['units'][IDX_S] = si.dimensionless\n", + " parameters['units'][IDX_T] = si.K\n", + " parameters['units'][:parameters['n_sd']] = [si.m] * parameters['n_sd']\n", + "\n", + "def integrate(parameters):\n", + " initial_condition = [None] * (parameters['n_sd'] + 2)\n", + " initial_condition[IDX_S] = parameters['initial relative humidity']\n", + " for i in range(parameters['n_sd']):\n", + " initial_condition[i] = scipy.optimize.root_scalar(\n", + " f=lambda radius, *params: equil_saturation(radius * si.m, *params).magnitude,\n", + " x0=(parameters['sampled_spectrum']['dry radius'][i] / si.m).magnitude,\n", + " args=(\n", + " parameters['kappa'],\n", + " parameters['sampled_spectrum']['dry radius'][i],\n", + " parameters['initial temperature']\n", + " )\n", + " ).root * si.m\n", + " initial_condition[IDX_T] = parameters['initial temperature']\n", + "\n", + " compute_derived_prameters(parameters, initial_condition)\n", + "\n", + " num_time_steps = 5\n", + " dt_out = 5 * si.s\n", + "\n", + " displacement = [(parameters[\"updraft\"] * i * dt_out).to(si.m).magnitude\n", + " for i in range(num_time_steps + 1)]\n", + "\n", + " integration_result = scipy.integrate.solve_ivp(\n", + " fun=funct,\n", + " t_span=(displacement[0], displacement[-1]),\n", + " y0=[x.to(parameters['units'][i]).magnitude for i, x in enumerate(initial_condition)],\n", + " t_eval=list(displacement),\n", + " args=(parameters,),\n", + " method='LSODA'\n", + " )\n", + "\n", + " return displacement * si.m, [\n", + " integration_result.y[i, :] * parameters['units'][i]\n", + " for i in range(integration_result.y.shape[0])\n", + " ]" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "9d5bcaaffa772ffd", + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-04T20:06:27.244924Z", + "start_time": "2024-04-04T20:06:27.240777Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "outputs": [], + "source": [ + "def plot(displacement, integration_result_with_units):\n", + " n_sd = len(integration_result_with_units) - 2\n", + " plot_variables = (\n", + " ((IDX_S,), \"saturation\", None),\n", + " (range(n_sd), \"drop radii\", si.um),\n", + " ((IDX_T,), \"temperature\", None),\n", + " )\n", + "\n", + " _, axs = pyplot.subplots(1, len(plot_variables), sharey=True)\n", + "\n", + " for i, (variables, name, unit) in enumerate(plot_variables):\n", + " for idx in variables:\n", + " axs[i].plot(\n", + " integration_result_with_units[idx],\n", + " displacement,\n", + " )\n", + " axs[i].set_title(f\"{name} vs. height\")\n", + " if unit:\n", + " axs[i].xaxis.set_units(unit)\n", + "\n", + " for subplot in axs:\n", + " subplot.grid()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "32d6bc6d68c8fe9b", + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-04T20:06:27.922964Z", + "start_time": "2024-04-04T20:06:27.894205Z" + }, + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "985ec0aebaf54e2ebc39697256f02e09", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(FloatSlider(value=1.28, description='kappa', max=3.84, min=-1.28), IntSlider(value=2, de…" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "b8e9c388433f451db7e05e3ba8d4ca25", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Output()" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "output = ipywidgets.Output()\n", + "\n", + "def parcel(kappa, updraft_m_s, n_sd):\n", + " with output:\n", + " IPython.display.clear_output(wait=True)\n", + "\n", + " parameters = {\n", + " 'updraft': updraft_m_s * si.m / si.s,\n", + " 'kappa': kappa,\n", + " 'total_aerosol_number': 1e6,\n", + " 'lognormal_geometric_stdev': 1.2,\n", + " 'lognormal_median_radius': .05 * si.um,\n", + " 'mass_of_dry_air': 1 * si.kg,\n", + " 'initial total pressure': 900 *si.hPa,\n", + " 'initial relative humidity': .75 * si.dimensionless,\n", + " 'initial temperature': 233 * si.K,\n", + " 'n_sd': n_sd\n", + " }\n", + "\n", + " parameters['sampled_spectrum'] = sample_spectrum(parameters)\n", + "\n", + " displacement, integration_result_with_units = integrate(parameters)\n", + " plot(displacement, integration_result_with_units)\n", + "\n", + " pyplot.show()\n", + "\n", + "ipywidgets.interact_manual(parcel, kappa=1.28, updraft_m_s=2, n_sd=32)\n", + "IPython.display.display(output)\n", + "# TODO: progress bar" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "141cc874f8f24b0f", + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d3ee69f4-185d-47f3-9556-61f579ba4c02", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}