PySDM_examples.Arabas_et_al_2023
box-model and 2D prescribed-flow immersion-freezing examples based on Arabas et al. 2023
aida.ipynb:
copula_hello.ipynb:
fig_2.ipynb:
fig_11.ipynb:
{
"cells": [
{
"cell_type": "markdown",
"id": "e75404edab7940f9",
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"source": [
" fig_A2.ipynb:
figs_3_and_7_and_8.ipynb:
figs_5_and_6.ipynb:
\n",
""
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "6d39ff74",
"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": "b19713b4",
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"import string\n",
"import subprocess\n",
"import platform\n",
"import pathlib\n",
"\n",
"import numpy as np\n",
"\n",
"from scipy.io import netcdf_file\n",
"from scipy.ndimage import uniform_filter1d\n",
"from matplotlib import pyplot\n",
"\n",
"from PySDM.exporters import NetCDFExporter, VTKExporter\n",
"from PySDM_examples.utils import ProgBarController\n",
"from open_atmos_jupyter_utils import show_plot\n",
"import PySDM.products as PySDM_products\n",
"from PySDM.physics import si\n",
"from PySDM import Formulae\n",
"from PySDM.initialisation import spectra\n",
"from PySDM.dynamics import Freezing\n",
"\n",
"from PySDM_examples.Arabas_et_al_2015 import Settings\n",
"from PySDM_examples.Szumowski_et_al_1998 import Simulation, Storage\n",
"from PySDM_examples.Arabas_et_al_2023.commons import FREEZING_CONSTANTS, LOGNORMAL_MODE_SURF_A, LOGNORMAL_SGM_G"
]
},
{
"cell_type": "markdown",
"id": "083e1887-fd91-4655-912e-43c8cf9f13d8",
"metadata": {},
"source": [
"## Simulations"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "b3ac5007",
"metadata": {},
"outputs": [],
"source": [
"lognormal_mode_A = LOGNORMAL_MODE_SURF_A\n",
"lognormal_sgm_g = LOGNORMAL_SGM_G\n",
"inp_frac = 1 / (270 + 45 + 1)\n",
"\n",
"conc_cld_unit = '1/cc'\n",
"conc_ice_unit = '1/l'\n",
"cool_rate_unit = 'K/min'\n",
"wall_time_unit = 'ms'\n",
"\n",
"N_REALISATIONS = 2"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "8ac0108f",
"metadata": {},
"outputs": [],
"source": [
"runs = (\n",
" {'settings': {'rhod_w_max': 2.0 * si.m/si.ssi.kg/si.m3, 'freezing_singular': False}},\n",
" {'settings': {'rhod_w_max': 2.0 * si.m/si.s*si.kg/si.m3, 'freezing_singular': True}},\n",
" {'settings': {'rhod_w_max': 1.0 * si.m/si.ssi.kg/si.m3, 'freezing_singular': False}},\n",
" {'settings': {'rhod_w_max': 1.0 * si.m/si.s*si.kg/si.m3, 'freezing_singular': True}},\n",
" {'settings': {'rhod_w_max': 0.5 * si.m/si.ssi.kg/si.m3, 'freezing_singular': False}},\n",
" {'settings': {'rhod_w_max': 0.5 * si.m/si.s*si.kg/si.m3, 'freezing_singular': True}},\n",
")\n",
"runs = tuple({'settings': {run['settings'], 'seed': seed}} for run in runs for seed in range(N_REALISATIONS))"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "6a3c69a9",
"metadata": {},
"outputs": [],
"source": [
"products = (\n",
" PySDM_products.DynamicWallTime(\n",
" 'Condensation', name='Condensation_wall_time', unit=wall_time_unit\n",
" ),\n",
" PySDM_products.DynamicWallTime(\n",
" 'Displacement', name='Displacement_wall_time', unit=wall_time_unit\n",
" ),\n",
" PySDM_products.DynamicWallTime(\n",
" 'Freezing', name='Freezing_wall_time', unit=wall_time_unit\n",
" ),\n",
" PySDM_products.DynamicWallTime(\n",
" 'EulerianAdvection', name='EulerianAdvection_wall_time', unit=wall_time_unit\n",
" ),\n",
" PySDM_products.ParticleConcentration(\n",
" radius_range=(-np.inf, 0*si.um), name='n_i', unit=conc_ice_unit, stp=True\n",
" ),\n",
" PySDM_products.ParticleConcentration(\n",
" radius_range=(1*si.um, np.inf), name='n_c', unit=conc_cld_unit, stp=True\n",
" ),\n",
" PySDM_products.CoolingRate(\n",
" unit=cool_rate_unit\n",
" ),\n",
" PySDM_products.IceNucleiConcentration(\n",
" name='n_inp', unit=conc_ice_unit, stp=True\n",
" ),\n",
" PySDM_products.FrozenParticleConcentration(\n",
" name='n_frozen_aerosols',\n",
" unit=conc_ice_unit,\n",
" count_activated=False,\n",
" count_unactivated=True,\n",
" stp=True\n",
" ),\n",
" PySDM_products.FrozenParticleConcentration(\n",
" name='n_frozen_droplets',\n",
" unit=conc_ice_unit,\n",
" count_activated=True,\n",
" count_unactivated=False,\n",
" stp=True\n",
" )\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "0b25f1e7-4887-4012-bf3b-13f6dc6d672f",
"metadata": {},
"outputs": [],
"source": [
"class SpinUp:\n",
" \"\"\" enables freezing dynamic after a given number of steps \"\"\"\n",
" def __init__(self, particulator, spin_up_steps):\n",
" self.spin_up_steps = spin_up_steps\n",
" particulator.observers.append(self)\n",
" self.particulator = particulator\n",
" self.set(Freezing, \"enable\", False)\n",
"\n",
" def notify(self):\n",
" if self.particulator.n_steps == self.spin_up_steps:\n",
" self.set(Freezing, \"enable\", True)\n",
"\n",
" def set(self, dynamic, attr, val):\n",
" name = dynamic.__name__\n",
" if name in self.particulator.dynamics:\n",
" setattr(self.particulator.dynamics[name], attr, val)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "944f6247",
"metadata": {},
"outputs": [],
"source": [
"for i, run in enumerate(runs):\n",
" folder = f\"output/rhod_w_max={run['settings']['rhod_w_max']}_singular={run['settings']['freezing_singular']}_seed={run['settings']['seed']}\"\n",
" os.makedirs(folder, exist_ok=True)\n",
" \n",
" run['ncfile'] = f'{folder}/out.nc'\n",
"\n",
" formulae = Formulae(\n",
" particle_shape_and_density=\"MixedPhaseSpheres\",\n",
" freezing_temperature_spectrum='Niemand_et_al_2012',\n",
" heterogeneous_ice_nucleation_rate='ABIFM',\n",
" constants=FREEZING_CONSTANTS[\"dust\"],\n",
" seed=run['settings']['seed'],\n",
" )\n",
" settings = Settings(formulae)\n",
" settings.dt = 2.5 * si.s\n",
" settings.output_interval = settings.dt * 12\n",
" settings.simulation_time = 6000 * si.second if 'CI' not in os.environ else 2 * settings.output_interval\n",
" settings.spin_up_time = 600 * si.second\n",
" settings.size = (1500, 500)\n",
" settings.n_sd_per_gridbox = 32\n",
" settings.grid = (60, 20)\n",
" settings.th_std0 -= 33.3 * si.kelvins\n",
" settings.initial_water_vapour_mixing_ratio -= 6.66 * si.grams / si.kilogram\n",
" \n",
" settings.processes['coalescence'] = False\n",
" settings.processes['sedimentation'] = False\n",
" settings.processes['freezing'] = True\n",
" settings.freezing_inp_spec = spectra.Lognormal(\n",
" norm_factor=1,\n",
" m_mode=lognormal_mode_A,\n",
" s_geom=lognormal_sgm_g\n",
" )\n",
" settings.freezing_inp_frac = inp_frac\n",
" settings.freezing_thaw = True\n",
"\n",
" settings.kappa = 0.61\n",
" settings.mode_1 = spectra.Lognormal(\n",
" norm_factor=(270 + 270/315) / si.centimetre3 / formulae.constants.rho_STP,\n",
" m_mode=0.03 * si.micrometre,\n",
" s_geom=1.28,\n",
" )\n",
" settings.mode_2 = spectra.Lognormal(\n",
" norm_factor=(45 + 45/315) / si.centimetre3 / formulae.constants.rho_STP,\n",
" m_mode=0.14 * si.micrometre,\n",
" s_geom=1.75,\n",
" )\n",
" settings.spectrum_per_mass_of_dry_air = spectra.Sum((settings.mode_1, settings.mode_2))\n",
" \n",
" for key, value in run['settings'].items(): \n",
" if key != 'seed':\n",
" assert hasattr(settings, key)\n",
" setattr(settings, key, value)\n",
"\n",
" storage = Storage()\n",
" simulation = Simulation(settings, storage, SpinUp=SpinUp)\n",
" simulation.reinit(products)\n",
"\n",
" vtk_exporter = VTKExporter(path=folder) \n",
" simulation.run(ProgBarController(f\"run {i+1}/{len(runs)}\"), vtk_exporter=vtk_exporter)\n",
" vtk_exporter.write_pvd()\n",
"\n",
" ncdf_exporter = NetCDFExporter(storage, settings, simulation, run['ncfile'])\n",
" ncdf_exporter.run(ProgBarController('netCDF'))"
]
},
{
"cell_type": "markdown",
"id": "27d1536b-6c7d-4e68-91e9-6a1fe2742026",
"metadata": {},
"source": [
"## Fig 11"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "241f48c4",
"metadata": {},
"outputs": [],
"source": [
"def label(settings_arg):\n",
" tmp = str({k.replace('condensation_', ''):\n",
" f\"{v:.1e}\" if isinstance(v, float) else\n",
" str(v).zfill(2) if isinstance(v, int) else\n",
" v for k, v in settings_arg.items()})\n",
" return tmp\\n",
" .replace('{', '')\\n",
" .replace('}', '')\\n",
" .replace(\"'\", '')\\n",
" .replace('rhod_w_max:', '$w_{max}\\approx$')\\n",
" .replace('e+00', r' m/s')\\n",
" .replace('5.0e-01', '0.5 m/s')\\n",
" .replace('freezing_singular: True', r'singular$\,\,\,$')\\n",
" .replace('freezing_singular: False', 'time-dep')\\n",
" .replace(', seed: 00', '')"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "bba6ffb0",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"$w_{max}\approx$ 0.5 m/s, singular$\,\,\,$, seed: 01: time=16.32ms\n",
"$w_{max}\approx$ 0.5 m/s, singular$\,\,\,$: time=16.44ms\n",
"$w_{max}\approx$ 0.5 m/s, time-dep, seed: 01: time=16.73ms\n",
"$w_{max}\approx$ 0.5 m/s, time-dep: time=16.81ms\n",
"$w_{max}\approx$ 1.0 m/s, singular$\,\,\,$, seed: 01: time=18.51ms\n",
"$w_{max}\approx$ 1.0 m/s, singular$\,\,\,$: time=18.69ms\n",
"$w_{max}\approx$ 1.0 m/s, time-dep, seed: 01: time=19.21ms\n",
"$w_{max}\approx$ 1.0 m/s, time-dep: time=19.20ms\n",
"$w_{max}\approx$ 2.0 m/s, singular$\,\,\,$, seed: 01: time=19.62ms\n",
"$w_{max}\approx$ 2.0 m/s, singular$\,\,\,$: time=19.45ms\n",
"$w_{max}\approx$ 2.0 m/s, time-dep, seed: 01: time=19.52ms\n",
"$w_{max}\approx$ 2.0 m/s, time-dep: time=27.82ms\n",
"$w_{max}\approx$ 0.5 m/s, singular$\,\,\,$, seed: 01: time=16.32ms\n",
"$w_{max}\approx$ 0.5 m/s, singular$\,\,\,$: time=16.44ms\n",
"$w_{max}\approx$ 0.5 m/s, time-dep, seed: 01: time=16.73ms\n",
"$w_{max}\approx$ 0.5 m/s, time-dep: time=16.81ms\n",
"$w_{max}\approx$ 1.0 m/s, singular$\,\,\,$, seed: 01: time=18.51ms\n",
"$w_{max}\approx$ 1.0 m/s, singular$\,\,\,$: time=18.69ms\n",
"$w_{max}\approx$ 1.0 m/s, time-dep, seed: 01: time=19.21ms\n",
"$w_{max}\approx$ 1.0 m/s, time-dep: time=19.20ms\n",
"$w_{max}\approx$ 2.0 m/s, singular$\,\,\,$, seed: 01: time=19.62ms\n",
"$w_{max}\approx$ 2.0 m/s, singular$\,\,\,$: time=19.45ms\n",
"$w_{max}\approx$ 2.0 m/s, time-dep, seed: 01: time=19.52ms\n",
"$w_{max}\approx$ 2.0 m/s, time-dep: time=27.82ms\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
"
\"), HTML(value=\"
1""" 2box-model and 2D prescribed-flow immersion-freezing examples based on 3[Arabas et al. 2023](https://arxiv.org/abs/2308.05015) 4 5aida.ipynb: 6.. include:: ./aida.ipynb.badges.md 7 8copula_hello.ipynb: 9.. include:: ./copula_hello.ipynb.badges.md 10 11fig_2.ipynb: 12.. include:: ./fig_2.ipynb.badges.md 13 14fig_11.ipynb: 15.. include:: ./figs_10_and_11_and_animations.ipynb 16 17fig_A2.ipynb: 18.. include:: ./fig_A2.ipynb.badges.md 19 20figs_3_and_7_and_8.ipynb: 21.. include:: ./figs_3_and_7_and_8.ipynb.badges.md 22 23figs_5_and_6.ipynb: 24.. include:: ./figs_5_and_6.ipynb.badges.md 25""" 26 27from .make_particulator import make_particulator 28from .plots import make_freezing_spec_plot, make_pdf_plot, make_temperature_plot 29from .run_simulation import run_simulation