PySDM_examples.Arabas_et_al_2025.run_simulation 

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