{
"cells": [
{
"cell_type": "markdown",
"id": "20b1a2f1",
"metadata": {},
"source": [
"[](https://github.com/open-atmos/PyPartMC/blob/main/examples/mie_optical.ipynb) \n",
"[](https://colab.research.google.com/github/open-atmos/PyPartMC/blob/main/examples/mie_optical.ipynb) \n",
"[](https://mybinder.org/v2/gh/open-atmos/PyPartMC.git/main?urlpath=lab/tree/examples/mie_optical.ipynb) \n",
"[](https://jupyterhub.arm.gov/hub/user-redirect/git-pull?repo=https%3A//github.com/open-atmos/PyPartMC&branch=main&urlPath=) (requires [logging in with ARM account](https://www.arm.gov/capabilities/computing-resources) and directing Jupyter to a notebook within the cloned repo)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "159edeb4",
"metadata": {},
"outputs": [],
"source": [
"# This file is a part of PyPartMC licensed under the GNU General Public License v3\n",
"# Copyright (C) 2023 University of Illinois Urbana-Champaign\n",
"# Authors:\n",
"# - https://github.com/compdyn/partmc/graphs/contributors\n",
"# - https://github.com/open-atmos/PyPartMC/graphs/contributors"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "4f8359c2",
"metadata": {},
"outputs": [],
"source": [
"import sys\n",
"import os\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('PyPartMC', 'git+https://github.com/bsumlin/PyMieScatt.git')\n",
"elif 'JUPYTER_IMAGE' in os.environ and '.arm.gov' in os.environ['JUPYTER_IMAGE']:\n",
" !pip --quiet install PyPartMC git+https://github.com/bsumlin/PyMieScatt.git open_atmos_jupyter_utils\n",
" _pypartmc_path = !pip show PyPartMC | fgrep Location | cut -f2 -d' '\n",
" sys.path.extend(_pypartmc_path if _pypartmc_path[0] not in sys.path else [])"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "b494ea6e",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"# https://pymiescatt.readthedocs.io/en/latest/ and https://github.com/bsumlin/PyMieScatt\n",
"import PyMieScatt as ps\n",
"from open_atmos_jupyter_utils import show_plot\n",
"import PyPartMC as ppmc\n",
"from PyPartMC import si"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "1747cc4c",
"metadata": {},
"outputs": [],
"source": [
"plt.rcParams.update({'font.size': 9})\n",
"plt.rcParams.update({'figure.figsize': (3.08,2.5)})\n",
"plt.rcParams.update({\"axes.grid\" : True})"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "5aae092f",
"metadata": {},
"outputs": [],
"source": [
"def aero_state_compute_optical(aero_state_optical):\n",
" wl = 550.0 # unit: nm\n",
" refr_shell = 1.52+0*1j\n",
" refr_core = 1.82+0.74*1j\n",
" diams_core = np.array(aero_state_optical.diameters(include=[\"BC\"])) * 1e9 # unit: nm\n",
" diams_total = np.array(aero_state_optical.diameters()) * 1e9 # unit: nm\n",
" \n",
" qsca_part = np.zeros(len(aero_state_optical))\n",
" qabs_part = np.zeros(len(aero_state_optical))\n",
" for i_part in range(len(aero_state_optical)):\n",
" val = ps.MieQCoreShell(refr_core,\n",
" refr_shell,\n",
" wl,\n",
" diams_core[i_part],\n",
" diams_total[i_part],\n",
" asDict=True\n",
" )\n",
" qsca_part[i_part] = val['Qsca']\n",
" qabs_part[i_part] = val['Qabs']\n",
" \n",
" cross_section = np.pi * (diams_total / 2 / 1e9)**2\n",
" \n",
" num_concs = aero_state_optical.num_concs\n",
" B_abs = np.sum(qabs_part * cross_section * num_concs)\n",
" B_sca = np.sum(qsca_part * cross_section * num_concs)\n",
"\n",
" return (diams_total, qsca_part, qabs_part, B_sca, B_abs)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "df8d1a84",
"metadata": {},
"outputs": [],
"source": [
"gas_data = ppmc.GasData((\"H2SO4\",\"HNO3\",\"HCl\",\"NH3\",\"NO\",\"NO2\", \"NO3\",\n",
" \"N2O5\", \"HONO\", \"HNO4\", \"O3\", \"O1D\", \"O3P\", \"OH\",\n",
" \"HO2\", \"H2O2\", \"CO\", \"SO2\", \"CH4\", \"C2H6\", \"CH3O2\", \n",
" \"ETHP\", \"HCHO\", \"CH3OH\", \"ANOL\", \"CH3OOH\", \"ETHOOH\",\n",
" \"ALD2\", \"HCOOH\", \"RCOOH\", \"C2O3\", \"PAN\", \"ARO1\", \"ARO2\",\n",
" \"ALK1\", \"OLE1\", \"API1\", \"API2\", \"LIM1\", \"LIM2\", \"PAR\", \"AONE\",\n",
" \"MGLY\", \"ETH\", \"OLET\", \"OLEI\", \"TOL\", \"XYL\", \"CRES\", \"TO2\",\n",
" \"CRO\", \"OPEN\", \"ONIT\", \"ROOH\", \"RO2\", \"ANO2\", \"NAP\", \"XO2\",\n",
" \"XPAR\", \"ISOP\", \"ISOPRD\", \"ISOPP\", \"ISOPN\", \"ISOPO2\", \"API\",\n",
" \"LIM\", \"DMS\", \"MSA\", \"DMSO\", \"DMSO2\", \"CH3SO2H\", \"CH3SCH2OO\", \n",
" \"CH3SO2\", \"CH3SO3\", \"CH3SO2OO\", \"CH3SO2CH2OO\", \"SULFHOX\"\n",
" ))"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "a2d7bad8",
"metadata": {},
"outputs": [],
"source": [
"env_state = ppmc.EnvState(\n",
" {\n",
" \"rel_humidity\": 0.95,\n",
" \"latitude\": 0,\n",
" \"longitude\": 0,\n",
" \"altitude\": 0 * si.m,\n",
" \"start_time\": 21600 * si.s,\n",
" \"start_day\": 200,\n",
" }\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "815283a0",
"metadata": {},
"outputs": [],
"source": [
"aero_data = ppmc.AeroData(\n",
" (\n",
" # density ions in soln (1) molecular weight kappa (1)\n",
" # | | | |\n",
" {\"SO4\": [1800 * si.kg / si.m**3, 1, 96.0 * si.g / si.mol, 0.00]},\n",
" {\"NO3\": [1800 * si.kg / si.m**3, 1, 62.0 * si.g / si.mol, 0.00]},\n",
" {\"Cl\": [2200 * si.kg / si.m**3, 1, 35.5 * si.g / si.mol, 0.00]},\n",
" {\"NH4\": [1800 * si.kg / si.m**3, 1, 18.0 * si.g / si.mol, 0.00]},\n",
" {\"MSA\": [1800 * si.kg / si.m**3, 0, 95.0 * si.g / si.mol, 0.53]},\n",
" {\"ARO1\": [1400 * si.kg / si.m**3, 0, 150.0 * si.g / si.mol, 0.10]},\n",
" {\"ARO2\": [1400 * si.kg / si.m**3, 0, 150.0 * si.g / si.mol, 0.10]},\n",
" {\"ALK1\": [1400 * si.kg / si.m**3, 0, 140.0 * si.g / si.mol, 0.10]},\n",
" {\"OLE1\": [1400 * si.kg / si.m**3, 0, 140.0 * si.g / si.mol, 0.10]},\n",
" {\"API1\": [1400 * si.kg / si.m**3, 0, 184.0 * si.g / si.mol, 0.10]},\n",
" {\"API2\": [1400 * si.kg / si.m**3, 0, 184.0 * si.g / si.mol, 0.10]},\n",
" {\"LIM1\": [1400 * si.kg / si.m**3, 0, 200.0 * si.g / si.mol, 0.10]},\n",
" {\"LIM2\": [1400 * si.kg / si.m**3, 0, 200.0 * si.g / si.mol, 0.10]},\n",
" {\"CO3\": [2600 * si.kg / si.m**3, 1, 60.0 * si.g / si.mol, 0.00]},\n",
" {\"Na\": [2200 * si.kg / si.m**3, 1, 23.0 * si.g / si.mol, 0.00]},\n",
" {\"Ca\": [2600 * si.kg / si.m**3, 1, 40.0 * si.g / si.mol, 0.00]},\n",
" {\"OIN\": [2600 * si.kg / si.m**3, 0, 1.0 * si.g / si.mol, 0.10]},\n",
" {\"OC\": [1400 * si.kg / si.m**3, 0, 1.0 * si.g / si.mol, 0.10]},\n",
" {\"BC\": [1800 * si.kg / si.m**3, 0, 1.0 * si.g / si.mol, 0.00]},\n",
" {\"H2O\": [1000 * si.kg / si.m**3, 0, 18.0 * si.g / si.mol, 0.00]},\n",
" )\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "c822823d",
"metadata": {},
"outputs": [],
"source": [
"gas_state = ppmc.GasState(gas_data)\n",
"\n",
"input_gas_state = (\n",
" {\"NO\": [0.1E+00]},\n",
")\n",
"\n",
"gas_state.mix_rats = input_gas_state"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "58706963",
"metadata": {},
"outputs": [],
"source": [
"times = [0 * si.s]\n",
"back_gas = [{\"time\": times},\n",
" {\"rate\": [1.5e-5 / si.s]},\n",
" {\"NO\": [0.1E+00]},\n",
" ]\n",
"\n",
"gas_emit_times = [0, 3600, 7200, 10800, 14400, 18000, 21600, 25200, 28800, 32400, 36000,\n",
" 39600, 43200, 46800, 50400, 54000, 57600, 61200, 64800, 68400, 72000,\n",
" 75600, 79200, 82800, 90000, 93600, 97200, 100800, 104400, 108000]\n",
"\n",
"gas_emit_rates = np.zeros(len(gas_emit_times))\n",
"gas_emit_rates[0:12] = .5\n",
"\n",
"SO2 = [4.234E-09, 5.481E-09, 5.089E-09, 5.199E-09, 5.221E-09, 5.284E-09, 5.244E-09,\n",
" 5.280E-09, 5.560E-09, 5.343E-09, 4.480E-09, 3.858E-09, 3.823E-09, 3.607E-09,\n",
" 3.533E-09, 3.438E-09, 2.866E-09, 2.667E-09, 2.636E-09, 2.573E-09, 2.558E-09,\n",
" 2.573E-09, 2.715E-09, 3.170E-09, 4.2344E-09, 5.481E-09, 5.089E-09, 5.199E-09,\n",
" 5.221E-09, 5.284E-09]\n",
"\n",
"emit_gas = [\n",
" {\"time\": gas_emit_times},\n",
" {\"rate\": list(gas_emit_rates)},\n",
" {\"SO2\": SO2},\n",
"]\n",
"\n",
"AERO_DIST_BACKGROUND = {\n",
" \"back_small\": {\n",
" \"mass_frac\": [{\"SO4\": [1]}, {\"OC\": [1.375]}, {\"NH4\": [0.375]}],\n",
" \"diam_type\": \"geometric\",\n",
" \"mode_type\": \"log_normal\",\n",
" \"num_conc\": 3.2e9 / si.m**3,\n",
" \"geom_mean_diam\": 0.02 * si.um,\n",
" \"log10_geom_std_dev\": 0.161,\n",
" },\n",
" \"back_large\": {\n",
" \"mass_frac\": [{\"SO4\": [1]}, {\"OC\": [1.375]}, {\"NH4\": [0.375]}],\n",
" \"diam_type\": \"geometric\",\n",
" \"mode_type\": \"log_normal\",\n",
" \"num_conc\": 2.9e9 / si.m**3,\n",
" \"geom_mean_diam\": 0.16 * si.um,\n",
" \"log10_geom_std_dev\": 0.217,\n",
" },\n",
"}\n",
"\n",
"AERO_DIST_EMIT = {\n",
" \"gasoline\": {\n",
" \"mass_frac\": [{\"OC\": [0.8]}, {\"BC\": [0.2]}],\n",
" \"diam_type\": \"geometric\",\n",
" \"mode_type\": \"log_normal\",\n",
" \"num_conc\": 5e7 / si.m**3,\n",
" \"geom_mean_diam\": 5e-8 * si.m,\n",
" \"log10_geom_std_dev\": 0.24,\n",
" },\n",
" \"diesel\": {\n",
" \"mass_frac\": [{\"OC\": [0.3]}, {\"BC\": [0.7]}],\n",
" \"diam_type\": \"geometric\",\n",
" \"mode_type\": \"log_normal\",\n",
" \"num_conc\": 1.6e8 / si.m**3,\n",
" \"geom_mean_diam\": 5e-8 * si.m,\n",
" \"log10_geom_std_dev\": 0.24,\n",
" },\n",
" \"cooking\": {\n",
" \"mass_frac\": [{\"OC\": [1]}],\n",
" \"diam_type\": \"geometric\",\n",
" \"mode_type\": \"log_normal\",\n",
" \"num_conc\": 9e6 / si.m**3,\n",
" \"geom_mean_diam\": 8.64e-8 * si.m,\n",
" \"log10_geom_std_dev\": 0.28,\n",
" },\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "c6a96b7d",
"metadata": {},
"outputs": [],
"source": [
"time_timeseries = list(np.linspace(0,24*3600,25))\n",
"pressure_timeseries = list(np.ones(25) * 1e5)\n",
"temp_timeseries = [290.016,292.5, 294.5, 296.112, 297.649, 299.049, 299.684, 299.509,299.002,\n",
" 298.432, 296.943, 295.153, 293.475, 292.466, 291.972, 291.96, 291.512,\n",
" 291.481, 290.5, 290.313, 290.317, 290.362, 290.245, 290.228, 291.466]\n",
"height_timeseries = [171.045, 228.210, 296.987, 366.002, 410.868, 414.272, 417.807,414.133,\n",
" 397.465, 376.864, 364.257, 352.119, 338.660, 322.028, 305.246, 258.497, \n",
" 240.478, 187.229, 145.851, 128.072, 110.679, 97.628, 93.034, 93.034, 93.034]"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "920e41e3",
"metadata": {},
"outputs": [],
"source": [
"scenario = ppmc.Scenario(\n",
" gas_data,\n",
" aero_data,\n",
" {\n",
" \"temp_profile\": [{\"time\": time_timeseries}, {\"temp\": temp_timeseries}],\n",
" \"pressure_profile\": [\n",
" {\"time\": time_timeseries},\n",
" {\"pressure\": pressure_timeseries},\n",
" ],\n",
" \"height_profile\": [{\"time\": time_timeseries}, {\"height\": height_timeseries}],\n",
" \"gas_emissions\": emit_gas,\n",
" \"gas_background\": back_gas,\n",
" \"aero_emissions\": [\n",
" {\"time\": [0 * si.s, 12 * 3600 * si.s]},\n",
" {\"rate\": [1 / si.s, 0 / si.s]},\n",
" {\"dist\": [[AERO_DIST_EMIT],[AERO_DIST_EMIT]]},\n",
" ],\n",
" \"aero_background\": [\n",
" {\"time\": [0 * si.s]},\n",
" {\"rate\": [1.5e-5 / si.s]},\n",
" {\"dist\": [[AERO_DIST_BACKGROUND]]},\n",
" ],\n",
" \"loss_function\": \"none\",\n",
" },\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "9781ca2f",
"metadata": {},
"outputs": [],
"source": [
"T_INITIAL = 0.0\n",
"scenario.init_env_state(env_state, T_INITIAL)\n",
"\n",
"AERO_DIST_INIT = [\n",
" {\n",
" \"init_small\": {\n",
" \"mass_frac\": [{\"SO4\": [1]}, {\"OC\": [1.375]}, {\"NH4\": [0.375]}],\n",
" \"diam_type\": \"geometric\",\n",
" \"mode_type\": \"log_normal\",\n",
" \"num_conc\": 3.2e9 / si.m**3,\n",
" \"geom_mean_diam\": 0.02 * si.um,\n",
" \"log10_geom_std_dev\": 0.161,\n",
" },\n",
" \"init_large\": {\n",
" \"mass_frac\": [{\"SO4\": [1]}, {\"OC\": [1.375]}, {\"NH4\": [0.375]}],\n",
" \"diam_type\": \"geometric\",\n",
" \"mode_type\": \"log_normal\",\n",
" \"num_conc\": 2.9e9 / si.m**3,\n",
" \"geom_mean_diam\": 0.16 * si.um,\n",
" \"log10_geom_std_dev\": 0.217,\n",
" },\n",
" }\n",
"]\n",
"\n",
"aero_dist_init = ppmc.AeroDist(aero_data, AERO_DIST_INIT)\n",
"\n",
"run_part_opt = ppmc.RunPartOpt(\n",
" {\n",
" \"output_prefix\": \"urban_plume\",\n",
" \"do_coagulation\": True,\n",
" \"coag_kernel\": \"brown\",\n",
" \"t_max\": 24 * 3600 * si.s,\n",
" \"del_t\": 60 * si.s,\n",
" }\n",
")\n",
"\n",
"N_PART = 1000\n",
"aero_state = ppmc.AeroState(aero_data, N_PART, 'nummass_source')\n",
"n_added = aero_state.dist_sample(\n",
" aero_dist_init,\n",
" sample_prop=1.0,\n",
" create_time=0.0,\n",
" allow_doubling=True,\n",
" allow_halving=True,\n",
")\n",
"camp_core = ppmc.CampCore()\n",
"photolysis = ppmc.Photolysis()"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "abf48304",
"metadata": {},
"outputs": [],
"source": [
"N_BLOCKS = int(run_part_opt.t_max / 3600)\n",
"N_STEPS = int(run_part_opt.t_max / run_part_opt.del_t)\n",
"N_STEPS_PER_BLOCK = int(N_STEPS / N_BLOCKS)\n",
"Bsca = np.zeros(N_BLOCKS+1)\n",
"Babs = np.zeros(N_BLOCKS+1)\n",
"time = np.zeros(N_BLOCKS+1)\n",
"i_time = 1\n",
"i_output = 1\n",
"last_output_time = 0.\n",
"last_progress_time = 0.\n",
"diameters, qsca, qabs, Bsca[0], Babs[0] = aero_state_compute_optical(aero_state)\n",
"\n",
"for i_block in range(1,N_BLOCKS+1):\n",
" i_next = int(N_STEPS_PER_BLOCK * i_block)\n",
" (last_output_time, last_progress_time, i_output) = ppmc.run_part_timeblock(\n",
" scenario,\n",
" env_state,\n",
" aero_data,\n",
" aero_state,\n",
" gas_data,\n",
" gas_state,\n",
" run_part_opt,\n",
" camp_core,\n",
" photolysis,\n",
" i_time,\n",
" i_next,\n",
" 0,\n",
" last_output_time,\n",
" last_progress_time,\n",
" i_output\n",
" )\n",
" time[i_block] = env_state.elapsed_time\n",
" diameters, qsca, qabs, Bsca[i_block], Babs[i_block] = aero_state_compute_optical(aero_state)\n",
" i_time = i_next + 1"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "ab9b78ef",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "7d24301b1afe4937ae36b7bc32a81049",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"HTML(value=\"./tmpwa8h81i9.pdf
\")"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.scatter(diameters * 1e-3, qabs, s=2)\n",
"plt.ylabel(r\"Absorption efficiency $Q_{abs}$\")\n",
"plt.xlabel(r\"Particle diameter $D_P\\;[\\mu m]$\")\n",
"plt.xscale('log')\n",
"show_plot()"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "7f3223ed",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "d6865f9f84874debb64e2c33cf190809",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"HTML(value=\"./tmp0ssq943b.pdf
\")"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(time, Bsca * 1e6, label=r'scattering')\n",
"plt.plot(time, Babs * 1e6, label=r'absorption')\n",
"plt.ylabel(r\"Optical coefficients (Mm$^{-1}$)\")\n",
"plt.xlabel(r\"time (s)\")\n",
"plt.xlim([0,24*3600])\n",
"plt.ylim([0,160])\n",
"plt.legend()\n",
"show_plot()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "89ed5091",
"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.11.3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}