{ "cells": [ { "cell_type": "markdown", "id": "8ee51bf1", "metadata": {}, "source": [ "# Advanced Physical Modelling\n", "\n", "In the previous tutorial, we simply modelled the chemistry of a static cloud for 1 Myr. This is unlikely to meet everybody's modelling needs and UCLCHEM is capable of modelling much more complex environments such as hot cores and shocks. In this tutorial, we model both a hot core and a shock to explore how these models work and to demonstrate the workflow that the UCLCHEM team normally follow." ] }, { "cell_type": "code", "execution_count": 1, "id": "cdd6fc4e", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T14:04:56.242049Z", "iopub.status.busy": "2026-01-23T14:04:56.241790Z", "iopub.status.idle": "2026-01-23T14:04:57.505137Z", "shell.execute_reply": "2026-01-23T14:04:57.504302Z" } }, "outputs": [], "source": [ "import uclchem\n", "import matplotlib.pyplot as plt\n", "import os" ] }, { "cell_type": "markdown", "id": "f16b89d9", "metadata": {}, "source": [ "## The Hot Core\n", "\n", "### Initial Conditions (Phase 1)\n", "UCLCHEM typically starts with the gas in atomic/ionic form with no molecules. However, this clearly is not appropriate when modelling an object such as a hot core. In these objects, the gas is already evolved and there should be molecules in the gas phase as well as ice mantles on the dust. To allow for this, one must provide some initial abundances to the model. There are many ways to do this but we typically chose to run a preliminary model to produce our abundances. In many UCLCHEM papers, we refer to the preliminary model as *phase 1* and the science model as *phase 2*. Phase 1 simply models a collapsing cloud and phase 2 models the object in question.\n", "\n", "To do this, we will use `uclchem.model.cloud()` to run a model where a cloud of gas collapses from a density of $10^2 cm^{-3}$ to our hot core density of $10^6 cm^{-3}$, keeping all other parameters constant. During this collapse, chemistry will occur and we can assume the final abundances of this model will be reasonable starting abundances for the hot core." ] }, { "cell_type": "code", "execution_count": 2, "id": "eb791a2d", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T14:04:57.507410Z", "iopub.status.busy": "2026-01-23T14:04:57.507105Z", "iopub.status.idle": "2026-01-23T14:05:06.344419Z", "shell.execute_reply": "2026-01-23T14:05:06.343338Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0]\n" ] } ], "source": [ "# set a parameter dictionary for cloud collapse model\n", "param_dict = {\n", " \"endAtFinalDensity\": False, # stop at finalTime\n", " \"freefall\": True, # increase density in freefall\n", " \"initialDens\": 1e2, # starting density\n", " \"finalDens\": 1e6, # final density\n", " \"initialTemp\": 10.0, # temperature of gas\n", " \"finalTime\": 6.0e6, # final time\n", " \"rout\": 0.1, # radius of cloud in pc\n", " \"baseAv\": 1.0, # visual extinction at cloud edge.\n", " \"abundSaveFile\": \"../examples/test-output/startcollapse.dat\", # save final abundances to file\n", " \"outputFile\": \"../examples/test-output/phase1-full.dat\",\n", "}\n", "if not os.path.exists(\"../examples/test-output/\"):\n", " os.makedirs(\"../examples/test-output/\")\n", "result = uclchem.model.cloud(param_dict=param_dict)\n", "print(result)" ] }, { "cell_type": "markdown", "id": "f01c2a85", "metadata": {}, "source": [ "With that done, we now have a file containing the final abundances of a cloud of gas after this collapse: `param_dict[\"abundSaveFile\"]` we can pass this to our hot core model to use those abundances as our initial abundances.\n", "\n", "### Running the Science Model (Phase 2)\n", "\n", "We need to change just a few things in `param_dict` to set up the hot core model. The key one is that UCLCHEM saves final abundances to `abundSaveFile` but loads them from `abundLoadFile` so we need to swap that key over to make the abundances we just produced our initial abundances.\n", "\n", "We also want to turn off freefall and change how long the model runs for.\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "922ee4f3", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T14:05:06.346131Z", "iopub.status.busy": "2026-01-23T14:05:06.345954Z", "iopub.status.idle": "2026-01-23T14:06:39.605448Z", "shell.execute_reply": "2026-01-23T14:06:39.604538Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " At T(=R1) and step size H(=R2), the \n", " corrector convergence failed repeatedly \n", " or with ABS(H) = HMIN. \n", "In the above message, R1 = 0.7282869807308D+13 R2 = 0.1314147544832D+02\n", " ISTATE -5 - shortening step at time 230010.00000000000 years\n" ] } ], "source": [ "# change other bits of input to set up phase 2\n", "param_dict[\"initialDens\"] = 1e6\n", "param_dict[\"finalTime\"] = 1e6\n", "param_dict[\"freefall\"] = False\n", "\n", "# freeze out is completely overwhelmed by thermal desorption\n", "# so turning it off has no effect on abundances but speeds up integrator.\n", "param_dict[\"freezeFactor\"] = 0.0\n", "\n", "param_dict[\"abstol_factor\"] = 1e-18\n", "param_dict[\"reltol\"] = 1e-12\n", "\n", "# pop is dangerous, it removes the original key so you can't rerun this cell.\n", "param_dict[\"abundLoadFile\"] = param_dict.pop(\"abundSaveFile\")\n", "param_dict[\"outputFile\"] = \"../examples/test-output/phase2-full.dat\"\n", "\n", "result = uclchem.model.hot_core(\n", " temp_indx=3, max_temperature=300.0, param_dict=param_dict\n", ")" ] }, { "cell_type": "markdown", "id": "a4a5756d", "metadata": {}, "source": [ "Note that we've changed made two changes to the parameters here which aren't strictly necessary but can be helpful in certain situations.\n", "\n", "Since the gas temperature increases throughout a hot core model, freeze out is much slower than thermal desorption for all but the first few time steps. Turning it off doesn't affect the abundances but will speed up the solution.\n", "\n", "We also change abstol and reltol here, largely to demonstrate their use. They control the integrator accuracy and whilst making them smaller does slow down successful runs, it can make runs complete that stall completely otherwise or give correct solutions where lower tolerances allow issues like element conservation failure to sneak in. If your code does not complete or element conservation fails, you can change them.\n", "\n", "### Checking the Result\n", "With a successful run, we can check the output. We first load the file and check the abundance conservation, then we can plot it up." ] }, { "cell_type": "code", "execution_count": 4, "id": "78b676dc", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T14:06:39.607182Z", "iopub.status.busy": "2026-01-23T14:06:39.607009Z", "iopub.status.idle": "2026-01-23T14:06:39.678235Z", "shell.execute_reply": "2026-01-23T14:06:39.677408Z" } }, "outputs": [ { "data": { "text/plain": [ "{'H': '0.000%', 'N': '0.000%', 'C': '0.000%', 'O': '0.000%'}" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "phase2_df = uclchem.analysis.read_output_file(\"../examples/test-output/phase2-full.dat\")\n", "uclchem.analysis.check_element_conservation(phase2_df)" ] }, { "cell_type": "code", "execution_count": 5, "id": "d38613ed", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T14:06:39.680028Z", "iopub.status.busy": "2026-01-23T14:06:39.679856Z", "iopub.status.idle": "2026-01-23T14:06:39.705145Z", "shell.execute_reply": "2026-01-23T14:06:39.704410Z" } }, "outputs": [ { "data": { "text/html": [ "
| \n", " | Time | \n", "Density | \n", "gasTemp | \n", "dustTemp | \n", "Av | \n", "radfield | \n", "zeta | \n", "point | \n", "H | \n", "H+ | \n", "... | \n", "@OCS | \n", "@C4N | \n", "@SIC3 | \n", "@SO2 | \n", "@S2 | \n", "@HS2 | \n", "@H2S2 | \n", "E- | \n", "BULK | \n", "SURFACE | \n", "
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | \n", "0.000000e+00 | \n", "1000000.0 | \n", "10.0 | \n", "10.0 | \n", "193.88 | \n", "1.0 | \n", "1.0 | \n", "1 | \n", "2.772050e-08 | \n", "1.276110e-08 | \n", "... | \n", "7.695200e-07 | \n", "1.030070e-08 | \n", "1.158650e-30 | \n", "1.098210e-12 | \n", "1.481850e-08 | \n", "1.379990e-08 | \n", "2.825310e-08 | \n", "1.303430e-08 | \n", "4.601750e-01 | \n", "4.981620e-06 | \n", "
| 1 | \n", "1.000000e-07 | \n", "1000000.0 | \n", "10.0 | \n", "10.0 | \n", "193.88 | \n", "1.0 | \n", "1.0 | \n", "1 | \n", "2.772050e-08 | \n", "1.276110e-08 | \n", "... | \n", "7.695200e-07 | \n", "1.030070e-08 | \n", "1.158650e-30 | \n", "1.098210e-12 | \n", "1.481850e-08 | \n", "1.379990e-08 | \n", "2.825310e-08 | \n", "1.303430e-08 | \n", "4.601750e-01 | \n", "4.981610e-06 | \n", "
| 2 | \n", "1.000000e-06 | \n", "1000000.0 | \n", "10.0 | \n", "10.0 | \n", "193.88 | \n", "1.0 | \n", "1.0 | \n", "1 | \n", "2.772050e-08 | \n", "1.276110e-08 | \n", "... | \n", "7.695200e-07 | \n", "1.030070e-08 | \n", "1.158650e-30 | \n", "1.098210e-12 | \n", "1.481850e-08 | \n", "1.379990e-08 | \n", "2.825310e-08 | \n", "1.303430e-08 | \n", "4.601750e-01 | \n", "4.981610e-06 | \n", "
| 3 | \n", "1.000000e-05 | \n", "1000000.0 | \n", "10.0 | \n", "10.0 | \n", "193.88 | \n", "1.0 | \n", "1.0 | \n", "1 | \n", "2.772050e-08 | \n", "1.276110e-08 | \n", "... | \n", "7.695200e-07 | \n", "1.030070e-08 | \n", "1.158650e-30 | \n", "1.098210e-12 | \n", "1.481850e-08 | \n", "1.379990e-08 | \n", "2.825310e-08 | \n", "1.303430e-08 | \n", "4.601750e-01 | \n", "4.981610e-06 | \n", "
| 4 | \n", "1.000000e-04 | \n", "1000000.0 | \n", "10.0 | \n", "10.0 | \n", "193.88 | \n", "1.0 | \n", "1.0 | \n", "1 | \n", "2.772050e-08 | \n", "1.276110e-08 | \n", "... | \n", "7.695200e-07 | \n", "1.030070e-08 | \n", "1.158650e-30 | \n", "1.098210e-12 | \n", "1.481850e-08 | \n", "1.379990e-08 | \n", "2.825310e-08 | \n", "1.303430e-08 | \n", "4.601750e-01 | \n", "4.981610e-06 | \n", "
| ... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "
| 278 | \n", "9.600000e+05 | \n", "1000000.0 | \n", "300.0 | \n", "300.0 | \n", "193.88 | \n", "1.0 | \n", "1.0 | \n", "1 | \n", "2.603290e-06 | \n", "1.992350e-13 | \n", "... | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "3.855020e-08 | \n", "8.944510e-29 | \n", "6.178940e-22 | \n", "
| 279 | \n", "9.700000e+05 | \n", "1000000.0 | \n", "300.0 | \n", "300.0 | \n", "193.88 | \n", "1.0 | \n", "1.0 | \n", "1 | \n", "2.606930e-06 | \n", "1.988970e-13 | \n", "... | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "3.839700e-08 | \n", "8.953440e-29 | \n", "6.178940e-22 | \n", "
| 280 | \n", "9.800000e+05 | \n", "1000000.0 | \n", "300.0 | \n", "300.0 | \n", "193.88 | \n", "1.0 | \n", "1.0 | \n", "1 | \n", "2.610520e-06 | \n", "1.985600e-13 | \n", "... | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "3.824470e-08 | \n", "8.962360e-29 | \n", "6.178940e-22 | \n", "
| 281 | \n", "9.900000e+05 | \n", "1000000.0 | \n", "300.0 | \n", "300.0 | \n", "193.88 | \n", "1.0 | \n", "1.0 | \n", "1 | \n", "2.614070e-06 | \n", "1.982230e-13 | \n", "... | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "3.809350e-08 | \n", "8.971280e-29 | \n", "6.178940e-22 | \n", "
| 282 | \n", "1.000000e+06 | \n", "1000000.0 | \n", "300.0 | \n", "300.0 | \n", "193.88 | \n", "1.0 | \n", "1.0 | \n", "1 | \n", "2.617550e-06 | \n", "1.978880e-13 | \n", "... | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "1.000000e-30 | \n", "3.794320e-08 | \n", "8.980210e-29 | \n", "6.178940e-22 | \n", "
283 rows × 343 columns
\n", "