PySDM_examples.Luettmer_homogeneous_freezing.simulation

  1import numpy as np
  2import PySDM.products as PySDM_products
  3from PySDM.builder import Builder
  4from PySDM.dynamics import (
  5    AmbientThermodynamics,
  6    Condensation,
  7    Freezing,
  8    VapourDepositionOnIce,
  9)
 10from PySDM.environments import Parcel
 11from PySDM.initialisation import discretise_multiplicities
 12from PySDM.initialisation.hygroscopic_equilibrium import equilibrate_dry_radii
 13
 14
 15class Simulation:
 16    def __init__(self, settings):
 17
 18        self.dt = settings.dt
 19
 20        formulae = settings.formulae
 21
 22        self.silent = settings.silent
 23
 24        env = Parcel(
 25            mixed_phase=True,
 26            dt=self.dt,
 27            mass_of_dry_air=settings.mass_of_dry_air,
 28            p0=settings.initial_pressure,
 29            initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio,
 30            T0=settings.initial_temperature,
 31            w=settings.w_updraft,
 32        )
 33
 34        builder = Builder(
 35            backend=settings.backend,
 36            n_sd=settings.n_sd,
 37            environment=env,
 38            dynamics=(
 39                [
 40                    AmbientThermodynamics(),
 41                    Condensation(adaptive=True),
 42                ]
 43                + (
 44                    [VapourDepositionOnIce(adaptive=True)]
 45                    if settings.deposition_enable
 46                    else []
 47                )
 48                + [
 49                    Freezing(
 50                        homogeneous_freezing=settings.hom_freezing_type,
 51                        immersion_freezing=None,
 52                    )
 53                ]
 54            ),
 55        )
 56
 57        self.n_sd = settings.n_sd
 58        self.multiplicities = discretise_multiplicities(
 59            settings.specific_concentration * env.mass_of_dry_air
 60        )
 61        self.r_wet = settings.r_wet
 62
 63        kappa = np.full_like(settings.r_wet, settings.kappa)
 64
 65        self.r_dry = equilibrate_dry_radii(
 66            r_wet=self.r_wet,
 67            environment=builder.particulator.environment,
 68            kappa=kappa,
 69        )
 70        v_dry = settings.formulae.trivia.volume(radius=self.r_dry)
 71        self.initial_mass = formulae.particle_shape_and_density.radius_to_mass(
 72            self.r_wet
 73        )
 74
 75        attributes = {
 76            "multiplicity": self.multiplicities,
 77            "dry volume": v_dry,
 78            "kappa times dry volume": kappa * v_dry,
 79            "signed water mass": self.initial_mass,
 80        }
 81        builder.request_attribute("temperature of last freezing")
 82        builder.request_attribute("supersaturation of last freezing")
 83        builder.request_attribute("radius")
 84        builder.request_attribute("wet to critical volume ratio")
 85
 86        products = [
 87            PySDM_products.ParcelDisplacement(name="z"),
 88            PySDM_products.Time(name="t"),
 89            PySDM_products.AmbientRelativeHumidity(name="RH", unit="%"),
 90            PySDM_products.AmbientRelativeHumidity(name="RH_ice", unit="%"),
 91            PySDM_products.AmbientTemperature(name="T"),
 92            PySDM_products.AmbientPressure(name="p", unit="hPa"),
 93            PySDM_products.WaterMixingRatio(name="LWC", radius_range=(0, np.inf)),
 94            PySDM_products.WaterMixingRatio(name="IWC", radius_range=(-np.inf, 0)),
 95            PySDM_products.AmbientWaterVapourMixingRatio(
 96                name="qv", var="water_vapour_mixing_ratio"
 97            ),
 98            PySDM_products.ParticleSpecificConcentration(
 99                name="ns", radius_range=(0, np.inf), unit="kg^-1"
100            ),
101            PySDM_products.ParticleSpecificConcentration(
102                name="ni", radius_range=(-np.inf, 0), unit="kg^-1"
103            ),
104            PySDM_products.MeanRadius(name="rs", radius_range=(0, np.inf)),
105            PySDM_products.MeanRadius(name="ri", radius_range=(-np.inf, 0)),
106        ]
107        self.output_product_list = [
108            "z",
109            "t",
110            "T",
111            "p",
112            "RH",
113            "RH_ice",
114            "LWC",
115            "IWC",
116            "qv",
117            "ns",
118            "ni",
119            "rs",
120            "ri",
121        ]
122
123        self.particulator = builder.build(attributes, products)
124
125        self.n_output = settings.n_output
126        if settings.n_output == 1:
127            self.n_substeps = 1
128        else:
129            self.n_substeps = int(self.n_output / self.dt)
130        self.t_max_duration = settings.t_max_duration
131
132    def save(self, output):
133        cell_id = 0
134
135        for key in self.output_product_list:
136            if key == "t":
137                output[key].append(self.particulator.products[key].get())
138            else:
139                output[key].append(self.particulator.products[key].get()[cell_id])
140
141        output["T_frz"] = self.particulator.attributes[
142            "temperature of last freezing"
143        ].data.tolist()
144        output["RHi_frz"] = self.particulator.attributes[
145            "supersaturation of last freezing"
146        ].data.tolist()
147        if not output["radius"]:
148            output["radius"] = self.particulator.attributes["radius"].data.tolist()
149            output["multiplicity"] = self.particulator.attributes[
150                "multiplicity"
151            ].data.tolist()
152
153    def run(self):
154
155        if not self.silent:
156            print("Starting simulation...")
157
158        output = {
159            "T_frz": [],
160            "RHi_frz": [],
161            "radius": [],
162            "multiplicity": [],
163        }
164        for key in self.output_product_list:
165            output[key] = []
166
167        self.save(output)
168
169        while True:
170
171            self.particulator.run(self.n_substeps)
172            self.save(output)
173
174            w_cr_v_ratio = self.particulator.attributes[
175                "wet to critical volume ratio"
176            ].data
177            sig_mass = self.particulator.attributes["signed water mass"].data
178            frozen = sig_mass < 0
179            unactivated = w_cr_v_ratio < 1
180            if any(frozen) and all(np.logical_or(frozen, unactivated)):
181                if not self.silent:
182                    print("all particles frozen or evaporated")
183                # Assert for water saturation
184                test_water_saturation = np.asarray(output["RH"])
185                # Sort out times before CCN activation & after first occurence of ice
186                test_water_saturation = np.where(
187                    np.asarray(output["rs"]) < 1e-6, 100.0, test_water_saturation
188                )
189                test_water_saturation = np.where(
190                    np.asarray(output["IWC"]) > 0.0, 100.0, test_water_saturation
191                )
192                if np.allclose(test_water_saturation, 100.0, rtol=5e-2) is False:
193                    print(
194                        "Warning: water saturation is too high outside "
195                        "of activation and mixed-phase environment"
196                    )
197
198                break
199            if output["t"][-1] >= self.t_max_duration * self.dt:
200                print("time exceeded")
201                break
202
203        return output