PySDM_examples.Abdul_Razzak_Ghan_2000.run_ARG_parcel

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