PySDM_examples.Alpert_and_Knopf_2016.simulation
1from typing import Union 2 3import matplotlib 4import numpy as np 5from matplotlib import pyplot 6from packaging import version 7 8from PySDM import Builder, Formulae 9from PySDM.backends import CPU 10from PySDM.dynamics import Freezing 11from PySDM.environments import Box 12from PySDM.initialisation import discretise_multiplicities 13from PySDM.initialisation.sampling import spectral_sampling 14from PySDM.physics import si 15from PySDM.products import IceWaterContent, TotalUnfrozenImmersedSurfaceArea 16 17 18class Simulation: 19 # note: dv and droplet_volume are dummy multipliers (multiplied and then divided by) 20 # will become used if coalescence or other processes are turned on 21 def __init__( 22 self, 23 *, 24 cases, 25 n_runs_per_case=10, 26 multiplicity=1, 27 time_step, 28 droplet_volume=1 * si.um**3, 29 heterogeneous_ice_nucleation_rate="Constant", 30 total_time: Union[None, float] = None, 31 temperature_range: Union[None, tuple] = None, 32 ): 33 self.cases = cases 34 self.n_runs_per_case = n_runs_per_case 35 self.multiplicity = multiplicity 36 self.volume = cases.volume 37 self.time_step = time_step 38 self.droplet_volume = droplet_volume 39 self.heterogeneous_ice_nucleation_rate = heterogeneous_ice_nucleation_rate 40 self.output = None 41 self.total_time = total_time 42 self.temperature_range = temperature_range 43 44 def run(self, keys): 45 self.output = {} 46 for key in keys: 47 case = self.cases[key] 48 49 assert (self.total_time is None) + (self.temperature_range is None) == 1 50 if self.total_time is not None: 51 total_time = self.total_time 52 else: 53 total_time = ( 54 np.diff(np.asarray(self.temperature_range)) / case["cooling_rate"] 55 ) 56 57 constants = None 58 if "J_het" not in case: 59 case["J_het"] = None 60 constants = {"ABIFM_C": case["ABIFM_c"], "ABIFM_M": case["ABIFM_m"]} 61 if "cooling_rate" not in case: 62 case["cooling_rate"] = 0 63 constants = {"J_HET": case["J_het"]} 64 65 self.output[key] = [] 66 for i in range(self.n_runs_per_case): 67 number_of_real_droplets = case["ISA"].norm_factor * self.volume 68 n_sd = number_of_real_droplets / self.multiplicity 69 np.testing.assert_approx_equal(n_sd, int(n_sd)) 70 n_sd = int(n_sd) 71 initial_temp = ( 72 self.temperature_range[1] if self.temperature_range else np.nan 73 ) 74 f_ufz, a_tot = simulation( 75 constants=constants, 76 seed=i, 77 n_sd=n_sd, 78 time_step=self.time_step, 79 volume=self.volume, 80 spectrum=case["ISA"], 81 droplet_volume=self.droplet_volume, 82 multiplicity=self.multiplicity, 83 total_time=total_time, 84 number_of_real_droplets=number_of_real_droplets, 85 cooling_rate=self.cases[key]["cooling_rate"], 86 heterogeneous_ice_nucleation_rate=self.heterogeneous_ice_nucleation_rate, 87 initial_temperature=initial_temp, 88 ) 89 self.output[key].append({"f_ufz": f_ufz, "A_tot": a_tot}) 90 91 def plot(self, ylim, grid=None): 92 pyplot.rc("font", size=10) 93 for key in self.output: 94 for run in range(self.n_runs_per_case): 95 time = self.time_step * np.arange(len(self.output[key][run]["f_ufz"])) 96 if self.cases[key]["cooling_rate"] == 0: 97 plot_x = time / si.min 98 plot_y = self.output[key][run]["f_ufz"] 99 else: 100 plot_x = ( 101 self.temperature_range[1] 102 - time * self.cases[key]["cooling_rate"] 103 ) 104 plot_y = 1 - np.asarray(self.output[key][run]["f_ufz"]) 105 pyplot.step( 106 plot_x, 107 plot_y, 108 label=self.cases.label(key) if run == 0 else None, 109 color=self.cases[key]["color"], 110 linewidth=0.666, 111 ) 112 key = None 113 if version.parse(matplotlib.__version__) >= version.parse("3.3.0"): 114 pyplot.gca().set_box_aspect(1) 115 pyplot.legend() 116 if grid is not None: 117 pyplot.grid(which=grid) 118 pyplot.ylim(ylim) 119 if self.temperature_range: 120 pyplot.xlim(*self.temperature_range) 121 pyplot.xlabel("T / K") 122 pyplot.ylabel("$f_{frz}$") 123 else: 124 pyplot.xlim(0, self.total_time / si.min) 125 pyplot.xlabel("t / min") 126 pyplot.ylabel("$f_{ufz}$") 127 pyplot.yscale("log") 128 129 def plot_j_het(self, variant: str, abifm_params_case: str, ylim=None): 130 assert variant in ("apparent", "actual") 131 132 formulae = Formulae( 133 particle_shape_and_density="MixedPhaseSpheres", 134 heterogeneous_ice_nucleation_rate="ABIFM", 135 constants={ 136 "ABIFM_M": self.cases[abifm_params_case]["ABIFM_m"], 137 "ABIFM_C": self.cases[abifm_params_case]["ABIFM_c"], 138 }, 139 ) 140 141 yunit = 1 / si.cm**2 / si.s 142 svp = formulae.saturation_vapour_pressure 143 plot_x = np.linspace(*self.temperature_range) * si.K 144 plot_y = formulae.heterogeneous_ice_nucleation_rate.j_het( 145 svp.pvs_ice(plot_x) / svp.pvs_water(plot_x) 146 ) 147 pyplot.grid() 148 pyplot.plot(plot_x, plot_y / yunit, color="red", label="ABIFM $J_{het}$") 149 150 for key in self.output: 151 for run in range(self.n_runs_per_case): 152 time = self.time_step * np.arange(len(self.output[key][run]["f_ufz"])) 153 if self.cases[key]["cooling_rate"] == 0: 154 raise NotImplementedError() 155 156 temperature = ( 157 self.temperature_range[1] - time * self.cases[key]["cooling_rate"] 158 ) 159 spec = self.cases[key]["ISA"] 160 161 particle_number = spec.norm_factor * self.volume 162 n_ufz = particle_number * np.asarray(self.output[key][run]["f_ufz"]) 163 n_frz = particle_number - n_ufz 164 165 j_het = np.diff(n_frz) / self.time_step 166 if variant == "apparent": 167 j_het /= n_ufz[:-1] * spec.m_mode 168 else: 169 a_tot = np.asarray(self.output[key][run]["A_tot"][:-1]) 170 j_het = np.divide( 171 j_het, a_tot, out=np.zeros_like(j_het), where=a_tot != 0 172 ) 173 174 pyplot.scatter( 175 temperature[:-1] + np.diff(temperature) / 2, 176 np.where(j_het != 0, j_het, np.nan) / yunit, 177 label=self.cases.label(key) if run == 0 else None, 178 color=self.cases[key]["color"], 179 ) 180 key = None 181 182 pyplot.yscale("log") 183 pyplot.xlabel("K") 184 pyplot.ylabel( 185 f"$J_{{het}}$, $J_{{het}}^{{{variant}}}$ / cm$^{{-2}}$ s$^{{-1}}$" 186 ) 187 pyplot.xlim(self.temperature_range) 188 if ylim is not None: 189 pyplot.ylim(ylim) 190 pyplot.legend() 191 if version.parse(matplotlib.__version__) >= version.parse("3.3.0"): 192 pyplot.gca().set_box_aspect(1) 193 194 195def simulation( 196 *, 197 constants, 198 seed, 199 n_sd, 200 time_step, 201 volume, 202 spectrum, 203 droplet_volume, 204 multiplicity, 205 total_time, 206 number_of_real_droplets, 207 cooling_rate=0, 208 heterogeneous_ice_nucleation_rate="Constant", 209 initial_temperature=np.nan, 210): 211 formulae = Formulae( 212 seed=seed, 213 heterogeneous_ice_nucleation_rate=heterogeneous_ice_nucleation_rate, 214 constants=constants, 215 particle_shape_and_density="MixedPhaseSpheres", 216 ) 217 builder = Builder( 218 n_sd=n_sd, 219 backend=CPU(formulae=formulae), 220 environment=Box(dt=time_step, dv=volume), 221 ) 222 builder.add_dynamic(Freezing(singular=False)) 223 224 if hasattr(spectrum, "s_geom") and spectrum.s_geom == 1: 225 _isa, _conc = np.full(n_sd, spectrum.m_mode), np.full( 226 n_sd, multiplicity / volume 227 ) 228 else: 229 _isa, _conc = spectral_sampling.ConstantMultiplicity(spectrum).sample(n_sd) 230 attributes = { 231 "multiplicity": discretise_multiplicities(_conc * volume), 232 "immersed surface area": _isa, 233 "volume": np.full(n_sd, droplet_volume), 234 } 235 np.testing.assert_almost_equal(attributes["multiplicity"], multiplicity) 236 products = ( 237 IceWaterContent(name="qi"), 238 TotalUnfrozenImmersedSurfaceArea(name="A_tot"), 239 ) 240 particulator = builder.build(attributes=attributes, products=products) 241 env = particulator.environment 242 243 env["T"] = initial_temperature 244 env["a_w_ice"] = np.nan 245 env["RH"] = 1 + np.finfo(float).eps 246 svp = particulator.formulae.saturation_vapour_pressure 247 248 cell_id = 0 249 f_ufz = [] 250 a_tot = [] 251 for i in range(int(total_time / time_step) + 1): 252 if cooling_rate != 0: 253 env["T"] -= np.full((1,), cooling_rate * time_step / 2) 254 env["a_w_ice"] = svp.pvs_ice(env["T"][0]) / svp.pvs_water(env["T"][0]) 255 particulator.run(0 if i == 0 else 1) 256 if cooling_rate != 0: 257 env["T"] -= np.full((1,), cooling_rate * time_step / 2) 258 259 ice_mass_per_volume = particulator.products["qi"].get()[cell_id] 260 ice_mass = ice_mass_per_volume * volume 261 ice_number = ice_mass / (formulae.constants.rho_w * droplet_volume) 262 unfrozen_fraction = 1 - ice_number / number_of_real_droplets 263 f_ufz.append(unfrozen_fraction) 264 a_tot.append(particulator.products["A_tot"].get()[cell_id]) 265 return f_ufz, a_tot
class
Simulation:
19class Simulation: 20 # note: dv and droplet_volume are dummy multipliers (multiplied and then divided by) 21 # will become used if coalescence or other processes are turned on 22 def __init__( 23 self, 24 *, 25 cases, 26 n_runs_per_case=10, 27 multiplicity=1, 28 time_step, 29 droplet_volume=1 * si.um**3, 30 heterogeneous_ice_nucleation_rate="Constant", 31 total_time: Union[None, float] = None, 32 temperature_range: Union[None, tuple] = None, 33 ): 34 self.cases = cases 35 self.n_runs_per_case = n_runs_per_case 36 self.multiplicity = multiplicity 37 self.volume = cases.volume 38 self.time_step = time_step 39 self.droplet_volume = droplet_volume 40 self.heterogeneous_ice_nucleation_rate = heterogeneous_ice_nucleation_rate 41 self.output = None 42 self.total_time = total_time 43 self.temperature_range = temperature_range 44 45 def run(self, keys): 46 self.output = {} 47 for key in keys: 48 case = self.cases[key] 49 50 assert (self.total_time is None) + (self.temperature_range is None) == 1 51 if self.total_time is not None: 52 total_time = self.total_time 53 else: 54 total_time = ( 55 np.diff(np.asarray(self.temperature_range)) / case["cooling_rate"] 56 ) 57 58 constants = None 59 if "J_het" not in case: 60 case["J_het"] = None 61 constants = {"ABIFM_C": case["ABIFM_c"], "ABIFM_M": case["ABIFM_m"]} 62 if "cooling_rate" not in case: 63 case["cooling_rate"] = 0 64 constants = {"J_HET": case["J_het"]} 65 66 self.output[key] = [] 67 for i in range(self.n_runs_per_case): 68 number_of_real_droplets = case["ISA"].norm_factor * self.volume 69 n_sd = number_of_real_droplets / self.multiplicity 70 np.testing.assert_approx_equal(n_sd, int(n_sd)) 71 n_sd = int(n_sd) 72 initial_temp = ( 73 self.temperature_range[1] if self.temperature_range else np.nan 74 ) 75 f_ufz, a_tot = simulation( 76 constants=constants, 77 seed=i, 78 n_sd=n_sd, 79 time_step=self.time_step, 80 volume=self.volume, 81 spectrum=case["ISA"], 82 droplet_volume=self.droplet_volume, 83 multiplicity=self.multiplicity, 84 total_time=total_time, 85 number_of_real_droplets=number_of_real_droplets, 86 cooling_rate=self.cases[key]["cooling_rate"], 87 heterogeneous_ice_nucleation_rate=self.heterogeneous_ice_nucleation_rate, 88 initial_temperature=initial_temp, 89 ) 90 self.output[key].append({"f_ufz": f_ufz, "A_tot": a_tot}) 91 92 def plot(self, ylim, grid=None): 93 pyplot.rc("font", size=10) 94 for key in self.output: 95 for run in range(self.n_runs_per_case): 96 time = self.time_step * np.arange(len(self.output[key][run]["f_ufz"])) 97 if self.cases[key]["cooling_rate"] == 0: 98 plot_x = time / si.min 99 plot_y = self.output[key][run]["f_ufz"] 100 else: 101 plot_x = ( 102 self.temperature_range[1] 103 - time * self.cases[key]["cooling_rate"] 104 ) 105 plot_y = 1 - np.asarray(self.output[key][run]["f_ufz"]) 106 pyplot.step( 107 plot_x, 108 plot_y, 109 label=self.cases.label(key) if run == 0 else None, 110 color=self.cases[key]["color"], 111 linewidth=0.666, 112 ) 113 key = None 114 if version.parse(matplotlib.__version__) >= version.parse("3.3.0"): 115 pyplot.gca().set_box_aspect(1) 116 pyplot.legend() 117 if grid is not None: 118 pyplot.grid(which=grid) 119 pyplot.ylim(ylim) 120 if self.temperature_range: 121 pyplot.xlim(*self.temperature_range) 122 pyplot.xlabel("T / K") 123 pyplot.ylabel("$f_{frz}$") 124 else: 125 pyplot.xlim(0, self.total_time / si.min) 126 pyplot.xlabel("t / min") 127 pyplot.ylabel("$f_{ufz}$") 128 pyplot.yscale("log") 129 130 def plot_j_het(self, variant: str, abifm_params_case: str, ylim=None): 131 assert variant in ("apparent", "actual") 132 133 formulae = Formulae( 134 particle_shape_and_density="MixedPhaseSpheres", 135 heterogeneous_ice_nucleation_rate="ABIFM", 136 constants={ 137 "ABIFM_M": self.cases[abifm_params_case]["ABIFM_m"], 138 "ABIFM_C": self.cases[abifm_params_case]["ABIFM_c"], 139 }, 140 ) 141 142 yunit = 1 / si.cm**2 / si.s 143 svp = formulae.saturation_vapour_pressure 144 plot_x = np.linspace(*self.temperature_range) * si.K 145 plot_y = formulae.heterogeneous_ice_nucleation_rate.j_het( 146 svp.pvs_ice(plot_x) / svp.pvs_water(plot_x) 147 ) 148 pyplot.grid() 149 pyplot.plot(plot_x, plot_y / yunit, color="red", label="ABIFM $J_{het}$") 150 151 for key in self.output: 152 for run in range(self.n_runs_per_case): 153 time = self.time_step * np.arange(len(self.output[key][run]["f_ufz"])) 154 if self.cases[key]["cooling_rate"] == 0: 155 raise NotImplementedError() 156 157 temperature = ( 158 self.temperature_range[1] - time * self.cases[key]["cooling_rate"] 159 ) 160 spec = self.cases[key]["ISA"] 161 162 particle_number = spec.norm_factor * self.volume 163 n_ufz = particle_number * np.asarray(self.output[key][run]["f_ufz"]) 164 n_frz = particle_number - n_ufz 165 166 j_het = np.diff(n_frz) / self.time_step 167 if variant == "apparent": 168 j_het /= n_ufz[:-1] * spec.m_mode 169 else: 170 a_tot = np.asarray(self.output[key][run]["A_tot"][:-1]) 171 j_het = np.divide( 172 j_het, a_tot, out=np.zeros_like(j_het), where=a_tot != 0 173 ) 174 175 pyplot.scatter( 176 temperature[:-1] + np.diff(temperature) / 2, 177 np.where(j_het != 0, j_het, np.nan) / yunit, 178 label=self.cases.label(key) if run == 0 else None, 179 color=self.cases[key]["color"], 180 ) 181 key = None 182 183 pyplot.yscale("log") 184 pyplot.xlabel("K") 185 pyplot.ylabel( 186 f"$J_{{het}}$, $J_{{het}}^{{{variant}}}$ / cm$^{{-2}}$ s$^{{-1}}$" 187 ) 188 pyplot.xlim(self.temperature_range) 189 if ylim is not None: 190 pyplot.ylim(ylim) 191 pyplot.legend() 192 if version.parse(matplotlib.__version__) >= version.parse("3.3.0"): 193 pyplot.gca().set_box_aspect(1)
Simulation( *, cases, n_runs_per_case=10, multiplicity=1, time_step, droplet_volume=9.999999999999999e-19, heterogeneous_ice_nucleation_rate='Constant', total_time: Optional[float] = None, temperature_range: Optional[tuple] = None)
22 def __init__( 23 self, 24 *, 25 cases, 26 n_runs_per_case=10, 27 multiplicity=1, 28 time_step, 29 droplet_volume=1 * si.um**3, 30 heterogeneous_ice_nucleation_rate="Constant", 31 total_time: Union[None, float] = None, 32 temperature_range: Union[None, tuple] = None, 33 ): 34 self.cases = cases 35 self.n_runs_per_case = n_runs_per_case 36 self.multiplicity = multiplicity 37 self.volume = cases.volume 38 self.time_step = time_step 39 self.droplet_volume = droplet_volume 40 self.heterogeneous_ice_nucleation_rate = heterogeneous_ice_nucleation_rate 41 self.output = None 42 self.total_time = total_time 43 self.temperature_range = temperature_range
def
run(self, keys):
45 def run(self, keys): 46 self.output = {} 47 for key in keys: 48 case = self.cases[key] 49 50 assert (self.total_time is None) + (self.temperature_range is None) == 1 51 if self.total_time is not None: 52 total_time = self.total_time 53 else: 54 total_time = ( 55 np.diff(np.asarray(self.temperature_range)) / case["cooling_rate"] 56 ) 57 58 constants = None 59 if "J_het" not in case: 60 case["J_het"] = None 61 constants = {"ABIFM_C": case["ABIFM_c"], "ABIFM_M": case["ABIFM_m"]} 62 if "cooling_rate" not in case: 63 case["cooling_rate"] = 0 64 constants = {"J_HET": case["J_het"]} 65 66 self.output[key] = [] 67 for i in range(self.n_runs_per_case): 68 number_of_real_droplets = case["ISA"].norm_factor * self.volume 69 n_sd = number_of_real_droplets / self.multiplicity 70 np.testing.assert_approx_equal(n_sd, int(n_sd)) 71 n_sd = int(n_sd) 72 initial_temp = ( 73 self.temperature_range[1] if self.temperature_range else np.nan 74 ) 75 f_ufz, a_tot = simulation( 76 constants=constants, 77 seed=i, 78 n_sd=n_sd, 79 time_step=self.time_step, 80 volume=self.volume, 81 spectrum=case["ISA"], 82 droplet_volume=self.droplet_volume, 83 multiplicity=self.multiplicity, 84 total_time=total_time, 85 number_of_real_droplets=number_of_real_droplets, 86 cooling_rate=self.cases[key]["cooling_rate"], 87 heterogeneous_ice_nucleation_rate=self.heterogeneous_ice_nucleation_rate, 88 initial_temperature=initial_temp, 89 ) 90 self.output[key].append({"f_ufz": f_ufz, "A_tot": a_tot})
def
plot(self, ylim, grid=None):
92 def plot(self, ylim, grid=None): 93 pyplot.rc("font", size=10) 94 for key in self.output: 95 for run in range(self.n_runs_per_case): 96 time = self.time_step * np.arange(len(self.output[key][run]["f_ufz"])) 97 if self.cases[key]["cooling_rate"] == 0: 98 plot_x = time / si.min 99 plot_y = self.output[key][run]["f_ufz"] 100 else: 101 plot_x = ( 102 self.temperature_range[1] 103 - time * self.cases[key]["cooling_rate"] 104 ) 105 plot_y = 1 - np.asarray(self.output[key][run]["f_ufz"]) 106 pyplot.step( 107 plot_x, 108 plot_y, 109 label=self.cases.label(key) if run == 0 else None, 110 color=self.cases[key]["color"], 111 linewidth=0.666, 112 ) 113 key = None 114 if version.parse(matplotlib.__version__) >= version.parse("3.3.0"): 115 pyplot.gca().set_box_aspect(1) 116 pyplot.legend() 117 if grid is not None: 118 pyplot.grid(which=grid) 119 pyplot.ylim(ylim) 120 if self.temperature_range: 121 pyplot.xlim(*self.temperature_range) 122 pyplot.xlabel("T / K") 123 pyplot.ylabel("$f_{frz}$") 124 else: 125 pyplot.xlim(0, self.total_time / si.min) 126 pyplot.xlabel("t / min") 127 pyplot.ylabel("$f_{ufz}$") 128 pyplot.yscale("log")
def
plot_j_het(self, variant: str, abifm_params_case: str, ylim=None):
130 def plot_j_het(self, variant: str, abifm_params_case: str, ylim=None): 131 assert variant in ("apparent", "actual") 132 133 formulae = Formulae( 134 particle_shape_and_density="MixedPhaseSpheres", 135 heterogeneous_ice_nucleation_rate="ABIFM", 136 constants={ 137 "ABIFM_M": self.cases[abifm_params_case]["ABIFM_m"], 138 "ABIFM_C": self.cases[abifm_params_case]["ABIFM_c"], 139 }, 140 ) 141 142 yunit = 1 / si.cm**2 / si.s 143 svp = formulae.saturation_vapour_pressure 144 plot_x = np.linspace(*self.temperature_range) * si.K 145 plot_y = formulae.heterogeneous_ice_nucleation_rate.j_het( 146 svp.pvs_ice(plot_x) / svp.pvs_water(plot_x) 147 ) 148 pyplot.grid() 149 pyplot.plot(plot_x, plot_y / yunit, color="red", label="ABIFM $J_{het}$") 150 151 for key in self.output: 152 for run in range(self.n_runs_per_case): 153 time = self.time_step * np.arange(len(self.output[key][run]["f_ufz"])) 154 if self.cases[key]["cooling_rate"] == 0: 155 raise NotImplementedError() 156 157 temperature = ( 158 self.temperature_range[1] - time * self.cases[key]["cooling_rate"] 159 ) 160 spec = self.cases[key]["ISA"] 161 162 particle_number = spec.norm_factor * self.volume 163 n_ufz = particle_number * np.asarray(self.output[key][run]["f_ufz"]) 164 n_frz = particle_number - n_ufz 165 166 j_het = np.diff(n_frz) / self.time_step 167 if variant == "apparent": 168 j_het /= n_ufz[:-1] * spec.m_mode 169 else: 170 a_tot = np.asarray(self.output[key][run]["A_tot"][:-1]) 171 j_het = np.divide( 172 j_het, a_tot, out=np.zeros_like(j_het), where=a_tot != 0 173 ) 174 175 pyplot.scatter( 176 temperature[:-1] + np.diff(temperature) / 2, 177 np.where(j_het != 0, j_het, np.nan) / yunit, 178 label=self.cases.label(key) if run == 0 else None, 179 color=self.cases[key]["color"], 180 ) 181 key = None 182 183 pyplot.yscale("log") 184 pyplot.xlabel("K") 185 pyplot.ylabel( 186 f"$J_{{het}}$, $J_{{het}}^{{{variant}}}$ / cm$^{{-2}}$ s$^{{-1}}$" 187 ) 188 pyplot.xlim(self.temperature_range) 189 if ylim is not None: 190 pyplot.ylim(ylim) 191 pyplot.legend() 192 if version.parse(matplotlib.__version__) >= version.parse("3.3.0"): 193 pyplot.gca().set_box_aspect(1)
def
simulation( *, constants, seed, n_sd, time_step, volume, spectrum, droplet_volume, multiplicity, total_time, number_of_real_droplets, cooling_rate=0, heterogeneous_ice_nucleation_rate='Constant', initial_temperature=nan):
196def simulation( 197 *, 198 constants, 199 seed, 200 n_sd, 201 time_step, 202 volume, 203 spectrum, 204 droplet_volume, 205 multiplicity, 206 total_time, 207 number_of_real_droplets, 208 cooling_rate=0, 209 heterogeneous_ice_nucleation_rate="Constant", 210 initial_temperature=np.nan, 211): 212 formulae = Formulae( 213 seed=seed, 214 heterogeneous_ice_nucleation_rate=heterogeneous_ice_nucleation_rate, 215 constants=constants, 216 particle_shape_and_density="MixedPhaseSpheres", 217 ) 218 builder = Builder( 219 n_sd=n_sd, 220 backend=CPU(formulae=formulae), 221 environment=Box(dt=time_step, dv=volume), 222 ) 223 builder.add_dynamic(Freezing(singular=False)) 224 225 if hasattr(spectrum, "s_geom") and spectrum.s_geom == 1: 226 _isa, _conc = np.full(n_sd, spectrum.m_mode), np.full( 227 n_sd, multiplicity / volume 228 ) 229 else: 230 _isa, _conc = spectral_sampling.ConstantMultiplicity(spectrum).sample(n_sd) 231 attributes = { 232 "multiplicity": discretise_multiplicities(_conc * volume), 233 "immersed surface area": _isa, 234 "volume": np.full(n_sd, droplet_volume), 235 } 236 np.testing.assert_almost_equal(attributes["multiplicity"], multiplicity) 237 products = ( 238 IceWaterContent(name="qi"), 239 TotalUnfrozenImmersedSurfaceArea(name="A_tot"), 240 ) 241 particulator = builder.build(attributes=attributes, products=products) 242 env = particulator.environment 243 244 env["T"] = initial_temperature 245 env["a_w_ice"] = np.nan 246 env["RH"] = 1 + np.finfo(float).eps 247 svp = particulator.formulae.saturation_vapour_pressure 248 249 cell_id = 0 250 f_ufz = [] 251 a_tot = [] 252 for i in range(int(total_time / time_step) + 1): 253 if cooling_rate != 0: 254 env["T"] -= np.full((1,), cooling_rate * time_step / 2) 255 env["a_w_ice"] = svp.pvs_ice(env["T"][0]) / svp.pvs_water(env["T"][0]) 256 particulator.run(0 if i == 0 else 1) 257 if cooling_rate != 0: 258 env["T"] -= np.full((1,), cooling_rate * time_step / 2) 259 260 ice_mass_per_volume = particulator.products["qi"].get()[cell_id] 261 ice_mass = ice_mass_per_volume * volume 262 ice_number = ice_mass / (formulae.constants.rho_w * droplet_volume) 263 unfrozen_fraction = 1 - ice_number / number_of_real_droplets 264 f_ufz.append(unfrozen_fraction) 265 a_tot.append(particulator.products["A_tot"].get()[cell_id]) 266 return f_ufz, a_tot