PySDM_examples.Murray_et_al_2010.simulation

  1import numpy as np
  2from matplotlib import pyplot
  3from matplotlib.ticker import MultipleLocator
  4
  5from PySDM import Particulator, Formulae
  6from PySDM.backends import CPU
  7from PySDM.dynamics import Freezing
  8from PySDM.environments import Box
  9from PySDM.physics import si
 10
 11
 12class Simulation:
 13    def __init__(
 14        self,
 15        *,
 16        cases,
 17        n_sd=1000,
 18        n_runs_per_case=1,
 19        temperature_step=0.01 * si.K,
 20        droplet_radius=10 * si.um,
 21        homogeneous_ice_nucleation_rate="KoopMurray2016",
 22        temperature_range=None,
 23    ):
 24        self.cases = cases
 25        self.n_sd = n_sd
 26        self.n_runs_per_case = n_runs_per_case
 27        self.temperature_step = temperature_step
 28        self.droplet_radius = droplet_radius
 29        self.homogeneous_ice_nucleation_rate = homogeneous_ice_nucleation_rate
 30        if temperature_range is None:
 31            temperature_range = [234.5 * si.K, 238 * si.K]
 32        self.temperature_range = temperature_range
 33        self.multiplicity = 1
 34        self.volume = 1 * si.cm**3
 35        self.output = None
 36
 37    def run(self, keys):
 38        self.output = {}
 39        for key in keys:
 40            case = self.cases[key]
 41            assert case["cooling_rate"] != 0
 42            total_time = (
 43                np.diff(np.asarray(self.temperature_range)).item()
 44                / case["cooling_rate"]
 45            )
 46            time_step = self.temperature_step / case["cooling_rate"]
 47
 48            self.output[key] = []
 49            for i in range(self.n_runs_per_case):
 50                T_frz = simulation(
 51                    seed=i,
 52                    n_sd=self.n_sd,
 53                    time_step=time_step,
 54                    volume=self.volume,
 55                    droplet_radius=self.droplet_radius,
 56                    cooling_rate=case["cooling_rate"],
 57                    multiplicity=self.multiplicity,
 58                    total_time=total_time,
 59                    homogeneous_ice_nucleation_rate=self.homogeneous_ice_nucleation_rate,
 60                    initial_temperature=self.temperature_range[1],
 61                )
 62                self.output[key].append({"T_frz": T_frz})
 63
 64    def plot_histogram(self, title=None):
 65
 66        pyplot.rc("font", size=12)
 67        T_bins = np.arange(self.temperature_range[0], self.temperature_range[1], 0.1)
 68
 69        for key in self.output:
 70            for run in range(self.n_runs_per_case):
 71                if run == 0:
 72                    cooling_rate = self.cases[key]["cooling_rate"]
 73                    label = (
 74                        "-"
 75                        + f"{cooling_rate*si.minute:.1f}"
 76                        + r" $\mathrm{K \, min^{-1}}$"
 77                    )
 78                else:
 79                    label = None
 80
 81                pyplot.hist(
 82                    self.output[key][run]["T_frz"],
 83                    bins=T_bins,
 84                    cumulative=-1,
 85                    density=True,
 86                    alpha=1.0,
 87                    histtype="step",
 88                    linewidth=1.5,
 89                    color=self.cases[key]["color"],
 90                    label=label,
 91                )
 92
 93        if title is None:
 94            title = (
 95                r"$r_\mathrm{drop}:$ "
 96                + f"{self.droplet_radius/si.micrometer:.2f}"
 97                + "µm"
 98                + r"   $\mathrm{N_{sd}}:$ "
 99                + str(self.n_sd)
100                + "   dT: "
101                + f"{self.temperature_step:.2f}"
102                + "K"
103            )
104        pyplot.title(title, pad=15)
105        pyplot.xlim(*self.temperature_range)
106        pyplot.ylim(-0.05, 1.05)
107        pyplot.xlabel("Temperature / K")
108        pyplot.ylabel("Fraction of droplets frozen")
109        pyplot.minorticks_on()
110        pyplot.tick_params(axis="both", which="both", top=True, right=True)
111        pyplot.gca().yaxis.set_minor_locator(MultipleLocator(0.1))
112        pyplot.tick_params(which="minor", length=3)
113        pyplot.tick_params(which="major", length=6)
114        pyplot.legend()
115
116
117def simulation(
118    *,
119    seed,
120    n_sd,
121    time_step,
122    volume,
123    droplet_radius,
124    cooling_rate,
125    multiplicity,
126    total_time,
127    homogeneous_ice_nucleation_rate="KoopMurray2016",
128    initial_temperature=None,
129):
130    formulae = Formulae(
131        seed=seed,
132        homogeneous_ice_nucleation_rate=homogeneous_ice_nucleation_rate,
133        particle_shape_and_density="MixedPhaseSpheres",
134        saturation_vapour_pressure="MurphyKoop2005",
135    )
136    environment = Box(dt=time_step, dv=volume, backend=CPU(formulae=formulae))
137
138    droplet_volume = formulae.trivia.volume(radius=droplet_radius)
139
140    attributes = {
141        "multiplicity": np.full(n_sd, multiplicity),
142        "signed water mass": np.full(n_sd, droplet_volume * formulae.constants.rho_w),
143    }
144    particulator = Particulator(
145        n_sd=n_sd,
146        dynamics=(
147            Freezing(homogeneous_freezing="time-dependent", immersion_freezing=None),
148        ),
149        environment=environment,
150        attributes=attributes,
151        requested_attributes=("temperature of last freezing", "volume"),
152    )
153    env = particulator.environment
154
155    env["T"] = initial_temperature
156    env["a_w_ice"] = np.nan
157    env["RH"] = 1 + np.finfo(float).eps
158    svp = particulator.formulae.saturation_vapour_pressure
159
160    cell_id = 0
161    for i in range(int(total_time / time_step) + 1):
162        env["RH_ice"] = svp.pvs_water(env["T"][cell_id]) / svp.pvs_ice(
163            env["T"][cell_id]
164        )
165        env["a_w_ice"] = svp.pvs_ice(env["T"][cell_id]) / svp.pvs_water(
166            env["T"][cell_id]
167        )
168        env["T"] -= np.full((1,), cooling_rate * time_step)
169        particulator.run(0 if i == 0 else 1)
170
171    assert all(particulator.attributes["signed water mass"].data < 0)
172    T_frz = particulator.attributes["temperature of last freezing"].data.tolist()
173    return T_frz