PySDM_examples.Arabas_et_al_2023.run_simulation 

 1import numpy as np
 2
 3
 4def update_thermo(particulator, T):
 5    env = particulator.environment
 6    svp = particulator.formulae.saturation_vapour_pressure
 7
 8    env["T"] = T
 9    env["a_w_ice"] = svp.pvs_ice(T) / svp.pvs_water(T)
10
11
12def run_simulation(particulator, temperature_profile, n_steps):
13    output = {
14        "products": {k: [] for k in particulator.products.keys()},
15        "attributes": [],
16    }
17    for step in range(n_steps + 1):
18        if step != 0:
19            update_thermo(
20                particulator, T=temperature_profile((step - 0.5) * particulator.dt)
21            )
22            particulator.run(step - particulator.n_steps)
23            update_thermo(particulator, T=temperature_profile(step * particulator.dt))
24            output["frozen"].append(particulator.attributes["volume"].to_ndarray() < 0)
25        else:
26            output["spectrum"] = {}
27            output["frozen"] = [np.full(particulator.n_sd, False)]
28            for k in ("multiplicity", "freezing temperature", "immersed surface area"):
29                if k in particulator.attributes:
30                    output["spectrum"][k] = particulator.attributes[k].to_ndarray()
31        for k, v in particulator.products.items():
32            output["products"][k].append(v.get() + 0)
33    return output
def update_thermo(particulator, T): 
 5def update_thermo(particulator, T):
 6    env = particulator.environment
 7    svp = particulator.formulae.saturation_vapour_pressure
 8
 9    env["T"] = T
10    env["a_w_ice"] = svp.pvs_ice(T) / svp.pvs_water(T)
def run_simulation(particulator, temperature_profile, n_steps): 
13def run_simulation(particulator, temperature_profile, n_steps):
14    output = {
15        "products": {k: [] for k in particulator.products.keys()},
16        "attributes": [],
17    }
18    for step in range(n_steps + 1):
19        if step != 0:
20            update_thermo(
21                particulator, T=temperature_profile((step - 0.5) * particulator.dt)
22            )
23            particulator.run(step - particulator.n_steps)
24            update_thermo(particulator, T=temperature_profile(step * particulator.dt))
25            output["frozen"].append(particulator.attributes["volume"].to_ndarray() < 0)
26        else:
27            output["spectrum"] = {}
28            output["frozen"] = [np.full(particulator.n_sd, False)]
29            for k in ("multiplicity", "freezing temperature", "immersed surface area"):
30                if k in particulator.attributes:
31                    output["spectrum"][k] = particulator.attributes[k].to_ndarray()
32        for k, v in particulator.products.items():
33            output["products"][k].append(v.get() + 0)
34    return output