Skip to main content
Version: Develop 🚧

Advanced Physical Modelling

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.

import uclchem
import matplotlib.pyplot as plt
import os

The Hot Core​

Initial Conditions (Phase 1)​

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.

To do this, we will use uclchem.model.cloud() to run a model where a cloud of gas collapses from a density of 102cm−310^2 cm^{-3} to our hot core density of 106cm−310^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.

# set a parameter dictionary for cloud collapse model
param_dict = {
    "endAtFinalDensity": False,  # stop at finalTime
    "freefall": True,  # increase density in freefall
    "initialDens": 1e2,  # starting density
    "finalDens": 1e6,  # final density
    "initialTemp": 10.0,  # temperature of gas
    "finalTime": 6.0e6,  # final time
    "rout": 0.1,  # radius of cloud in pc
    "baseAv": 1.0,  # visual extinction at cloud edge.
    "abundSaveFile": "../examples/test-output/startcollapse.dat",  # save final abundances to file
    "outputFile": "../examples/test-output/phase1-full.dat",
}
if not os.path.exists("../examples/test-output/"):
    os.makedirs("../examples/test-output/")
result = uclchem.model.cloud(param_dict=param_dict)
print(result)

[0]

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.

Running the Science Model (Phase 2)​

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.

We also want to turn off freefall and change how long the model runs for.

# change other bits of input to set up phase 2
param_dict["initialDens"] = 1e6
param_dict["finalTime"] = 1e6
param_dict["freefall"] = False

# freeze out is completely overwhelmed by thermal desorption
# so turning it off has no effect on abundances but speeds up integrator.
param_dict["freezeFactor"] = 0.0

param_dict["abstol_factor"] = 1e-18
param_dict["reltol"] = 1e-12

# pop is dangerous, it removes the original key so you can't rerun this cell.
param_dict["abundLoadFile"] = param_dict.pop("abundSaveFile")
param_dict["outputFile"] = "../examples/test-output/phase2-full.dat"

result = uclchem.model.hot_core(
    temp_indx=3, max_temperature=300.0, param_dict=param_dict
)

At T(=R1) and step size H(=R2), the
corrector convergence failed repeatedly
or with ABS(H) = HMIN.

In the above message, R1 = 0.7284566946643D+13 R2 = 0.7805477683448D+00 ISTATE -5 - shortening step at time 230524.27046337884 years

Note that we've changed made two changes to the parameters here which aren't strictly necessary but can be helpful in certain situations.

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.

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.

Checking the Result​

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.

phase2_df = uclchem.analysis.read_output_file("../examples/test-output/phase2-full.dat")
uclchem.analysis.check_element_conservation(phase2_df)

{'H': '0.002%', 'N': '0.000%', 'C': '0.000%', 'O': '0.000%'}

phase2_df
TimeDensitygasTempdustTempavradfieldzetapointH#H...HOSO+#HS2@HS2H2S2+H2S2#H2S2@H2S2E-BULKSURFACE
00.000000e+001000000.010.010.0193.881.01.012.786290e-085.448210e-12...2.062680e-261.029490e-213.223450e-088.196880e-208.349650e-187.682020e-144.626350e-081.292190e-084.601420e-014.981620e-06
11.000000e-071000000.010.010.0193.881.01.012.786290e-085.448210e-12...2.062680e-262.028010e-193.223450e-088.196880e-208.349650e-187.682050e-144.626350e-081.292190e-084.601410e-014.981620e-06
21.000000e-061000000.010.010.0193.881.01.012.786290e-085.448170e-12...2.062680e-262.018750e-183.223450e-088.196880e-208.349650e-187.682310e-144.626350e-081.292190e-084.601410e-014.981620e-06
31.000000e-051000000.010.010.0193.881.01.012.786290e-085.447770e-12...2.062680e-262.017810e-173.223450e-088.196880e-208.349650e-187.684920e-144.626350e-081.292190e-084.601410e-014.981620e-06
41.000000e-041000000.010.010.0193.881.01.012.786290e-085.443610e-12...2.062680e-262.017640e-163.223450e-088.196880e-208.349650e-187.710980e-144.626350e-081.292190e-084.601410e-014.981620e-06
..................................................................
1899.601000e+051000000.0300.0300.0193.881.01.012.624020e-062.543690e-30...6.347290e-161.000000e-301.000000e-305.819460e-142.228320e-081.000000e-301.000000e-302.302190e-088.300000e-297.936150e-23
1909.701000e+051000000.0300.0300.0193.881.01.012.628700e-062.565080e-30...6.164570e-161.000000e-301.000000e-305.607780e-142.209860e-081.000000e-301.000000e-302.292690e-088.300000e-297.936150e-23
1919.801000e+051000000.0300.0300.0193.881.01.012.633290e-062.586490e-30...5.982640e-161.000000e-301.000000e-305.402890e-142.191550e-081.000000e-301.000000e-302.283210e-088.300000e-297.936150e-23
1929.901000e+051000000.0300.0300.0193.881.01.012.637800e-062.607900e-30...5.801690e-161.000000e-301.000000e-305.204660e-142.173400e-081.000000e-301.000000e-302.273750e-088.300000e-297.936150e-23
1931.000000e+061000000.0300.0300.0193.881.01.012.642210e-062.629310e-30...5.621920e-161.000000e-301.000000e-305.012950e-142.155400e-081.000000e-301.000000e-302.264300e-088.300000e-297.936150e-23

194 rows × 343 columns

species = ["CO", "H2O", "CH3OH", "#CO", "#H2O", "#CH3OH", "@H2O", "@CO", "@CH3OH"]
fig, [ax, ax2] = plt.subplots(1, 2, figsize=(16, 9))
ax = uclchem.analysis.plot_species(ax, phase2_df, species)
settings = ax.set(
    yscale="log",
    xlim=(1e2, 1e6),
    ylim=(1e-10, 1e-2),
    xlabel="Time / years",
    ylabel="Fractional Abundance",
    xscale="log",
)

ax2.plot(phase2_df["Time"], phase2_df["Density"], color="black")
ax2.set(xscale="log")
ax3 = ax2.twinx()
ax3.plot(phase2_df["Time"], phase2_df["gasTemp"], color="red")
ax2.set(xlabel="Time / year", ylabel="Density")
ax3.set(ylabel="Temperature", facecolor="red", xlim=(1e2, 1e6))
ax3.tick_params(axis="y", colors="red")

CO H2O CH3OH #CO #H2O #CH3OH @H2O @CO @CH3OH

png

Here, we see the value of running a collapse phase before the science run. Having run a collapse, we start this model with well developed ices and having material in the surface and bulk allows us to properly model the effect of warm up in a hot core. For example, the @CO abundance is ∼10−4\sim10^{-4} and #CO is ∼10−6\sim10^{-6}. As the gas warms to around 30K, the #CO abundance drops drastically as CO's binding energy is such that it is efficiently desorbed from the surface at this temperature. However, the rest of the CO is trapped in the bulk, surrounded by more strongly bound H2O molecules. Thus, the @CO abundance stays high until the gas reaches around 130K, when the H2O molecules are released along with the entire bulk.

Shocks​

Essentially the same process should be followed for shocks. Let's run a C-type and J-type shock through a gas of density 104cm−310^4 cm^{-3}. Again, we first run a simple cloud model to obtain some reasonable starting abundances, then we can run the shocks.

# set a parameter dictionary for phase 1 collapse model

param_dict = {
    "endAtFinalDensity": False,  # stop at finalTime
    "freefall": True,  # increase density in freefall
    "initialDens": 1e2,  # starting density
    "finalDens": 1e4,  # final density
    "initialTemp": 10.0,  # temperature of gas
    "finalTime": 6.0e6,  # final time
    "rout": 0.1,  # radius of cloud in pc
    "baseAv": 1.0,  # visual extinction at cloud edge.
    "abundSaveFile": "../examples/test-output/shockstart.dat",
}
result = uclchem.model.cloud(param_dict=param_dict)

C-shock​

We'll first run a c-shock. We'll run a 40 km s −1^{-1} shock through a gas of density 10410^4 cm −3^{-3}, using the abundances we just produced. Note that c-shock is the only model which returns an additional output in its result list. Not only is the first element the success flag indicating whether UCLCHEM completed, the second element is the dissipation time of the shock. We'll use that time to make our plots look nicer, cutting to a reasonable time. You can also obtain it from uclchem.utils.cshock_dissipation_time().

# change other bits of input to set up phase 2
param_dict["initialDens"] = 1e4
param_dict["finalTime"] = 1e6
if "abundSaveFile" in param_dict:
    param_dict.pop("abundSaveFile")
param_dict["abundLoadFile"] = "../examples/test-output/shockstart.dat"
param_dict["outputFile"] = "../examples/test-output/cshock.dat"


result = uclchem.model.cshock(shock_vel=40, param_dict=param_dict)
result, dissipation_time = result

Cannot have freefall on during cshock setting freefall=0 and continuing

The code completes fine. We do get a couple of warnings though. First, we're informed that freefall must be set to False for the C-shock model. Then we get a few integrator warnings. These are not important and can be ignored as long as the element conservation looks ok. However, it is an indication that the integrator did struggle with these ODEs under these conditions.

phase2_df = uclchem.analysis.read_output_file("../examples/test-output/cshock.dat")
uclchem.analysis.check_element_conservation(phase2_df)

{'H': '0.005%', 'N': '1.257%', 'C': '1.458%', 'O': '1.302%'}

species = ["CO", "H2O", "CH3OH", "NH3", "$CO", "$H2O", "$CH3OH", "$NH3"]

fig, [ax, ax2] = plt.subplots(1, 2, figsize=(16, 9))
ax = uclchem.analysis.plot_species(ax, phase2_df, species)
settings = ax.set(
    yscale="log",
    xlim=(1, 20 * dissipation_time),
    ylim=(1e-10, 1e-2),
    xlabel="Time / years",
    ylabel="Fractional Abundance",
    xscale="log",
)

ax2.plot(phase2_df["Time"], phase2_df["Density"], color="black")
ax2.set(xscale="log")
ax3 = ax2.twinx()
ax3.plot(phase2_df["Time"], phase2_df["gasTemp"], color="red")
ax2.set(xlabel="Time / year", ylabel="Density")
ax3.set(ylabel="Temperature", facecolor="red", xlim=(1, 20 * dissipation_time))
ax3.tick_params(axis="y", colors="red")

CO H2O CH3OH NH3 COCO H2O CH3OHCH3OH NH3

png

J-shock​

Running a j-shock is a simple case of changing function. We'll run a 10 km s −1^{-1} shock through a gas of density 10310^3 cm −3^{-3} gas this time. Note that nothing stops us using the intial abundances we produced for the c-shock. UCLCHEM will not check that the initial density matches the density of the abundLoadFile. It may not always be a good idea to do this but we should remember the intial abundances really are just a rough approximation.

param_dict["initialDens"] = 1e3
param_dict["freefall"] = False  # lets remember to turn it off this time
param_dict["reltol"] = 1e-12

shock_vel = 10.0

if "abundSaveFile" in param_dict:
    param_dict.pop("abundSaveFile")
param_dict["abundLoadFile"] = "../examples/test-output/shockstart.dat"
param_dict["outputFile"] = "../examples/test-output/jshock.dat"

result = uclchem.model.jshock(shock_vel=shock_vel, param_dict=param_dict)
print(result)

[0]

This time, we've turned off the freefall option and made reltol a little more stringent. The j-shock ends up running a bit slower but we get no warnings on this run.

phase2_df = uclchem.analysis.read_output_file(param_dict["outputFile"])
uclchem.analysis.check_element_conservation(phase2_df)

{'H': '0.103%', 'N': '1.581%', 'C': '1.322%', 'O': '1.578%'}

species = ["CO", "H2O", "CH3OH", "NH3", "$CO", "$H2O", "$CH3OH", "$NH3"]

fig, [ax, ax2] = plt.subplots(1, 2, figsize=(16, 9))
ax = uclchem.analysis.plot_species(ax, phase2_df, species)
settings = ax.set(
    yscale="log",
    xlim=(1e-7, 1e6),
    ylim=(1e-10, 1e-2),
    xlabel="Time / years",
    ylabel="Fractional Abundance",
    xscale="log",
)

ax2.plot(phase2_df["Time"], phase2_df["Density"], color="black")
ax2.set(xscale="log", yscale="log")
ax3 = ax2.twinx()
ax3.plot(phase2_df["Time"], phase2_df["gasTemp"], color="red")
ax2.set(xlabel="Time / year", ylabel="Density")
ax3.set(ylabel="Temperature", facecolor="red", xlim=(1e-7, 1e6))
ax3.tick_params(axis="y", colors="red")

CO H2O CH3OH NH3 COCO H2O CH3OHCH3OH NH3

png

That's everything! We've run various science models using reasonable starting abundances that we produced by running a simple UCLCHEM model beforehand. One benefit of this method is that the abundances are consistent with the network. If we start with arbitrary, perhaps observationally motivated, abundances, it would be possible to initiate the model in a state our network could never produce.

However, one should be aware of the limitations of this method. A freefall collapse from low density to high is not really how a molecular cloud forms and so the abundances are only approximately similar to values they'd truly have in a real cloud. Testing whether your results are sensitive to things like the time you run the preliminary for or the exact density is a good way to make sure these approximations are not problematic.

Bear in mind that you can use abundSaveFile and abundLoadFile in the same model run. This lets you chain model runs together. For example, you could run a c-shock from a cloud model as we did here and then a j-shock with the c-shock's abundances as the initial abundances.

Last updated on by Gijs Vermariën