diff --git a/examples/additive_coag_comparison.ipynb b/examples/additive_coag_comparison.ipynb
index b29d201e..481db287 100644
--- a/examples/additive_coag_comparison.ipynb
+++ b/examples/additive_coag_comparison.ipynb
@@ -52,13 +52,7 @@
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from open_atmos_jupyter_utils import show_plot\n",
- "from open_atmos_jupyter_utils.show_anim import show_anim \n",
- "import PySDM as pysdm\n",
- "from PySDM import environments\n",
- "from PySDM.initialisation import spectra\n",
- "from PySDM.dynamics.collisions import collision_kernels\n",
- "from PySDM import products\n",
- "import PyPartMC as ppmc"
+ "from open_atmos_jupyter_utils.show_anim import show_anim \n"
]
},
{
@@ -68,4190 +62,471 @@
"metadata": {},
"outputs": [],
"source": [
- "N_PART = 2**16\n",
- "VOLUME_M3 = 1e6\n",
- "DT_SEC = 1.\n",
- "T_MAX_SEC = 3600\n",
- "NUM_CONC_PER_M3 = 2**23\n",
- "DIAM_AT_MEAN_VOL_M = 2*30.531e-6\n",
- "ADDITIVE_KERNEL_COEFF = 1000\n",
- "AMBIENT_T_K = 288\n",
- "AMBIENT_P_Pa = 1e5\n",
- "AMBIENT_RH = 0.999\n",
+ "from collections import namedtuple\n",
"\n",
- "# Plotting parameters\n",
- "N_BINS = 128\n",
- "N_BIN_EDGES = N_BINS + 1\n",
- "R_BINS_MIN = 1e-5\n",
- "R_BINS_MAX = 5e-3\n",
- "PLOT_TIME_STEP = 120 # seconds\n",
- "BINNING_METHOD = 'Mass Density (g/m^3 dlnr)' # 'Number Concentration (m^-3)' or 'Mass Density (g/m^3 dlnr)'\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 5,
- "id": "7f986747",
- "metadata": {},
- "outputs": [
- {
- "name": "stderr",
- "output_type": "stream",
- "text": [
- "/Users/emmaware/Library/Python/3.9/lib/python/site-packages/PySDM/backends/numba.py:46: UserWarning: Disabling Numba threading due to ARM64 CPU (atomics do not work yet)\n",
- " warnings.warn(\n"
- ]
- }
- ],
- "source": [
- "radius_bins_edges=np.logspace(np.log10(R_BINS_MIN), np.log10(R_BINS_MAX), num=N_BIN_EDGES, endpoint=True)\n",
- "\n",
- "\n",
- "builder = pysdm.Builder(\n",
- " N_PART,\n",
- " backend=pysdm.backends.CPU(),\n",
- " environment=environments.Box(dt=DT_SEC * pysdm.physics.si.s, \n",
- " dv=VOLUME_M3 * pysdm.physics.si.m**3)\n",
- ")\n",
- "trivia = builder.formulae.trivia\n",
- "spectrum = spectra.Exponential(\n",
- " norm_factor=NUM_CONC_PER_M3 *VOLUME_M3 / pysdm.physics.si.m**3, #TODO: check\n",
- " scale=trivia.volume(radius=DIAM_AT_MEAN_VOL_M / 2 * pysdm.physics.si.m)\n",
- ")\n",
- "builder.add_dynamic(pysdm.dynamics.Coalescence(\n",
- " collision_kernel=collision_kernels.Golovin(b=ADDITIVE_KERNEL_COEFF)\n",
- "))\n",
- "if BINNING_METHOD == 'Number Concentration (m^-3)':\n",
- " binning = (products.ParticleSizeSpectrumPerVolume(\n",
- " radius_bins_edges=radius_bins_edges,\n",
- " name='spectrum',),)\n",
- "elif BINNING_METHOD == 'Mass Density (g/m^3 dlnr)':\n",
- " binning = (products.ParticleVolumeVersusRadiusLogarithmSpectrum(\n",
- " radius_bins_edges=radius_bins_edges,\n",
- " name='spectrum',),)\n",
- "\n",
- "particulator = builder.build(\n",
- " attributes=builder.particulator.environment.init_attributes(\n",
- " spectral_discretisation=pysdm.initialisation.sampling.spectral_sampling.Logarithmic(\n",
- " spectrum,\n",
- " )\n",
- " ),\n",
- " products=binning\n",
- ")\n",
- "# print(np.amax(particulator.attributes[\"multiplicity\"].to_ndarray()))\n",
- "particulator.environment[\"rhod\"] = 1. # TODO!\n",
- "pysdm_output = {}\n",
- "for step in range(int(T_MAX_SEC // PLOT_TIME_STEP)+1):\n",
- " if step != 0:\n",
- " particulator.run(int(PLOT_TIME_STEP // DT_SEC))\n",
- " if BINNING_METHOD == 'Number Concentration (m^-3)':\n",
- " pysdm_output[step] = particulator.products['spectrum'].get()\n",
- " #next two lines convert from dr to dlnr\n",
- " pysdm_output[step] *= (radius_bins_edges[1:]-radius_bins_edges[:-1])\n",
- " pysdm_output[step] /= (np.log(radius_bins_edges[1:]) - np.log(radius_bins_edges[:-1]))\n",
- " elif BINNING_METHOD == 'Mass Density (g/m^3 dlnr)': \n",
- " pysdm_output[step] = particulator.products['spectrum'].get()[0]\n",
- " pysdm_output[step][:] *= pysdm.physics.constants_defaults.rho_w * pysdm.physics.si.kilogram / pysdm.physics.si.metre**3 \n",
- " pysdm_output[step][:] /= pysdm.physics.si.g"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 6,
- "id": "6c9e2cbd",
- "metadata": {},
- "outputs": [],
- "source": [
- "collision_kernel=collision_kernels.Golovin(b=ADDITIVE_KERNEL_COEFF)\n",
- "\n",
- "def analytical_solution_golovin(t):\n",
- " X0 = trivia.volume(DIAM_AT_MEAN_VOL_M / 2)\n",
- " volume_bins_edges = trivia.volume(radius_bins_edges)\n",
- " dm = np.diff(volume_bins_edges)\n",
- " dr = np.diff(radius_bins_edges)\n",
- " pdf_m_x = volume_bins_edges[:-1] + dm / 2\n",
- " pdf_r_x = radius_bins_edges[:-1] + dr / 2\n",
- " pdf_m_y = NUM_CONC_PER_M3*VOLUME_M3*collision_kernel.analytic_solution(x = pdf_m_x,t = t,x_0=X0,N_0=NUM_CONC_PER_M3)\n",
- " pdf_r_y = pdf_m_y * dm / dr * pdf_r_x\n",
- " y_true = (\n",
- " pdf_r_y\n",
- " * trivia.volume(radius=pdf_r_x)\n",
- " * pysdm.physics.constants_defaults.rho_w\n",
- " / VOLUME_M3*pysdm.physics.si.metres**3\n",
- " * pysdm.physics.si.kilograms\n",
- " / pysdm.physics.si.grams\n",
- " )\n",
- "\n",
- " return y_true\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 7,
- "id": "6c8b3b0c",
- "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": "29bbee5e1a65441ea3671b41ce9ce5e2",
- "version_major": 2,
- "version_minor": 0
- },
- "text/plain": [
- "HBox(children=(HTML(value=\"./tmpbn0w2ohc.pdf
\"), HTML(value…"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "for step, output in pysdm_output.items():\n",
- " if step*PLOT_TIME_STEP%1200!=0:\n",
- " continue\n",
- " plt.plot(\n",
- " ((radius_bins_edges[1:] + radius_bins_edges[:-1]) / 2),\n",
- " output,\n",
- " label=f\"t: {step*PLOT_TIME_STEP} sec\",\n",
- " marker='x'\n",
- " )\n",
- "plt.xscale(\"log\")\n",
- "\n",
- "if BINNING_METHOD == 'Number Concentration (m^-3)':\n",
- " plt.yscale(\"log\")\n",
- " plt.ylim([1e2,5e10])\n",
- "plt.ylabel(BINNING_METHOD)\n",
- "plt.xlim([R_BINS_MIN,R_BINS_MAX])\n",
- "plt.xlabel(\"particle radius [m]\")\n",
- "plt.legend()\n",
- "\n",
- "show_plot()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 8,
- "id": "df8d1a84",
- "metadata": {},
- "outputs": [],
- "source": [
- "gas_data = ppmc.GasData((\"H2SO4\",\"HNO3\",\"HCl\",\"NH3\",\"NO\",\"NO2\"))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 9,
- "id": "a2d7bad8",
- "metadata": {},
- "outputs": [],
- "source": [
- "env_state = ppmc.EnvState(\n",
- " {\n",
- " \"rel_humidity\": AMBIENT_RH,\n",
- " \"latitude\": 40,\n",
- " \"longitude\": 0,\n",
- " \"altitude\": 0 * ppmc.si.m,\n",
- " \"start_time\": 0 * ppmc.si.s,\n",
- " \"start_day\": 1,\n",
- " \"additive_kernel_coefficient\": ADDITIVE_KERNEL_COEFF,\n",
- " }\n",
- ")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 10,
- "id": "815283a0",
- "metadata": {},
- "outputs": [],
- "source": [
- "aero_data = ppmc.AeroData(\n",
- " (\n",
- " {\"H2O\": [1000 * ppmc.si.kg / ppmc.si.m**3, 0, 18.0 * ppmc.si.g / ppmc.si.mol, 0.00]},\n",
- " )\n",
- ")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 11,
- "id": "c822823d",
- "metadata": {},
- "outputs": [],
- "source": [
- "gas_state = ppmc.GasState(gas_data)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 12,
- "id": "58706963",
- "metadata": {},
- "outputs": [],
- "source": [
- "times = [0 * ppmc.si.s]\n",
- "back_gas = [{\"time\": times},\n",
- " {\"rate\": [0 / ppmc.si.s]},\n",
- " {\"NO\": [0.0E+00]},\n",
- " ]\n",
- "\n",
- "gas_emit_times = [0]\n",
- "\n",
- "gas_emit_rates = np.ones(len(gas_emit_times))\n",
- "\n",
- "emit_gas = [\n",
- " {\"time\": gas_emit_times},\n",
- " {\"rate\": list(gas_emit_rates)},\n",
- " {\"NO\": [0]},\n",
- "]\n",
- "\n",
- "AERO_DIST_BACKGROUND = {\n",
- " \"back_small\": {\n",
- " \"mass_frac\": [{\"H2O\": [1]}],\n",
- " \"diam_type\": \"geometric\",\n",
- " \"mode_type\": \"log_normal\",\n",
- " \"num_conc\": 0 / ppmc.si.m**3,\n",
- " \"geom_mean_diam\": 0.02 * ppmc.si.um,\n",
- " \"log10_geom_std_dev\": 0.161,\n",
- " },\n",
- "}\n",
- "\n",
- "AERO_DIST_EMIT = {\n",
- " \"gasoline\": {\n",
- " \"mass_frac\": [{\"H2O\": [1.0]}],\n",
- " \"diam_type\": \"geometric\",\n",
- " \"mode_type\": \"log_normal\",\n",
- " \"num_conc\": 0 / ppmc.si.m**3,\n",
- " \"geom_mean_diam\": 5e-8 * ppmc.si.m,\n",
- " \"log10_geom_std_dev\": 0.24,\n",
- " },\n",
- "}"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 13,
- "id": "c6a96b7d",
- "metadata": {},
- "outputs": [],
- "source": [
- "time_timeseries = [0]\n",
- "pressure_timeseries = [AMBIENT_P_Pa * ppmc.si.Pa]\n",
- "temp_timeseries = [AMBIENT_T_K * ppmc.si.K]\n",
- "height_timeseries = [1000]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 14,
- "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 * ppmc.si.s]},\n",
- " {\"rate\": [0 / ppmc.si.s]},\n",
- " {\"dist\": [[AERO_DIST_EMIT]]},\n",
- " ],\n",
- " \"aero_background\": [\n",
- " {\"time\": [0 * ppmc.si.s]},\n",
- " {\"rate\": [0 / ppmc.si.s]},\n",
- " {\"dist\": [[AERO_DIST_BACKGROUND]]},\n",
- " ],\n",
- " \"loss_function\": \"none\",\n",
- " },\n",
- ")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 15,
- "id": "6722ba83",
- "metadata": {},
- "outputs": [],
- "source": [
- "T_INITIAL = 0.0\n",
- "scenario.init_env_state(env_state, T_INITIAL)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 16,
- "id": "9781ca2f",
- "metadata": {},
- "outputs": [],
- "source": [
- "AERO_DIST_INIT = [\n",
- " {\n",
- " \"init1\": {\n",
- " \"mass_frac\": [{\"H2O\": [1]}],\n",
- " \"diam_type\": \"geometric\",\n",
- " \"mode_type\": \"exp\",\n",
- " \"num_conc\": NUM_CONC_PER_M3 / ppmc.si.m**3,\n",
- " \"diam_at_mean_vol\": DIAM_AT_MEAN_VOL_M * ppmc.si.m,\n",
- " },\n",
- " }\n",
- "]\n",
- "\n",
- "aero_dist_init = ppmc.AeroDist(aero_data, AERO_DIST_INIT)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 17,
- "id": "d8d9c8fd",
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "65338"
- ]
- },
- "execution_count": 17,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
- "source": [
- "run_part_opt = ppmc.RunPartOpt(\n",
- " {\n",
- " \"output_prefix\": 'additive',\n",
- " \"do_coagulation\": True,\n",
- " \"coag_kernel\": \"additive\",\n",
- " \"t_max\": T_MAX_SEC * ppmc.si.s,\n",
- " \"del_t\": DT_SEC * ppmc.si.s,\n",
- " }\n",
- ")\n",
- "\n",
- "aero_state = ppmc.AeroState(aero_data, N_PART, \"flat\")\n",
- "aero_state.dist_sample(aero_dist_init)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 18,
- "id": "5c8aef2b",
- "metadata": {},
- "outputs": [],
- "source": [
- "camp_core = ppmc.CampCore()\n",
- "photolysis = ppmc.Photolysis()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 19,
- "id": "0f7c29fe",
- "metadata": {},
- "outputs": [],
- "source": [
- "N_STEPS = int(run_part_opt.t_max / run_part_opt.del_t)\n",
- "num_conc = np.zeros(N_STEPS + 1)\n",
- "num_conc[0] = aero_state.total_num_conc\n",
- "mass_conc = np.zeros(N_STEPS + 1)\n",
- "mass_conc[0] = aero_state.total_mass_conc\n",
- "time = np.zeros(N_STEPS + 1)\n",
- "gas_mix_rat = np.zeros((N_STEPS + 1, gas_state.n_spec))\n",
- "gas_mix_rat[0, :] = gas_state.mix_rats\n",
- "\n",
- "height = np.zeros((N_STEPS + 1))\n",
- "temperature = np.zeros((N_STEPS + 1))\n",
- "rh = np.zeros((N_STEPS + 1))\n",
- "\n",
- "height[0] = env_state.height\n",
- "temperature[0] = env_state.temp\n",
- "rh[0] = env_state.rh\n",
- "\n",
- "diam_grid = ppmc.BinGrid(N_BINS, \"log\", (R_BINS_MIN)*2, (R_BINS_MAX)*2)\n",
- "rad_grid = ppmc.BinGrid(N_BINS, \"log\", R_BINS_MIN, R_BINS_MAX)\n",
- "dists = []\n",
- "diameters = aero_state.diameters()\n",
- "num_concs = aero_state.num_concs\n",
- "masses = aero_state.masses()\n",
- "mass_concs = np.array(num_concs) * np.array(masses)\n",
- "dists_mass = []\n",
- "dists_dmlnr = []\n",
- "dists.append(ppmc.histogram_1d(rad_grid, np.array(diameters)/2, num_concs))\n",
- "dists_mass.append(ppmc.histogram_1d(rad_grid, np.array(diameters)/2, mass_concs / ppmc.si.g))\n",
- "\n",
- "for i_time in range(1,N_STEPS + 1):\n",
- " ppmc.run_part_timestep(\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",
- " T_INITIAL,\n",
- " )\n",
- " num_conc[i_time] = aero_state.total_num_conc\n",
- " mass_conc[i_time] = aero_state.total_mass_conc\n",
- " time[i_time] = env_state.elapsed_time\n",
- " gas_mix_rat[i_time, :] = gas_state.mix_rats\n",
- " height[i_time] = env_state.height\n",
- " temperature[i_time] = env_state.temp\n",
- " rh[i_time] = env_state.rh\n",
- " \n",
- " if np.mod(i_time * run_part_opt.del_t, PLOT_TIME_STEP) == 0:\n",
- " diameters = aero_state.diameters()\n",
- " num_concs = aero_state.num_concs\n",
- " masses = aero_state.masses()\n",
- " mass_concs = np.array(num_concs) * np.array(masses)\n",
- " dists.append(ppmc.histogram_1d(rad_grid, np.array(diameters)/2, num_concs))\n",
- " dists_mass.append(ppmc.histogram_1d(rad_grid, np.array(diameters)/2, mass_concs / ppmc.si.g))\n",
- "\n",
- "if BINNING_METHOD == 'Number Concentration (m^-3)':\n",
- " PyPartDists = dists\n",
- "elif BINNING_METHOD == 'Mass Density (g/m^3 dlnr)': \n",
- " PyPartDists = dists_mass\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 20,
- "id": "82c0cebb",
- "metadata": {},
- "outputs": [],
- "source": [
- "plt.rcParams.update({'font.size': 9})\n",
- "plt.rcParams.update({'figure.figsize': (7.08,4.5)})\n",
- "plt.rcParams.update({\"axes.grid\" : True})"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 21,
- "id": "c85a622a",
- "metadata": {},
- "outputs": [],
- "source": [
- "def set_tickmarks(axes, n_ticks):\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_ticks))"
- ]
- },
- {
- "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": "ab198f13c1b148e7b6db5718e4ee4eb2",
- "version_major": 2,
- "version_minor": 0
- },
- "text/plain": [
- "HBox(children=(HTML(value=\"./tmpi0rk2n92.pdf
\"), HTML(value…"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "plt.plot(time, mass_conc, \"b\", label=\"mass conc\")\n",
- "plt.ylabel(\"Mass concentration (kg m$^{-3}$)\", color='b')\n",
- "plt.xlabel(\"Time (s)\")\n",
- "set_tickmarks(plt.gca(), 5)\n",
- "plt.twinx()\n",
- "plt.plot(time, num_conc, \"g\", label=\"num conc\")\n",
- "plt.xticks(np.linspace(0, time[-1], 5))\n",
- "plt.xlim([time[0],time[-1]])\n",
- "plt.ylim([min(num_conc)/2,max(num_conc)*1.01])\n",
- "set_tickmarks(plt.gca(), 5)\n",
- "plt.ylabel(r\"Number concentration ($\\#$ m$^{-3}$)\", color='g')\n",
- "show_plot()"
+ "settings = {\n",
+ " 'N_PART': 2**16,\n",
+ " 'VOLUME_M3': 1e6,\n",
+ " 'DT_SEC': 1.,\n",
+ " 'BINNING_METHOD': 'Mass Density (g/m^3 dlnr)', # 'Number Concentration (m^-3)' or 'Mass Density (g/m^3 dlnr)'\n",
+ " \n",
+ " # Plotting parameters\n",
+ " 'N_BINS': 128,\n",
+ " 'R_BINS_MIN': 1e-5,\n",
+ " 'R_BINS_MAX': 5e-3,\n",
+ " 'NUM_CONC_PER_M3': 2**23,\n",
+ " 'DIAM_AT_MEAN_VOL_M': 2*30.531e-6,\n",
+ " 'ADDITIVE_KERNEL_COEFF': 1500,\n",
+ " 'T_MAX_SEC': 3600,\n",
+ " 'PLOT_TIME_STEP': 120 # seconds\n",
+ "}\n",
+ "settings['N_BIN_EDGES'] = settings[\"N_BINS\"] + 1\n",
+ "\n",
+ "with open('setup.json', 'w', encoding='UTF-8') as f:\n",
+ " json.dump(settings, f)\n",
+ "\n",
+ "settings = namedtuple(\"Settings\", settings.keys())(*settings.values()) # ensure immutable\n",
+ "\n"
]
},
{
"cell_type": "code",
- "execution_count": 23,
- "id": "03ea41f1",
+ "execution_count": 5,
+ "id": "092168b8",
"metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "Examples/testfunctions.jl\n",
- "Examples/DSDvis.jl\n",
- "src/SDfunc/coalescence.jl\n"
- ]
- }
- ],
+ "outputs": [],
"source": [
- "BASE_URL = (\n",
- " # replace main with commit id for previous commit\n",
- " 'https://raw.githubusercontent.com/emmacware/Droplets.jl/refs/heads/dev/' \n",
- ")\n",
- "for path in (\n",
- " 'Examples/testfunctions.jl', 'Examples/DSDvis.jl','src/SDfunc/coalescence.jl'\n",
- "):\n",
- " print(path)\n",
- " with open('Droplets.jl-' + path.replace('/','-'), 'w', encoding='utf-8') as fout:\n",
- " with urllib.request.urlopen(BASE_URL + path) as fin:\n",
- " fout.write(fin.read().decode('utf-8'))"
+ "class AnalyticalSoln:\n",
+ " def __init__(self, settings):\n",
+ " \n",
+ " from PySDM.builder import Builder\n",
+ " from PySDM.physics import si\n",
+ " from PySDM.environments import Box\n",
+ " from PySDM.backends import CPU\n",
+ " from PySDM.dynamics.collisions import collision_kernels \n",
+ " \n",
+ " builder = Builder(\n",
+ " settings.N_PART,\n",
+ " backend=CPU(),\n",
+ " environment=Box(dt=settings.DT_SEC * si.s, dv=settings.VOLUME_M3 * si.m**3)\n",
+ " )\n",
+ "\n",
+ " X0 = builder.formulae.trivia.volume(settings.DIAM_AT_MEAN_VOL_M / 2)\n",
+ " radius_bins_edges = np.logspace(\n",
+ " np.log10(settings.R_BINS_MIN),\n",
+ " np.log10(settings.R_BINS_MAX),\n",
+ " num=settings.N_BIN_EDGES,\n",
+ " endpoint=True\n",
+ " )\n",
+ " volume_bins_edges = builder.formulae.trivia.volume(radius_bins_edges)\n",
+ " dm = np.diff(volume_bins_edges)\n",
+ " dr = np.diff(radius_bins_edges)\n",
+ " pdf_m_x = volume_bins_edges[:-1] + dm / 2\n",
+ " pdf_r_x = radius_bins_edges[:-1] + dr / 2\n",
+ "\n",
+ " self.settings = settings\n",
+ " self.X0 = X0\n",
+ " self.dm = dm\n",
+ " self.dr = dr\n",
+ " self.pdf_m_x = pdf_m_x\n",
+ " self.pdf_r_x = pdf_r_x\n",
+ " \n",
+ " self.builder = builder\n",
+ " self.collision_kernel=collision_kernels.Golovin(b=settings.ADDITIVE_KERNEL_COEFF)\n",
+ " self.times = np.linspace(0, settings.T_MAX_SEC, int(settings.T_MAX_SEC // settings.PLOT_TIME_STEP)+1)\n",
+ "\n",
+ "\n",
+ " def __call__(self): \n",
+ " from PySDM.physics import si, constants_defaults\n",
+ " from PySDM.builder import Builder\n",
+ "\n",
+ " output = {}\n",
+ "\n",
+ " for step in range(int(self.settings.T_MAX_SEC // self.settings.PLOT_TIME_STEP)+1):\n",
+ " t = step*self.settings.PLOT_TIME_STEP or 1e-10 #For 0 time\n",
+ "\n",
+ " pdf_m_y = self.settings.NUM_CONC_PER_M3*self.settings.VOLUME_M3*self.collision_kernel.analytic_solution(x = self.pdf_m_x,t = t,x_0=self.X0,N_0=self.settings.NUM_CONC_PER_M3)\n",
+ " pdf_r_y = pdf_m_y * self.dm / self.dr * self.pdf_r_x\n",
+ " y_true = (\n",
+ " pdf_r_y\n",
+ " * self.builder.formulae.trivia.volume(radius=self.pdf_r_x)\n",
+ " * constants_defaults.rho_w\n",
+ " / self.settings.VOLUME_M3*si.metres**3\n",
+ " * si.kilograms\n",
+ " / si.grams\n",
+ " )\n",
+ " output[step] = y_true\n",
+ "\n",
+ " return output"
]
},
{
"cell_type": "code",
- "execution_count": 24,
- "id": "f830dd13",
+ "execution_count": 6,
+ "id": "29491329",
"metadata": {},
"outputs": [],
"source": [
- "SETUP = {\n",
- " \"n_0\": NUM_CONC_PER_M3,\n",
- " \"Ns\": N_PART,\n",
- " \"V\": VOLUME_M3,\n",
- " \"dt\": DT_SEC,\n",
- " \"t_max\": T_MAX_SEC,\n",
- " \"d_mean\": DIAM_AT_MEAN_VOL_M,\n",
- " \"b\": ADDITIVE_KERNEL_COEFF,\n",
- " \"num_bins\" : N_BINS,\n",
- " \"num_bin_edges\": N_BIN_EDGES,\n",
- " \"radius_bins_log_low\":R_BINS_MIN, \n",
- " \"radius_bins_log_high\":R_BINS_MAX,\n",
- " \"smooth\": False,\n",
- " \"smooth_scope\": 2,\n",
- " \"init_random_seed\": 30,\n",
- " \"coag_threading\": \"Serial()\",\n",
- " \"output_steps\": list(range(0, T_MAX_SEC + PLOT_TIME_STEP, PLOT_TIME_STEP)),\n",
- " \"init_method\": \"init_ξ_const\",\n",
- " \"binning_method\": BINNING_METHOD,\n",
- "}\n",
+ "class RunPySDM:\n",
+ " def __init__(self, settings):\n",
+ " import PySDM\n",
+ " from PySDM.physics import si\n",
"\n",
+ " builder = PySDM.Builder(\n",
+ " settings.N_PART,\n",
+ " backend=PySDM.backends.CPU(),\n",
+ " environment=PySDM.environments.Box(dt=settings.DT_SEC * si.s, dv=settings.VOLUME_M3 * si.m**3)\n",
+ " )\n",
+ " trivia = builder.formulae.trivia\n",
+ " spectrum = PySDM.initialisation.spectra.Exponential(\n",
+ " norm_factor=settings.NUM_CONC_PER_M3 * settings.VOLUME_M3 / si.m**3,\n",
+ " scale=trivia.volume(radius=settings.DIAM_AT_MEAN_VOL_M / 2 * si.m)\n",
+ " )\n",
+ " builder.add_dynamic(PySDM.dynamics.Coalescence(\n",
+ " collision_kernel=PySDM.dynamics.collisions.collision_kernels.Golovin(b=settings.ADDITIVE_KERNEL_COEFF)\n",
+ " ))\n",
"\n",
- "with open('setup.json', 'w', encoding='UTF-8') as f:\n",
- " json.dump(SETUP, f)"
+ " self.radius_bins_edges=np.logspace(\n",
+ " np.log10(settings.R_BINS_MIN),\n",
+ " np.log10(settings.R_BINS_MAX),\n",
+ " num=settings.N_BIN_EDGES,\n",
+ " endpoint=True\n",
+ " )\n",
+ "\n",
+ " self.particulator = builder.build(\n",
+ " attributes=builder.particulator.environment.init_attributes(\n",
+ " spectral_discretisation=PySDM.initialisation.sampling.spectral_sampling.Logarithmic(\n",
+ " spectrum,\n",
+ " )\n",
+ " ),\n",
+ " products={\n",
+ " 'Number Concentration (m^-3)': (PySDM.products.ParticleSizeSpectrumPerVolume(\n",
+ " radius_bins_edges=self.radius_bins_edges,\n",
+ " name='spectrum',),\n",
+ " ),\n",
+ " 'Mass Density (g/m^3 dlnr)': (PySDM.products.ParticleVolumeVersusRadiusLogarithmSpectrum(\n",
+ " radius_bins_edges=self.radius_bins_edges,\n",
+ " name='spectrum',),\n",
+ " )\n",
+ " }[settings.BINNING_METHOD]\n",
+ " )\n",
+ " self.settings = settings\n",
+ " \n",
+ " def __call__(self, n_realisations=1):\n",
+ " from PySDM.physics import si, constants_defaults\n",
+ "\n",
+ " assert n_realisations == 1, \"Only one realisation is supported\"\n",
+ "\n",
+ " pysdm_output = {}\n",
+ " for step in range(int(self.settings.T_MAX_SEC // self.settings.PLOT_TIME_STEP)+1):\n",
+ " if step != 0:\n",
+ " self.particulator.run(int(self.settings.PLOT_TIME_STEP // self.settings.DT_SEC))\n",
+ " if self.settings.BINNING_METHOD == 'Number Concentration (m^-3)':\n",
+ " pysdm_output[step] = self.particulator.products['spectrum'].get()\n",
+ " #next two lines convert from dr to dlnr\n",
+ " pysdm_output[step] *= (self.radius_bins_edges[1:]-self.radius_bins_edges[:-1])\n",
+ " pysdm_output[step] /= (np.log(self.radius_bins_edges[1:]) - np.log(self.radius_bins_edges[:-1]))\n",
+ " elif self.settings.BINNING_METHOD == 'Mass Density (g/m^3 dlnr)': \n",
+ " pysdm_output[step] = self.particulator.products['spectrum'].get()[0]\n",
+ " pysdm_output[step][:] *= constants_defaults.rho_w * si.kilogram / si.metre**3 \n",
+ " pysdm_output[step][:] /= si.g\n",
+ " return pysdm_output"
]
},
{
"cell_type": "code",
- "execution_count": 25,
- "id": "47b994f3",
+ "execution_count": 7,
+ "id": "80e1dcdd",
"metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "Overwriting script.jl\n"
- ]
- }
- ],
+ "outputs": [],
"source": [
- "%%writefile script.jl\n",
+ "class RunPartMC:\n",
+ " def __init__(self, settings):\n",
+ " import PyPartMC as ppmc\n",
+ " gas_data = ppmc.GasData((\"Background\",\"NO2\"))\n",
+ " env_state = ppmc.EnvState(\n",
+ " {\n",
+ " \"rel_humidity\": .99,\n",
+ " \"latitude\": 40,\n",
+ " \"longitude\": 0,\n",
+ " \"altitude\": 0 * ppmc.si.m,\n",
+ " \"start_time\": 0 * ppmc.si.s,\n",
+ " \"start_day\": 1,\n",
+ " }\n",
+ " )\n",
+ " env_state.additive_kernel_coefficient = settings.ADDITIVE_KERNEL_COEFF\n",
+ " assert env_state.additive_kernel_coefficient == settings.ADDITIVE_KERNEL_COEFF\n",
+ "\n",
+ " aero_data = ppmc.AeroData(\n",
+ " (\n",
+ " {\"H2O\": [1000 * ppmc.si.kg / ppmc.si.m**3, 0, 18.0 * ppmc.si.g / ppmc.si.mol, 0.00]},\n",
+ " )\n",
+ " )\n",
+ " gas_state = ppmc.GasState(gas_data)\n",
+ " placeholder_aero = {\n",
+ " \"mass_frac\": [{\"H2O\": [1]}],\n",
+ " \"diam_type\": \"geometric\",\n",
+ " \"mode_type\": \"log_normal\",\n",
+ " \"num_conc\": 0 / ppmc.si.m**3,\n",
+ " \"geom_mean_diam\": 0.02 * ppmc.si.um,\n",
+ " \"log10_geom_std_dev\": 0.161,\n",
+ " }\n",
+ " placeholder_gas = [\n",
+ " {\"time\": [0 * ppmc.si.s]},\n",
+ " {\"rate\": [0 / ppmc.si.s]},\n",
+ " {\"dist\": [[{\"placeholder\": placeholder_aero,}]]},\n",
+ " ]\n",
+ " \n",
+ " pressure_timeseries = [100000]\n",
+ " temp_timeseries = [300]\n",
+ " height_timeseries = [1000]\n",
+ "\n",
+ " scenario = ppmc.Scenario(\n",
+ " gas_data,\n",
+ " aero_data,\n",
+ " {\n",
+ " \"temp_profile\": [{\"time\": [0]}, {\"temp\": temp_timeseries}],\n",
+ " \"pressure_profile\": [\n",
+ " {\"time\": [0]},\n",
+ " {\"pressure\": pressure_timeseries},\n",
+ " ],\n",
+ " \"height_profile\": [{\"time\": [0]}, {\"height\": height_timeseries}],\n",
+ " \"gas_emissions\": placeholder_gas,\n",
+ " \"gas_background\": placeholder_gas,\n",
+ " \"aero_emissions\": placeholder_gas,\n",
+ " \"aero_background\":placeholder_gas,\n",
+ " \"loss_function\": \"none\",\n",
+ " },\n",
+ " )\n",
"\n",
- "using Pkg\n",
- "Pkg.add([\"Combinatorics\",\"Distributions\",\"Random\",\"JSON\",\"DelimitedFiles\",\"CPUTime\"])\n",
- "using Random\n",
- "using Combinatorics\n",
- "using Distributions\n",
- "using CPUTime\n",
- "include(\"Droplets.jl-src-SDfunc-coalescence.jl\")\n",
- "include(\"Droplets.jl-Examples-DSDvis.jl\")\n",
- "include(\"Droplets.jl-Examples-testfunctions.jl\")\n",
- "using DelimitedFiles\n",
- "using JSON\n",
+ " T_INITIAL = 0.0\n",
+ " scenario.init_env_state(env_state, T_INITIAL)\n",
"\n",
- " \n",
- "FT = Float64\n",
- "setup = JSON.parsefile(\"./setup.json\")\n",
+ " AERO_DIST_INIT = [\n",
+ " {\n",
+ " \"init1\": {\n",
+ " \"mass_frac\": [{\"H2O\": [1]}],\n",
+ " \"diam_type\": \"geometric\",\n",
+ " \"mode_type\": \"exp\",\n",
+ " \"num_conc\": settings.NUM_CONC_PER_M3 / ppmc.si.m**3,\n",
+ " \"diam_at_mean_vol\": settings.DIAM_AT_MEAN_VOL_M * ppmc.si.m,\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\": 'additive',\n",
+ " \"do_coagulation\": True,\n",
+ " \"coag_kernel\": \"additive\",\n",
+ " \"t_max\": settings.T_MAX_SEC * ppmc.si.s,\n",
+ " \"del_t\": settings.DT_SEC * ppmc.si.s,\n",
+ " }\n",
+ " )\n",
+ "\n",
+ " aero_state = ppmc.AeroState(aero_data, settings.N_PART, \"flat\")\n",
+ " aero_state.dist_sample(aero_dist_init)\n",
+ "\n",
+ " camp_core = ppmc.CampCore()\n",
+ " photolysis = ppmc.Photolysis()\n",
+ "\n",
+ " self.run_part_args = (\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",
+ " )\n",
+ "\n",
+ " self.aero_state = aero_state\n",
+ " self.gas_state = gas_state\n",
+ " self.rad_grid = ppmc.BinGrid(settings.N_BINS, \"log\", settings.R_BINS_MIN, settings.R_BINS_MAX)\n",
+ " self.settings = settings\n",
+ " \n",
+ " def __call__(self, n_realisations=1):\n",
+ " import PyPartMC as ppmc\n",
+ " assert n_realisations == 1, \"Only one realisation is supported\"\n",
+ "\n",
+ " N_STEPS = int(self.settings.T_MAX_SEC / self.settings.DT_SEC)\n",
+ "\n",
+ " dists = []\n",
+ " diameters = self.aero_state.diameters()\n",
+ " num_concs = self.aero_state.num_concs\n",
+ " masses = self.aero_state.masses()\n",
+ " mass_concs = np.array(num_concs) * np.array(masses)\n",
+ " dists_mass = []\n",
+ " dists.append(ppmc.histogram_1d(self.rad_grid, np.array(diameters)/2, num_concs))\n",
+ " dists_mass.append(ppmc.histogram_1d(self.rad_grid, np.array(diameters)/2, mass_concs / ppmc.si.g))\n",
"\n",
- "Ns::Int = setup[\"Ns\"]\n",
- "FT = Float64\n",
- "num_bins::Int = setup[\"num_bins\"]\n",
- "radius_bins_edges = 10 .^ range(log10(setup[\"radius_bins_log_low\"]), log10(setup[\"radius_bins_log_high\"]), length=setup[\"num_bin_edges\"]) \n",
- "R0 = FT(setup[\"d_mean\"])/2\n",
+ " last_output_time = 0.\n",
+ " last_progress_time = 0.\n",
+ " i_output = 1\n",
+ " t_initial = 0\n",
+ " for i_time in range(1,N_STEPS + 1):\n",
+ " (last_output_time, last_progress_time, i_output) = ppmc.run_part_timestep(\n",
+ " *self.run_part_args,\n",
+ " i_time,\n",
+ " t_initial,\n",
+ " last_output_time,\n",
+ " last_progress_time,\n",
+ " i_output\n",
+ " )\n",
+ " \n",
+ " if np.mod(i_time * self.settings.DT_SEC, self.settings.PLOT_TIME_STEP) == 0:\n",
+ " diameters = self.aero_state.diameters()\n",
+ " num_concs = self.aero_state.num_concs\n",
+ " masses = self.aero_state.masses()\n",
+ " mass_concs = np.array(num_concs) * np.array(masses)\n",
+ " dists.append(ppmc.histogram_1d(self.rad_grid, np.array(diameters)/2, num_concs))\n",
+ " dists_mass.append(ppmc.histogram_1d(self.rad_grid, np.array(diameters)/2, mass_concs / ppmc.si.g))\n",
"\n",
- "if setup[\"binning_method\"] == \"Number Concentration (m^-3)\"\n",
- " binning_method = number_density\n",
- "elseif setup[\"binning_method\"] == \"Mass Density (g/m^3 dlnr)\"\n",
- " binning_method = mass_density_lnr\n",
- "end\n",
+ " if self.settings.BINNING_METHOD == 'Number Concentration (m^-3)':\n",
+ " PyPartDists = dists\n",
+ " elif self.settings.BINNING_METHOD == 'Mass Density (g/m^3 dlnr)': \n",
+ " PyPartDists = dists_mass\n",
+ " return PyPartDists\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "f7d54f29",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "class RunDropletsJL:\n",
+ " def __init__(self, settings):\n",
+ " BASE_URL = (\n",
+ " # replace main with commit id for previous commit\n",
+ " 'https://raw.githubusercontent.com/emmacware/Droplets.jl/refs/heads/main/' \n",
+ " )\n",
+ " for path in (\n",
+ " 'Examples/testfunctions.jl', 'Examples/DSDvis.jl','src/SDfunc/coalescence.jl'\n",
+ " ):\n",
+ " with open('Droplets.jl-' + path.replace('/','-'), 'w', encoding='utf-8') as fout:\n",
+ " with urllib.request.urlopen(BASE_URL + path) as fin:\n",
+ " fout.write(fin.read().decode('utf-8'))\n",
+ " \n",
+ " code_to_write = \"\"\"\n",
+ " using Pkg\n",
+ " Pkg.add([\"Combinatorics\", \"Distributions\", \"Random\", \"JSON\", \"DelimitedFiles\", \"CPUTime\", \"Plots\"])\n",
+ " using Random,Combinatorics,Distributions,CPUTime\n",
+ " include(\"Droplets.jl-src-SDfunc-coalescence.jl\")\n",
+ " include(\"Droplets.jl-Examples-DSDvis.jl\")\n",
+ " include(\"Droplets.jl-Examples-testfunctions.jl\")\n",
+ " using DelimitedFiles,JSON,Plots\n",
+ " \n",
+ " FT = Float64\n",
+ " setup = JSON.parsefile(\"./setup.json\")\n",
"\n",
+ " Ns::Int = setup[\"N_PART\"]\n",
+ " FT = Float64\n",
+ " num_bins::Int = setup[\"N_BINS\"]\n",
+ " radius_bins_edges = 10 .^ range(log10(setup[\"R_BINS_MIN\"]), log10(setup[\"R_BINS_MAX\"]), length=setup[\"N_BIN_EDGES\"]) \n",
+ " R0 = FT(setup[\"DIAM_AT_MEAN_VOL_M\"])/2\n",
+ " T_MAX_SEC = setup[\"T_MAX_SEC\"]\n",
+ " PLOT_TIME_STEP = setup[\"PLOT_TIME_STEP\"]\n",
+ " output_steps = collect(0:PLOT_TIME_STEP:T_MAX_SEC)\n",
"\n",
- "coagsettings = coag_settings{FT}(Ns=Ns,Δt=setup[\"dt\"],ΔV=setup[\"V\"],\n",
- " golovin_kernel_coeff=FT(setup[\"b\"]), n0=FT(setup[\"n_0\"]),R0=R0)\n",
+ " if setup[\"BINNING_METHOD\"] == \"Number Concentration (m^-3)\"\n",
+ " binning_method = number_density\n",
+ " elseif setup[\"BINNING_METHOD\"] == \"Mass Density (g/m^3 dlnr)\"\n",
+ " binning_method = mass_density_lnr\n",
+ " end\n",
"\n",
- "runsettings=run_settings{FT}(num_bins=num_bins,radius_bins_edges=radius_bins_edges,\n",
- " smooth=setup[\"smooth\"],coag_threading=eval(Meta.parse(setup[\"coag_threading\"])),\n",
- " output_steps=setup[\"output_steps\"],init_method=eval(Meta.parse(setup[\"init_method\"])),binning_method = binning_method)\n",
"\n",
- "drops = runsettings.init_method(coagsettings)\n",
- "bin,timing = coag_runtime(1,drops,coagsettings,runsettings)\n",
+ " coagsettings = coag_settings{FT}(Ns=Ns,Δt=setup[\"DT_SEC\"],ΔV=setup[\"VOLUME_M3\"],\n",
+ " golovin_kernel_coeff=FT(setup[\"ADDITIVE_KERNEL_COEFF\"]), n0=FT(setup[\"NUM_CONC_PER_M3\"]),R0=R0)\n",
+ " runsettings=run_settings{FT}(num_bins=num_bins,radius_bins_edges=radius_bins_edges,\n",
+ " smooth = false,output_steps=output_steps,init_method=init_logarithmic,binning_method = binning_method)\n",
"\n",
+ " drops = runsettings.init_method(coagsettings)\n",
+ " bin,timing = coag_runtime(1,drops,coagsettings,runsettings)\n",
"\n",
- "json_string = JSON.json(bin)\n",
"\n",
- "open(\"output.json\", \"w\") do file\n",
- " write(file, json_string)\n",
- "end\n",
+ " json_string = JSON.json(bin)\n",
+ " open(\"output.json\", \"w\") do file\n",
+ " write(file, json_string)\n",
+ " end\n",
+ " \"\"\"\n",
"\n",
- "println(\"output written to output.json\")"
+ " with open('script.jl', 'w') as file:\n",
+ " file.write(code_to_write)\n",
+ " def __call__(self, n_realisations=1):\n",
+ " subprocess.run([\"julia\", \"script.jl\"], check=True)\n",
+ " with open('output.json', 'r', encoding='utf8') as file:\n",
+ " Droplets = json.load(file)\n",
+ " return Droplets"
]
},
{
"cell_type": "code",
- "execution_count": 26,
- "id": "21550883",
+ "execution_count": 9,
+ "id": "b33d3ccb",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
- "The latest version of Julia in the `release` channel is 1.11.2+0.aarch64.apple.darwin14. You currently have `1.10.4+0.aarch64.apple.darwin14` installed. Run:\n",
- "\n",
- " juliaup update\n",
- "\n",
- "in your terminal shell to install Julia 1.11.2+0.aarch64.apple.darwin14 and update the `release` channel to that version.\n",
- " Updating registry at `~/.julia/registries/General.toml`\n",
- " Resolving package versions...\n",
- " No Changes to `~/.julia/environments/v1.10/Project.toml`\n",
- " No Changes to `~/.julia/environments/v1.10/Manifest.toml`\n"
+ "/Users/emmaware/Library/Python/3.9/lib/python/site-packages/PySDM/backends/numba.py:46: UserWarning: Disabling Numba threading due to ARM64 CPU (atomics do not work yet)\n",
+ " warnings.warn(\n",
+ "/Users/emmaware/Library/Python/3.9/lib/python/site-packages/PySDM/backends/numba.py:46: UserWarning: Disabling Numba threading due to ARM64 CPU (atomics do not work yet)\n",
+ " warnings.warn(\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
- "Running simulation...\n",
- "simtime =2.814713\n",
- "coal_func_time =2.5472279999999996\n",
- "output written to output.json\n"
+ "Running analytical\n",
+ "Running partmc\n",
+ "Running pysdm\n",
+ "Running droplets\n"
]
},
{
- "data": {
- "text/plain": [
- "CompletedProcess(args=['julia', 'script.jl'], returncode=0)"
- ]
- },
- "execution_count": 26,
- "metadata": {},
- "output_type": "execute_result"
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ " Resolving package versions...\n",
+ " No Changes to `~/.julia/environments/v1.11/Project.toml`\n",
+ " No Changes to `~/.julia/environments/v1.11/Manifest.toml`\n"
+ ]
+ },
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Running simulation...\n",
+ "simtime =3.083128\n",
+ "coal_func_time =2.790316\n"
+ ]
}
],
"source": [
- "subprocess.run([\"julia\", \"script.jl\"], check=True)\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 27,
- "id": "3aae69e6",
- "metadata": {},
- "outputs": [],
- "source": [
- "with open('output.json', 'r', encoding='utf8') as file:\n",
- " Droplets = json.load(file)"
+ "models = {\n",
+ " 'analytical': AnalyticalSoln(settings),\n",
+ " 'partmc': RunPartMC(settings),\n",
+ " 'pysdm': RunPySDM(settings),\n",
+ " 'droplets':RunDropletsJL(settings),\n",
+ "}\n",
+ "output = {}\n",
+ "for key, model in models.items():\n",
+ " print(f\"Running {key}\")\n",
+ " output[key] = model()\n"
]
},
{
"cell_type": "code",
- "execution_count": 28,
- "id": "829dd055",
+ "execution_count": 10,
+ "id": "40666710",
"metadata": {},
"outputs": [
{
@@ -4260,12 +535,12 @@
"\n",
"\n",
- "