diff --git a/examples/particle_simulation_with_camp.ipynb b/examples/particle_simulation_with_camp.ipynb
new file mode 100644
index 00000000..19dfed13
--- /dev/null
+++ b/examples/particle_simulation_with_camp.ipynb
@@ -0,0 +1,5718 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "7968866b-24b8-4ed6-bc13-6035310cb4f5",
+ "metadata": {},
+ "source": [
+ "[](https://github.com/open-atmos/PyPartMC/blob/main/examples/particle_simulation_with_camp.ipynb) \n",
+ "[](https://colab.research.google.com/github/open-atmos/PyPartMC/blob/main/examples/particle_simulation_with_camp.ipynb) \n",
+ "[](https://mybinder.org/v2/gh/open-atmos/PyPartMC.git/main?urlpath=lab/tree/examples/particle_simulation_with_camp.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) 2024 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')\n",
+ "elif 'JUPYTER_IMAGE' in os.environ and '.arm.gov' in os.environ['JUPYTER_IMAGE']:\n",
+ " !pip --quiet install PyPartMC 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 json\n",
+ "import urllib\n",
+ "from collections import defaultdict\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\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": "0ec2454c-337d-4dc4-8ceb-5692948a9ef0",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N_PART = 100"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "6ff265d8-b7f8-4f54-bfa0-d14ce0f4535c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "PATHS = [\n",
+ " 'cb05_abs_tol.json',\n",
+ " 'cb05_species.json',\n",
+ " 'custom_species.json',\n",
+ " 'aerosol_phases.json',\n",
+ " 'partitioning_species_params.json',\n",
+ " 'cb05_mechanism.noR142R143.json',\n",
+ " 'tsigaridis_2_product_SOA_scheme-mechanism.json',\n",
+ " 'tsigaridis_2_product_SOA_scheme-species.json',\n",
+ "]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "cc3bdd53-5da9-45f4-9a45-b6f1cc457fb5",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "cb05_abs_tol.json\n",
+ "cb05_species.json\n",
+ "custom_species.json\n",
+ "aerosol_phases.json\n",
+ "partitioning_species_params.json\n",
+ "cb05_mechanism.noR142R143.json\n",
+ "tsigaridis_2_product_SOA_scheme-mechanism.json\n",
+ "tsigaridis_2_product_SOA_scheme-species.json\n"
+ ]
+ }
+ ],
+ "source": [
+ "CAMP_URL = (\n",
+ " 'https://raw.githubusercontent.com/open-atmos/camp/refs/heads/main'\n",
+ " '/data/CAMP_v1_paper/modal/monarch_box_modal/'\n",
+ ")\n",
+ "for path in PATHS:\n",
+ " print(path)\n",
+ " with open(path, 'w', encoding='utf-8') as fout:\n",
+ " with urllib.request.urlopen(CAMP_URL + path.replace('-', '/')) as fin:\n",
+ " json.dump(\n",
+ " json.loads(fin.read().decode('utf-8').replace('TO' + 'DO', 'TO-DO')),\n",
+ " fout,\n",
+ " indent=4\n",
+ " )"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "16377896",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "with open('config.json', 'w', encoding='utf-8') as f:\n",
+ " json.dump({\n",
+ " \"camp-files\" : [\n",
+ " \"aerosol_representation.json\",\n",
+ " *PATHS\n",
+ " ]\n",
+ " }, f, indent=4)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "5e4e5c4b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "with open('aerosol_representation.json', 'w', encoding='utf-8') as f:\n",
+ " json.dump(\n",
+ " {\n",
+ " \"camp-data\" : [\n",
+ " {\n",
+ " \"name\" : \"PartMC single particle\",\n",
+ " \"type\" : \"AERO_REP_SINGLE_PARTICLE\",\n",
+ " \"maximum computational particles\" : N_PART * 2\n",
+ " }\n",
+ " ]\n",
+ " }, f, indent=4\n",
+ " )"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "c205c307-9d78-4cc9-840e-df1b87f5798f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "T_INITIAL = 0.0\n",
+ "RUN_PART_ARGS = {}\n",
+ "RUN_PART_ARGS['env_state'] = ppmc.EnvState({\n",
+ " \"rel_humidity\": 0.,\n",
+ " \"latitude\": 0,\n",
+ " \"longitude\": 0,\n",
+ " \"altitude\": 0 * si.m,\n",
+ " \"start_time\": 21600 * si.s,\n",
+ " \"start_day\": 200,\n",
+ "})"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "id": "de15cde5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "RUN_PART_ARGS['camp_core'] = ppmc.CampCore(\"config.json\")\n",
+ "RUN_PART_ARGS['photolysis'] = ppmc.Photolysis(RUN_PART_ARGS['camp_core'])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "id": "6de20b3f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "RUN_PART_ARGS['gas_data'] = ppmc.GasData(RUN_PART_ARGS['camp_core'])\n",
+ "RUN_PART_ARGS['aero_data'] = ppmc.AeroData({}, RUN_PART_ARGS['camp_core'])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "id": "c822823d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "RUN_PART_ARGS['gas_state'] = ppmc.GasState(RUN_PART_ARGS['gas_data'])\n",
+ "RUN_PART_ARGS['gas_state'].mix_rats = (\n",
+ " {\"NO\": [0.1]},\n",
+ " {\"NO2\": [1.]},\n",
+ " {\"O3\": [5.0e1]},\n",
+ " {\"H2O2\": [1.1]},\n",
+ " {\"CO\": [2.1e2]},\n",
+ " {\"SO2\": [0.8]},\n",
+ " {\"NH3\": [0.5]},\n",
+ " {\"HCL\": [0.7]},\n",
+ " {\"CH4\": [2.2e3]},\n",
+ " {\"ETHA\": [1.]},\n",
+ " {\"FORM\": [1.2]},\n",
+ " {\"MEOH\": [1.2e-1]},\n",
+ " {\"MEPX\": [0.5]},\n",
+ " {\"ALD2\": [1.]},\n",
+ " {\"PAR\": [2.]},\n",
+ " {\"ETH\": [0.2]},\n",
+ " {\"OLE\": [2.3e-2]},\n",
+ " {\"IOLE\": [3.1e-4]},\n",
+ " {\"TOL\": [0.1]},\n",
+ " {\"XYL\": [0.1]},\n",
+ " {\"NTR\": [0.1]},\n",
+ " {\"PAN\": [0.8]},\n",
+ " {\"AACD\": [0.2]},\n",
+ " {\"ROOH\": [2.5e-2]},\n",
+ " {\"ISOP\": [5.]},\n",
+ " {\"O2\": [2.095e8]},\n",
+ " {\"N2\": [7.8e8]},\n",
+ " {\"H2\": [5.6e2]},\n",
+ " {\"M\": [1e9]}\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "id": "920e41e3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "TIME_TIMESERIES = [0, 86400]\n",
+ "RUN_PART_ARGS['scenario'] = ppmc.Scenario(\n",
+ " RUN_PART_ARGS['gas_data'],\n",
+ " RUN_PART_ARGS['aero_data'],\n",
+ " {\n",
+ " \"temp_profile\": [{\"time\": TIME_TIMESERIES}, {\"temp\": [290.016, 290.016]}],\n",
+ " \"pressure_profile\": [\n",
+ " {\"time\": TIME_TIMESERIES},\n",
+ " {\"pressure\": [1e5, 1e5]},\n",
+ " ],\n",
+ " \"height_profile\": [{\"time\": TIME_TIMESERIES}, {\"height\": [1., 1.]}],\n",
+ " \"gas_emissions\": [\n",
+ " {\"time\": [0, 43200]},\n",
+ " {\"rate\": [1., 0.]},\n",
+ " {\"SO2\": [1.06E-9, 1.06E-9]},\n",
+ " {\"NO2\": [7.56E-12, 7.56E-12]},\n",
+ " {\"NO\": [1.44E-10, 1.44E-10]},\n",
+ " {\"CO\": [1.96E-9, 1.96E-9]},\n",
+ " {\"ALD2\": [4.25E-12, 4.25E-12]},\n",
+ " {\"FORM\": [1.02E-11, 1.02E-11]},\n",
+ " {\"ETH\": [4.62E-11, 4.62E-11]},\n",
+ " {\"IOLE\": [1.49E-11, 1.49E-11]},\n",
+ " {\"OLE\": [1.49E-11, 1.49E-11]},\n",
+ " {\"TOL\": [1.53E-11, 1.53E-11]},\n",
+ " {\"XYL\": [1.40E-11, 1.40E-11]},\n",
+ " {\"PAR\": [4.27E-10, 4.27E-10]},\n",
+ " {\"ISOP\": [6.03E-12, 6.03E-12]},\n",
+ " {\"MEOH\": [5.92E-13, 5.92E-13]},\n",
+ " ],\n",
+ " \"gas_background\": [\n",
+ " {\"time\": [0 * si.s]},\n",
+ " {\"rate\": [0 / si.s]},\n",
+ " {\"NO\": [0.1]},\n",
+ " {\"NO2\": [1.]},\n",
+ " {\"O3\": [5.0E+01]},\n",
+ " {\"H2O2\": [1.1]},\n",
+ " {\"CO\": [2.1E+02]},\n",
+ " {\"SO2\": [0.8]},\n",
+ " {\"NH3\": [0.5]},\n",
+ " {\"HCL\": [0.7]},\n",
+ " {\"CH4\": [2.2E+03]},\n",
+ " {\"ETHA\": [1.]},\n",
+ " {\"FORM\": [1.2]},\n",
+ " {\"MEOH\": [1.2E-1]},\n",
+ " {\"MEPX\": [0.5]},\n",
+ " {\"ALD2\": [1.]},\n",
+ " {\"PAR\": [2.]},\n",
+ " {\"ETH\": [0.2]},\n",
+ " {\"OLE\": [2.3E-2]},\n",
+ " {\"IOLE\": [3.1E-4]},\n",
+ " {\"TOL\": [0.1]},\n",
+ " {\"XYL\": [0.1]},\n",
+ " {\"NTR\": [0.1]},\n",
+ " {\"PAN\": [0.8]},\n",
+ " {\"AACD\": [0.2]},\n",
+ " {\"ROOH\": [2.5E-2]},\n",
+ " {\"ISOP\": [5.]},\n",
+ " {\"O2\": [2.095E+08]},\n",
+ " {\"N2\": [7.8E+08]},\n",
+ " {\"H2\": [5.6E+02]},\n",
+ " {\"M\": [1.0E+09]}\n",
+ " ],\n",
+ " \"aero_emissions\": [\n",
+ " {\"time\": [0 * si.s]},\n",
+ " {\"rate\": [0 / si.s]},\n",
+ " {\"dist\": [[{\n",
+ " \"gasoline\": {\n",
+ " \"mass_frac\": [{\"organic_matter.POA\": [1]}],\n",
+ " \"diam_type\": \"geometric\",\n",
+ " \"mode_type\": \"log_normal\",\n",
+ " \"num_conc\": 0.0 / si.m**3,\n",
+ " \"geom_mean_diam\": 5e-8 * si.m,\n",
+ " \"log10_geom_std_dev\": 0.24,\n",
+ " },\n",
+ " }]]},\n",
+ " ],\n",
+ " \"aero_background\": [\n",
+ " {\"time\": [0 * si.s]},\n",
+ " {\"rate\": [0 / si.s]},\n",
+ " {\"dist\": [[{\n",
+ " \"back_small\": {\n",
+ " \"mass_frac\": [{\"organic_matter.POA\": [1]}],\n",
+ " \"diam_type\": \"geometric\",\n",
+ " \"mode_type\": \"log_normal\",\n",
+ " \"num_conc\": 0 / si.m**3,\n",
+ " \"geom_mean_diam\": 0.02 * si.um,\n",
+ " \"log10_geom_std_dev\": 0.161,\n",
+ " },\n",
+ " }]]},\n",
+ " ],\n",
+ " \"loss_function\": \"none\",\n",
+ " },\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "id": "6722ba83",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "RUN_PART_ARGS['scenario'].init_env_state(RUN_PART_ARGS['env_state'], T_INITIAL)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "id": "96ba54c5-dc33-4e22-b774-390cec4be69f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "AERO_DIST_INIT = ppmc.AeroDist(RUN_PART_ARGS['aero_data'], [\n",
+ " {\n",
+ " \"init_small\": {\n",
+ " \"mass_frac\": [{\"organic_matter.POA\": [1]}],\n",
+ " \"diam_type\": \"geometric\",\n",
+ " \"mode_type\": \"log_normal\",\n",
+ " \"num_conc\": 3.2e9 / si.m**3,\n",
+ " \"geom_mean_diam\": 2.0e-8 * si.m,\n",
+ " \"log10_geom_std_dev\": 0.161,\n",
+ " },\n",
+ " \"init_large\": {\n",
+ " \"mass_frac\": [{\"organic_matter.POA\": [1]}],\n",
+ " \"diam_type\": \"geometric\",\n",
+ " \"mode_type\": \"log_normal\",\n",
+ " \"num_conc\": 2.9e9 / si.m**3,\n",
+ " \"geom_mean_diam\": 1.16e-7 * si.m,\n",
+ " \"log10_geom_std_dev\": 0.217,\n",
+ " },\n",
+ " \"init_coarse\": {\n",
+ " \"mass_frac\": [{\"organic_matter.POA\": [1]}],\n",
+ " \"diam_type\": \"geometric\",\n",
+ " \"mode_type\": \"log_normal\",\n",
+ " \"num_conc\": 0.3e6 / si.m**3,\n",
+ " \"geom_mean_diam\": 1.8e-6 * si.m,\n",
+ " \"log10_geom_std_dev\": 0.38021124171160603,\n",
+ " }\n",
+ " }\n",
+ " ]\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "id": "d8d9c8fd",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "52"
+ ]
+ },
+ "execution_count": 16,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "RUN_PART_ARGS['run_part_opt'] = ppmc.RunPartOpt({\n",
+ " \"output_prefix\": \"urban_plume\",\n",
+ " \"do_coagulation\": False,\n",
+ " \"t_max\": 86400 * si.s,\n",
+ " \"del_t\": 6 * si.s,\n",
+ " \"do_camp_chem\": True,\n",
+ "})\n",
+ "\n",
+ "RUN_PART_ARGS['aero_state'] = ppmc.AeroState(\n",
+ " RUN_PART_ARGS['aero_data'],\n",
+ " N_PART,\n",
+ " 'nummass_source',\n",
+ " RUN_PART_ARGS['camp_core']\n",
+ ")\n",
+ "\n",
+ "RUN_PART_ARGS['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",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "id": "5f3665a5-be3c-49e4-b1d0-47f3e2b82c30",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "DIAM_GRID = ppmc.BinGrid(30, \"log\", 1e-9, 1e-5)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "id": "0f7c29fe",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "OUTPUT = defaultdict(list)\n",
+ "last = defaultdict(float)\n",
+ "i_output = 1\n",
+ "for i_time in range(\n",
+ " 0,\n",
+ " int(RUN_PART_ARGS['run_part_opt'].t_max / RUN_PART_ARGS['run_part_opt'].del_t) + 1\n",
+ "):\n",
+ " if i_time != 0:\n",
+ " (last[\"output_time\"], last[\"progress_time\"], i_output) = ppmc.run_part_timestep(\n",
+ " *[RUN_PART_ARGS[key] for key in (\n",
+ " 'scenario', 'env_state', 'aero_data',\n",
+ " 'aero_state', 'gas_data', 'gas_state',\n",
+ " 'run_part_opt', 'camp_core', 'photolysis'\n",
+ " )],\n",
+ " i_time, T_INITIAL, last[\"output_time\"], last[\"progress_time\"], i_output,\n",
+ " )\n",
+ " for obj, keys in (\n",
+ " (RUN_PART_ARGS['aero_state'], ('total_num_conc', 'total_mass_conc')),\n",
+ " (RUN_PART_ARGS['gas_state'], ('mix_rats',)),\n",
+ " (RUN_PART_ARGS['env_state'], ('elapsed_time', 'height', 'temp', 'rh')),\n",
+ " ):\n",
+ " for key in keys:\n",
+ " OUTPUT[key].append(getattr(obj, key))\n",
+ " if np.mod(i_time * RUN_PART_ARGS['run_part_opt'].del_t, 3600.0) == 0:\n",
+ " OUTPUT['dists'].append(ppmc.histogram_1d(\n",
+ " DIAM_GRID,\n",
+ " RUN_PART_ARGS['aero_state'].dry_diameters,\n",
+ " RUN_PART_ARGS['aero_state'].num_concs\n",
+ " ))\n",
+ "for key in OUTPUT:\n",
+ " OUTPUT[key] = np.asarray(OUTPUT[key])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "id": "82c0cebb",
+ "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": 20,
+ "id": "47474cd3",
+ "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": "1d36a40d0404460cabbc36821e3a0a2f",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./tmpagzby5sf.pdf
\"), HTML(value…"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.plot(OUTPUT['elapsed_time'], OUTPUT['temp'], 'r')\n",
+ "plt.ylabel('Temperature (K)', color='r')\n",
+ "plt.ylim([275, 300])\n",
+ "plt.xticks(np.linspace(0, OUTPUT['elapsed_time'][-1], 5))\n",
+ "plt.xlim([OUTPUT['elapsed_time'][i] for i in (0, -1)])\n",
+ "plt.xlabel('Time (s)')\n",
+ "plt.twinx()\n",
+ "plt.plot(OUTPUT['elapsed_time'], OUTPUT['rh'] * 100, 'g')\n",
+ "plt.ylabel('Relative humidity (%)', color='g')\n",
+ "plt.ylim([0, 100])\n",
+ "show_plot()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "id": "c85a622a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def set_tickmarks(axes, n_yticks):\n",
+ " ylims = axes.get_ylim()\n",
+ " if np.log10(ylims[0]) > 1:\n",
+ " val = -int(np.ceil(np.abs(np.log10(ylims[0])))) + 1\n",
+ " else:\n",
+ " val = int(np.ceil(np.abs(np.log10(ylims[0])))) + 1 \n",
+ " ymin = round(ylims[0] - .1 * ylims[0], val)\n",
+ " ymax = round(ylims[1] + .1 * ylims[1], val)\n",
+ " plt.ylim([ymin, ymax])\n",
+ " plt.yticks(np.linspace(ymin, ymax, n_yticks))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "id": "8e1b89e0",
+ "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": "f8f137632cc94ef3bb77e01d28b2e574",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./tmpze5vo3wo.pdf
\"), HTML(value…"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.plot(OUTPUT['elapsed_time'], OUTPUT['total_mass_conc'], \"b\", label=\"mass conc\")\n",
+ "plt.ylabel(\"Mass concentration (kg m$^{-3}$)\", color='b')\n",
+ "plt.xlabel(\"Time (s)\")\n",
+ "n_ticks = 5\n",
+ "set_tickmarks(plt.gca(), n_ticks)\n",
+ "plt.twinx()\n",
+ "plt.plot(OUTPUT['elapsed_time'], OUTPUT['total_num_conc'], \"g\", label=\"num conc\")\n",
+ "plt.xticks(np.linspace(0, OUTPUT['elapsed_time'][-1], n_ticks))\n",
+ "plt.xlim([OUTPUT['elapsed_time'][i] for i in (0, -1)])\n",
+ "set_tickmarks(plt.gca(), n_ticks)\n",
+ "plt.ylabel(r\"Number concentration ($\\#$ m$^{-3}$)\", color='g')\n",
+ "show_plot()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "id": "386aa7c2",
+ "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": "dc4d39e89e0942268c8739bf581a6e95",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./tmpju1qj1fm.pdf
\"), HTML(value…"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "for i_spec, spec in enumerate([\"O3\"]):\n",
+ " i_spec = RUN_PART_ARGS['gas_data'].spec_by_name(spec)\n",
+ " l, = plt.plot(\n",
+ " OUTPUT['elapsed_time'],\n",
+ " OUTPUT['mix_rats'][:, i_spec],\n",
+ " label=spec\n",
+ " )\n",
+ "plt.xlabel(\"Time (s)\")\n",
+ "plt.ylabel(\"Mixing ratio (ppb)\")\n",
+ "plt.xticks(np.linspace(0, OUTPUT['elapsed_time'][-1], 5))\n",
+ "plt.legend()\n",
+ "show_plot()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 24,
+ "id": "faa1de28",
+ "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": "8d1a210b526145e79fb582b7d3409185",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./tmp6_g1f0zl.pdf
\"), HTML(value…"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "for hour in (0, 6, 12, 24):\n",
+ " plt.plot(DIAM_GRID.centers, OUTPUT['dists'][hour], label=f'$t = {hour}$ h')\n",
+ "plt.xscale(\"log\")\n",
+ "plt.xlabel(\"Dry diameter (m)\")\n",
+ "plt.ylabel(r\"Number concentration $n(D)$ ($\\#$ m$^{-3}$)\")\n",
+ "plt.ylim(bottom=0)\n",
+ "plt.legend()\n",
+ "plt.xlim([DIAM_GRID.edges[i] for i in (0, -1)])\n",
+ "show_plot()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3ddc93ff-798a-48e6-b60e-3c50dfe5f906",
+ "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.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/src/aero_data.F90 b/src/aero_data.F90
index edf84f95..3baae8a8 100644
--- a/src/aero_data.F90
+++ b/src/aero_data.F90
@@ -36,6 +36,19 @@ subroutine f_aero_data_from_json(ptr_c) bind(C)
call spec_file_read_aero_data(file, ptr_f)
end subroutine
+ subroutine f_aero_data_from_camp(ptr_c, camp_core_ptr_c) bind(C)
+ type(aero_data_t), pointer :: ptr_f => null()
+ type(c_ptr), intent(in) :: ptr_c
+ type(c_ptr), intent(in) :: camp_core_ptr_c
+ type(camp_core_t), pointer :: camp_core_ptr_f => null()
+
+ call c_f_pointer(ptr_c, ptr_f)
+ call c_f_pointer(camp_core_ptr_c, camp_core_ptr_f)
+
+ call aero_data_initialize(ptr_f, camp_core_ptr_f)
+
+ end subroutine
+
subroutine f_aero_data_spec_by_name(ptr_c, value, name_data, name_size) bind(C)
type(aero_data_t), pointer :: ptr_f => null()
type(c_ptr), intent(in) :: ptr_c
diff --git a/src/aero_data.hpp b/src/aero_data.hpp
index 4143959f..de880a9a 100644
--- a/src/aero_data.hpp
+++ b/src/aero_data.hpp
@@ -8,12 +8,14 @@
#include "pmc_resource.hpp"
#include "json_resource.hpp"
+#include "camp_core.hpp"
#include "pybind11/stl.h"
#include "aero_data_parameters.hpp"
extern "C" void f_aero_data_ctor(void *ptr) noexcept;
extern "C" void f_aero_data_dtor(void *ptr) noexcept;
extern "C" void f_aero_data_from_json(const void *ptr) noexcept;
+extern "C" void f_aero_data_from_camp(const void *ptr, const void *camp_core_ptr) noexcept;
extern "C" void f_aero_data_spec_by_name(const void *ptr, int *value, const char *name_data, const int *name_size) noexcept;
extern "C" void f_aero_data_len(const void *ptr, int *len) noexcept;
extern "C" void f_aero_data_n_source(const void *ptr, int *len) noexcept;
@@ -38,6 +40,12 @@ extern "C" void f_aero_data_spec_name_by_index(const void *ptr, const int *i_spe
struct AeroData {
PMCResource ptr;
+ AeroData(const nlohmann::json &json, const CampCore &camp_core) :
+ ptr(f_aero_data_ctor, f_aero_data_dtor)
+ {
+ f_aero_data_from_camp(this->ptr.f_arg(), camp_core.ptr.f_arg());
+ }
+
AeroData(const nlohmann::json &json) :
ptr(f_aero_data_ctor, f_aero_data_dtor)
{
diff --git a/src/aero_state.F90 b/src/aero_state.F90
index 231730e4..68e50eea 100644
--- a/src/aero_state.F90
+++ b/src/aero_state.F90
@@ -638,4 +638,18 @@ subroutine f_aero_state_sample_particles(ptr_c, aero_state_to_ptr_c, &
end subroutine
+ subroutine f_aero_state_initialize(ptr_c, aero_data_ptr_c, camp_core_ptr_c) bind(C)
+ type(c_ptr) :: ptr_c, aero_data_ptr_c, camp_core_ptr_c
+ type(aero_state_t), pointer :: ptr_f => null()
+ type(aero_data_t), pointer :: aero_data_ptr_f => null()
+ type(camp_core_t), pointer :: camp_core_ptr_f => null()
+
+ call c_f_pointer(ptr_c, ptr_f)
+ call c_f_pointer(aero_data_ptr_c, aero_data_ptr_f)
+ call c_f_pointer(camp_core_ptr_c, camp_core_ptr_f)
+
+ call aero_state_initialize(ptr_f, aero_data_ptr_f, camp_core_ptr_f)
+
+ end subroutine
+
end module
diff --git a/src/aero_state.hpp b/src/aero_state.hpp
index 53f63606..2ba0ce98 100644
--- a/src/aero_state.hpp
+++ b/src/aero_state.hpp
@@ -9,6 +9,7 @@
#include "pmc_resource.hpp"
#include "aero_data.hpp"
#include "aero_dist.hpp"
+#include "camp_core.hpp"
#include "aero_particle.hpp"
#include "env_state.hpp"
#include "bin_grid.hpp"
@@ -208,6 +209,12 @@ extern "C" void f_aero_state_sample_particles(
const double *sample_prob
) noexcept;
+extern "C" void f_aero_state_initialize(
+ void *ptr_c,
+ const void *aero_data_ptr_c,
+ const void *camp_core_ptr_c
+) noexcept;
+
template
auto pointer_vec_magic(arr_t &data_vec, const arg_t &arg) {
std::vector pointer_vec(data_vec.size());
@@ -269,6 +276,44 @@ struct AeroState {
{
}
+ AeroState(
+ std::shared_ptr aero_data,
+ const double &n_part,
+ const bpstd::string_view &weight,
+ const CampCore &camp_core
+ ):
+ ptr(f_aero_state_ctor, f_aero_state_dtor),
+ aero_data(aero_data)
+ {
+ static const std::map weight_c{
+ //{"none", '-'},
+ {"flat", 'f'},
+ {"flat_source", 'F'},
+ //{"power", 'p'},
+ //{"power_source", 'P'},
+ {"nummass", 'n'},
+ {"nummass_source", 'N'},
+ };
+
+ if (weight_c.find(weight) == weight_c.end()) {
+ std::ostringstream msg;
+ msg << "unknown weighting scheme '" << weight << "', valid options are: ";
+ auto index = 0;
+ for (auto const& pair: weight_c)
+ msg << (!index++ ? "" : ", ") << pair.first;
+ throw std::runtime_error(msg.str());
+ }
+
+ f_aero_state_init(
+ ptr.f_arg(),
+ aero_data->ptr.f_arg(),
+ &n_part,
+ &weight_c.at(weight)
+ );
+ f_aero_state_initialize(ptr.f_arg_non_const(), aero_data->ptr.f_arg(),
+ camp_core.ptr.f_arg());
+ }
+
static std::size_t __len__(const AeroState &self) {
int len;
f_aero_state_len(
diff --git a/src/camp_core.F90 b/src/camp_core.F90
index 58d2f30c..870324f6 100644
--- a/src/camp_core.F90
+++ b/src/camp_core.F90
@@ -20,6 +20,24 @@ subroutine f_camp_core_ctor(ptr_c) bind(C)
ptr_c = c_loc(ptr_f)
end subroutine
+ subroutine f_camp_core_initialize(ptr_c, prefix_data, prefix_size) bind(C)
+ type(camp_core_t), pointer :: ptr_f => null()
+ type(c_ptr), intent(out) :: ptr_c
+
+ character(kind=c_char), dimension(*), intent(in) :: prefix_data
+ integer(c_int), intent(in) :: prefix_size
+ character(len=prefix_size) :: prefix
+ integer :: i
+
+ do i=1, prefix_size
+ prefix(i:i) = prefix_data(i)
+ end do
+
+ ptr_f => camp_core_t(prefix)
+ call ptr_f%initialize()
+ ptr_c = c_loc(ptr_f)
+ end subroutine
+
subroutine f_camp_core_dtor(ptr_c) bind(C)
type(camp_core_t), pointer :: ptr_f => null()
type(c_ptr), intent(in) :: ptr_c
diff --git a/src/camp_core.hpp b/src/camp_core.hpp
index 53cd0186..bfd4ea28 100644
--- a/src/camp_core.hpp
+++ b/src/camp_core.hpp
@@ -11,6 +11,8 @@
extern "C" void f_camp_core_ctor(void *ptr) noexcept;
extern "C" void f_camp_core_dtor(void *ptr) noexcept;
+extern "C" void f_camp_core_initialize(void *ptr, const char *prefix,
+ const int *prefix_size) noexcept;
struct CampCore {
PMCResource ptr;
@@ -19,4 +21,11 @@ struct CampCore {
ptr(f_camp_core_ctor, f_camp_core_dtor)
{
}
+
+ CampCore(const std::string &config_name) :
+ ptr(f_camp_core_ctor, f_camp_core_dtor)
+ {
+ const int config_name_size = config_name.size();
+ f_camp_core_initialize(ptr.f_arg_non_const(), config_name.c_str(), &config_name_size);
+ }
};
diff --git a/src/gas_data.F90 b/src/gas_data.F90
index c0dab53e..d106c01e 100644
--- a/src/gas_data.F90
+++ b/src/gas_data.F90
@@ -46,6 +46,16 @@ subroutine f_gas_data_from_json(ptr_c) bind(C)
call spec_file_read_gas_data(nofile, ptr_f)
end subroutine
+ subroutine f_gas_data_from_camp(ptr_c, camp_core_ptr_c) bind(C)
+ type(gas_data_t), pointer :: ptr_f => null()
+ type(c_ptr), intent(in) :: ptr_c, camp_core_ptr_c
+ type(camp_core_t), pointer :: camp_core_ptr_f => null()
+
+ call c_f_pointer(ptr_c, ptr_f)
+ call c_f_pointer(camp_core_ptr_c, camp_core_ptr_f)
+ call gas_data_initialize(ptr_f, camp_core_ptr_f)
+ end subroutine
+
subroutine f_gas_data_spec_by_name(ptr_c, value, name_data, name_size) &
bind(C)
type(gas_data_t), pointer :: ptr_f => null()
diff --git a/src/gas_data.hpp b/src/gas_data.hpp
index 7bcd0103..d101afbe 100644
--- a/src/gas_data.hpp
+++ b/src/gas_data.hpp
@@ -9,6 +9,7 @@
#include "json_resource.hpp"
#include "pmc_resource.hpp"
#include "gas_data_parameters.hpp"
+#include "camp_core.hpp"
#include "pybind11/stl.h"
#include "pybind11_json/pybind11_json.hpp"
@@ -16,6 +17,7 @@ extern "C" void f_gas_data_ctor(void *ptr) noexcept;
extern "C" void f_gas_data_dtor(void *ptr) noexcept;
extern "C" void f_gas_data_len(const void *ptr, int *len) noexcept;
extern "C" void f_gas_data_from_json(const void *ptr) noexcept;
+extern "C" void f_gas_data_from_camp(const void *ptr, const void *camp_core_ptr) noexcept;
extern "C" void f_gas_data_to_json(const void *ptr) noexcept;
extern "C" void f_gas_data_spec_by_name(const void *ptr, int *value, const char *name_data,
const int *name_size) noexcept;
@@ -26,6 +28,12 @@ struct GasData {
PMCResource ptr;
const nlohmann::json json;
+ GasData(const CampCore &CampCore) :
+ ptr(f_gas_data_ctor, f_gas_data_dtor)
+ {
+ f_gas_data_from_camp(this->ptr.f_arg(), CampCore.ptr.f_arg());
+ }
+
GasData(const pybind11::tuple &tpl) :
ptr(f_gas_data_ctor, f_gas_data_dtor),
json(tpl)
diff --git a/src/photolysis.F90 b/src/photolysis.F90
index c628b524..89526931 100644
--- a/src/photolysis.F90
+++ b/src/photolysis.F90
@@ -6,6 +6,7 @@
module PyPartMC_photolysis
use iso_c_binding
+ use camp_camp_core
use pmc_photolysis
implicit none
@@ -19,6 +20,20 @@ subroutine f_photolysis_ctor(ptr_c) bind(C)
ptr_c = c_loc(ptr_f)
end subroutine
+ subroutine f_photolysis_create(ptr_c, camp_core_ptr_c) bind(C)
+ type(photolysis_t), pointer :: ptr_f => null()
+ type(camp_core_t), pointer :: camp_core_ptr_f => null()
+ type(c_ptr), intent(inout) :: ptr_c
+ type(c_ptr), intent(in) :: camp_core_ptr_c
+
+ call c_f_pointer(ptr_c, ptr_f)
+ call c_f_pointer(camp_core_ptr_c, camp_core_ptr_f)
+
+ ptr_f => photolysis_t(camp_core_ptr_f)
+
+ ptr_c = c_loc(ptr_f)
+ end subroutine
+
subroutine f_photolysis_dtor(ptr_c) bind(C)
type(photolysis_t), pointer :: ptr_f => null()
type(c_ptr), intent(in) :: ptr_c
diff --git a/src/photolysis.hpp b/src/photolysis.hpp
index f0095762..54fa830e 100644
--- a/src/photolysis.hpp
+++ b/src/photolysis.hpp
@@ -8,10 +8,11 @@
#include "json_resource.hpp"
#include "pmc_resource.hpp"
+#include "camp_core.hpp"
extern "C" void f_photolysis_ctor(void *ptr) noexcept;
extern "C" void f_photolysis_dtor(void *ptr) noexcept;
-
+extern "C" void f_photolysis_create(const void *ptr, const void *camp_core) noexcept;
struct Photolysis {
PMCResource ptr;
@@ -19,4 +20,10 @@ struct Photolysis {
ptr(f_photolysis_ctor, f_photolysis_dtor)
{
}
+
+ Photolysis(const CampCore &camp_core) :
+ ptr(f_photolysis_ctor, f_photolysis_dtor)
+ {
+ f_photolysis_create(this->ptr.f_arg(), camp_core.ptr.f_arg());
+ }
};
diff --git a/src/pypartmc.cpp b/src/pypartmc.cpp
index 99c6f9af..d01adfd5 100644
--- a/src/pypartmc.cpp
+++ b/src/pypartmc.cpp
@@ -104,6 +104,7 @@ PYBIND11_MODULE(_PyPartMC, m) {
same, but without the \c _a suffix.
)pbdoc"
)
+ .def(py::init())
.def(py::init())
.def("spec_by_name", AeroData::spec_by_name,
"Returns the number of the species in AeroData with the given name")
@@ -229,6 +230,8 @@ PYBIND11_MODULE(_PyPartMC, m) {
)pbdoc"
)
.def(py::init, const double, const std::string>())
+ .def(py::init, const double, const std::string,
+ const CampCore&>())
.def("__len__", AeroState::__len__,
"returns current number of particles")
.def_property_readonly("total_num_conc", AeroState::total_num_conc,
@@ -307,6 +310,7 @@ PYBIND11_MODULE(_PyPartMC, m) {
is gas_state%%mix_rat(i).
)pbdoc"
)
+ .def(py::init())
.def(py::init())
.def("__len__", GasData::__len__,
"returns number of gas species")
@@ -357,6 +361,7 @@ PYBIND11_MODULE(_PyPartMC, m) {
)pbdoc"
)
.def(py::init<>())
+ .def(py::init())
;
py::class_(m,
@@ -366,6 +371,7 @@ PYBIND11_MODULE(_PyPartMC, m) {
)pbdoc"
)
.def(py::init<>())
+ .def(py::init())
;
py::class_(m,
diff --git a/src/run_part.F90 b/src/run_part.F90
index a971aa08..aaae1e71 100644
--- a/src/run_part.F90
+++ b/src/run_part.F90
@@ -8,6 +8,8 @@ module PyPartMC_run_part
use iso_c_binding
use pmc_run_part
+ use camp_camp_core
+ use pmc_photolysis
implicit none
@@ -61,6 +63,10 @@ subroutine f_run_part( &
call c_f_pointer(camp_core_ptr_c, camp_core_ptr_f)
call c_f_pointer(photolysis_ptr_c, photolysis_ptr_f)
+ if (run_part_opt_ptr_f%do_camp_chem) then
+ call camp_core_ptr_f%solver_initialize()
+ end if
+
call run_part( &
scenario_ptr_f, &
env_state_ptr_f, &
@@ -148,6 +154,9 @@ subroutine f_run_part_timestep( &
if (env_state_ptr_f%elapsed_time < run_part_opt_ptr_f%del_t) then
call mosaic_init(env_state_ptr_f, aero_data_ptr_f, run_part_opt_ptr_f%del_t, &
run_part_opt_ptr_f%do_optical)
+ if (run_part_opt_ptr_f%do_camp_chem) then
+ call camp_core_ptr_f%solver_initialize()
+ end if
if (run_part_opt_ptr_f%t_output > 0) then
call output_state(run_part_opt_ptr_f%output_prefix, &
run_part_opt_ptr_f%output_type, aero_data_ptr_f, aero_state_ptr_f, gas_data_ptr_f, &
@@ -156,6 +165,7 @@ subroutine f_run_part_timestep( &
run_part_opt_ptr_f%do_optical, run_part_opt_ptr_f%uuid)
end if
end if
+
call run_part_timestep(scenario_ptr_f, env_state_ptr_f, aero_data_ptr_f, aero_state_ptr_f, &
gas_data_ptr_f, gas_state_ptr_f, run_part_opt_ptr_f, camp_core_ptr_f, photolysis_ptr_f, &
i_time, t_start, last_output_time, &
@@ -241,6 +251,9 @@ subroutine f_run_part_timeblock( &
if (env_state_ptr_f%elapsed_time < run_part_opt_ptr_f%del_t) then
call mosaic_init(env_state_ptr_f, aero_data_ptr_f, run_part_opt_ptr_f%del_t, &
run_part_opt_ptr_f%do_optical)
+ if (run_part_opt_ptr_f%do_camp_chem) then
+ call camp_core_ptr_f%solver_initialize()
+ end if
if (run_part_opt_ptr_f%t_output > 0) then
call output_state(run_part_opt_ptr_f%output_prefix, &
run_part_opt_ptr_f%output_type, aero_data_ptr_f, aero_state_ptr_f, gas_data_ptr_f, &
diff --git a/tests/test_photolysis.py b/tests/test_photolysis.py
new file mode 100644
index 00000000..841364b9
--- /dev/null
+++ b/tests/test_photolysis.py
@@ -0,0 +1,11 @@
+import PyPartMC as ppmc
+
+
+class TestPhotolysis:
+ @staticmethod
+ def test_ctor_without_camp():
+ _ = ppmc.Photolysis()
+
+ @staticmethod
+ def test_ctor_with_camp():
+ _ = ppmc.Photolysis(ppmc.CampCore())