PySDM_examples.featured.Abdul_Razzak_Ghan_2000.run_ARG_parcel

  1from collections import namedtuple
  2
  3import numpy as np
  4from PySDM_examples.featured.Abdul_Razzak_Ghan_2000.aerosol import (
  5    CONSTANTS_ARG,
  6    AerosolARG,
  7)
  8
  9from PySDM import Particulator, Formulae
 10from PySDM import products as PySDM_products
 11from PySDM.backends import CPU
 12from PySDM.dynamics import AmbientThermodynamics, Condensation
 13from PySDM.environments import Parcel
 14from PySDM.initialisation.hygroscopic_equilibrium import equilibrate_wet_radii
 15from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity
 16from PySDM.physics import si
 17
 18
 19def run_parcel(
 20    w,
 21    sol2,
 22    N2,
 23    rad2,
 24    n_sd_per_mode,
 25    *,
 26    RH0=1.0,
 27    T0=294 * si.K,
 28    p0=1e5 * si.Pa,
 29    n_steps=50,
 30    mass_of_dry_air=1e3 * si.kg,
 31    dt=2 * si.s,
 32):
 33    products = (
 34        PySDM_products.WaterMixingRatio(unit="g/kg", name="liquid water mixing ratio"),
 35        PySDM_products.PeakSaturation(name="S max"),
 36        PySDM_products.AmbientRelativeHumidity(name="RH"),
 37        PySDM_products.ParcelDisplacement(name="z"),
 38    )
 39
 40    formulae = Formulae(constants=CONSTANTS_ARG)
 41    const = formulae.constants
 42    pv0 = RH0 * formulae.saturation_vapour_pressure.pvs_water(T0)
 43
 44    env = Parcel(
 45        dt=dt,
 46        mass_of_dry_air=mass_of_dry_air,
 47        p0=p0,
 48        initial_water_vapour_mixing_ratio=const.eps * pv0 / (p0 - pv0),
 49        w=w,
 50        T0=T0,
 51        backend=CPU(formulae),
 52    )
 53
 54    aerosol = AerosolARG(
 55        M2_sol=sol2, M2_N=N2, M2_rad=rad2, water_molar_volume=const.Mv / const.rho_w
 56    )
 57    n_sd = n_sd_per_mode * len(aerosol.modes)
 58
 59    attributes = {
 60        k: np.empty(0) for k in ("dry volume", "kappa times dry volume", "multiplicity")
 61    }
 62    for i, mode in enumerate(aerosol.modes):
 63        kappa, spectrum = mode["kappa"]["CompressedFilmOvadnevaite"], mode["spectrum"]
 64        r_dry, concentration = ConstantMultiplicity(spectrum).sample_deterministic(
 65            n_sd_per_mode
 66        )
 67        v_dry = formulae.trivia.volume(radius=r_dry)
 68        specific_concentration = concentration / env.backend.formulae.constants.rho_STP
 69        attributes["multiplicity"] = np.append(
 70            attributes["multiplicity"],
 71            specific_concentration * env.mass_of_dry_air,
 72        )
 73        attributes["dry volume"] = np.append(attributes["dry volume"], v_dry)
 74        attributes["kappa times dry volume"] = np.append(
 75            attributes["kappa times dry volume"], v_dry * kappa
 76        )
 77
 78    r_wet = equilibrate_wet_radii(
 79        r_dry=formulae.trivia.radius(volume=attributes["dry volume"]),
 80        environment=env,
 81        kappa_times_dry_volume=attributes["kappa times dry volume"],
 82    )
 83    attributes["volume"] = env.backend.formulae.trivia.volume(radius=r_wet)
 84
 85    particulator = Particulator(
 86        n_sd=n_sd,
 87        environment=env,
 88        dynamics=[AmbientThermodynamics(), Condensation()],
 89        requested_attributes=("critical saturation",),
 90        attributes=attributes,
 91        products=products,
 92    )
 93
 94    output = {product.name: [] for product in particulator.products.values()}
 95    output_attributes = {
 96        "multiplicity": tuple([] for _ in range(particulator.n_sd)),
 97        "volume": tuple([] for _ in range(particulator.n_sd)),
 98        "critical volume": tuple([] for _ in range(particulator.n_sd)),
 99        "critical saturation": tuple([] for _ in range(particulator.n_sd)),
100    }
101
102    for _ in range(n_steps):
103        particulator.advance(n_steps=1)
104        for product in particulator.products.values():
105            value = product.get()
106            output[product.name].append(value[0])
107        for key, attr in output_attributes.items():
108            attr_data = particulator.attributes[key].to_ndarray()
109            for drop_id in range(particulator.n_sd):
110                attr[drop_id].append(attr_data[drop_id])
111
112    error = np.zeros(len(aerosol.modes))
113    activated_fraction_S = np.zeros(len(aerosol.modes))
114    activated_fraction_V = np.zeros(len(aerosol.modes))
115    for j, mode in enumerate(aerosol.modes):
116        activated_drops_j_S = 0
117        activated_drops_j_V = 0
118        RHmax = np.nanmax(np.asarray(output["RH"]))
119        for i, volume in enumerate(output_attributes["volume"]):
120            if j * n_sd_per_mode <= i < (j + 1) * n_sd_per_mode:
121                if output_attributes["critical saturation"][i][-1] < RHmax:
122                    activated_drops_j_S += output_attributes["multiplicity"][i][-1]
123                if output_attributes["critical volume"][i][-1] < volume[-1]:
124                    activated_drops_j_V += output_attributes["multiplicity"][i][-1]
125        Nj = np.asarray(output_attributes["multiplicity"])[
126            j * n_sd_per_mode : (j + 1) * n_sd_per_mode, -1
127        ]
128        max_multiplicity_j = np.max(Nj)
129        sum_multiplicity_j = np.sum(Nj)
130        error[j] = max_multiplicity_j / sum_multiplicity_j
131        activated_fraction_S[j] = activated_drops_j_S / sum_multiplicity_j
132        activated_fraction_V[j] = activated_drops_j_V / sum_multiplicity_j
133
134    Output = namedtuple(
135        "Output",
136        [
137            "profile",
138            "attributes",
139            "aerosol",
140            "activated_fraction_S",
141            "activated_fraction_V",
142            "error",
143        ],
144    )
145    return Output(
146        profile=output,
147        attributes=output_attributes,
148        aerosol=aerosol,
149        activated_fraction_S=activated_fraction_S,
150        activated_fraction_V=activated_fraction_V,
151        error=error,
152    )