PySDM_examples.Murray_et_al_2010.simulation

  1import numpy as np
  2from matplotlib import pyplot
  3from matplotlib.ticker import MultipleLocator
  4
  5from PySDM import Builder, 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)) / case["cooling_rate"]
 44            )
 45            time_step = self.temperature_step / case["cooling_rate"]
 46
 47            self.output[key] = []
 48            for i in range(self.n_runs_per_case):
 49                T_frz = simulation(
 50                    seed=i,
 51                    n_sd=self.n_sd,
 52                    time_step=time_step,
 53                    volume=self.volume,
 54                    droplet_radius=self.droplet_radius,
 55                    cooling_rate=case["cooling_rate"],
 56                    multiplicity=self.multiplicity,
 57                    total_time=total_time,
 58                    homogeneous_ice_nucleation_rate=self.homogeneous_ice_nucleation_rate,
 59                    initial_temperature=self.temperature_range[1],
 60                )
 61                self.output[key].append({"T_frz": T_frz})
 62
 63    def plot_histogram(self, title=None):
 64
 65        pyplot.rc("font", size=12)
 66        T_bins = np.arange(self.temperature_range[0], self.temperature_range[1], 0.1)
 67
 68        for key in self.output:
 69            for run in range(self.n_runs_per_case):
 70                if run == 0:
 71                    cooling_rate = self.cases[key]["cooling_rate"]
 72                    label = (
 73                        "-"
 74                        + f"{cooling_rate*si.minute:.1f}"
 75                        + r" $\mathrm{K \, min^{-1}}$"
 76                    )
 77                else:
 78                    label = None
 79
 80                pyplot.hist(
 81                    self.output[key][run]["T_frz"],
 82                    bins=T_bins,
 83                    cumulative=-1,
 84                    density=True,
 85                    alpha=1.0,
 86                    histtype="step",
 87                    linewidth=1.5,
 88                    color=self.cases[key]["color"],
 89                    label=label,
 90                )
 91
 92        if title is None:
 93            title = (
 94                r"$r_\mathrm{drop}:$ "
 95                + f"{self.droplet_radius/si.micrometer:.2f}"
 96                + "µm"
 97                + r"   $\mathrm{N_{sd}}:$ "
 98                + str(self.n_sd)
 99                + "   dT: "
100                + f"{self.temperature_step:.2f}"
101                + "K"
102            )
103        pyplot.title(title, pad=15)
104        pyplot.xlim(*self.temperature_range)
105        pyplot.ylim(-0.05, 1.05)
106        pyplot.xlabel("Temperature / K")
107        pyplot.ylabel("Fraction of droplets frozen")
108        pyplot.minorticks_on()
109        pyplot.tick_params(axis="both", which="both", top=True, right=True)
110        pyplot.gca().yaxis.set_minor_locator(MultipleLocator(0.1))
111        pyplot.tick_params(which="minor", length=3)
112        pyplot.tick_params(which="major", length=6)
113        pyplot.legend()
114
115
116def simulation(
117    *,
118    seed,
119    n_sd,
120    time_step,
121    volume,
122    droplet_radius,
123    cooling_rate,
124    multiplicity,
125    total_time,
126    homogeneous_ice_nucleation_rate="KoopMurray2016",
127    initial_temperature=None,
128):
129    formulae = Formulae(
130        seed=seed,
131        homogeneous_ice_nucleation_rate=homogeneous_ice_nucleation_rate,
132        particle_shape_and_density="MixedPhaseSpheres",
133        saturation_vapour_pressure="MurphyKoop2005",
134    )
135    builder = Builder(
136        n_sd=n_sd,
137        backend=CPU(formulae=formulae),
138        environment=Box(dt=time_step, dv=volume),
139        dynamics=(
140            Freezing(homogeneous_freezing="time-dependent", immersion_freezing=None),
141        ),
142    )
143
144    builder.request_attribute("temperature of last freezing")
145    builder.request_attribute("volume")
146    droplet_volume = formulae.trivia.volume(radius=droplet_radius)
147
148    attributes = {
149        "multiplicity": np.full(n_sd, multiplicity),
150        "signed water mass": np.full(n_sd, droplet_volume * formulae.constants.rho_w),
151    }
152    particulator = builder.build(attributes=attributes)
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