PySDM_examples.deJong_Mackay_et_al_2023.simulation_0D
1from collections import namedtuple 2 3import numpy as np 4 5from PySDM.backends import CPU 6from PySDM.builder import Builder 7from PySDM.dynamics import Coalescence, Collision 8from PySDM.environments import Box 9from PySDM.formulae import Formulae 10from PySDM.initialisation.sampling.spectral_sampling import ( 11 ConstantMultiplicity, 12 Logarithmic, 13) 14from PySDM.physics import si 15from PySDM.products.collision.collision_rates import ( 16 BreakupRatePerGridbox, 17 CoalescenceRatePerGridbox, 18 CollisionRateDeficitPerGridbox, 19 CollisionRatePerGridbox, 20) 21from PySDM.products.size_spectral import ( 22 NumberSizeSpectrum, 23 ParticleVolumeVersusRadiusLogarithmSpectrum, 24) 25 26 27def run_box_breakup( 28 settings, steps=None, backend_class=CPU, sample_in_radius=False, return_nv=False 29): 30 builder = Builder( 31 n_sd=settings.n_sd, 32 backend=backend_class(settings.formulae), 33 environment=Box(dv=settings.dv, dt=settings.dt), 34 ) 35 builder.particulator.environment["rhod"] = 1.0 36 attributes = {} 37 if sample_in_radius: 38 diams, attributes["multiplicity"] = Logarithmic(settings.spectrum).sample( 39 settings.n_sd 40 ) 41 radii = diams / 2 42 attributes["volume"] = Formulae().trivia.volume(radius=radii) 43 else: 44 attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity( 45 settings.spectrum 46 ).sample(settings.n_sd) 47 breakup = Collision( 48 collision_kernel=settings.kernel, 49 coalescence_efficiency=settings.coal_eff, 50 breakup_efficiency=settings.break_eff, 51 fragmentation_function=settings.fragmentation, 52 adaptive=settings.adaptive, 53 warn_overflows=settings.warn_overflows, 54 ) 55 builder.add_dynamic(breakup) 56 products = ( 57 ParticleVolumeVersusRadiusLogarithmSpectrum( 58 radius_bins_edges=settings.radius_bins_edges, name="dv/dlnr" 59 ), 60 NumberSizeSpectrum(radius_bins_edges=settings.radius_bins_edges, name="N(v)"), 61 CollisionRatePerGridbox(name="cr"), 62 CollisionRateDeficitPerGridbox(name="crd"), 63 CoalescenceRatePerGridbox(name="cor"), 64 BreakupRatePerGridbox(name="br"), 65 ) 66 core = builder.build(attributes, products) 67 68 if steps is None: 69 steps = settings.output_steps 70 y = np.ndarray((len(steps), len(settings.radius_bins_edges) - 1)) 71 if return_nv: 72 y2 = np.ndarray((len(steps), len(settings.radius_bins_edges) - 1)) 73 else: 74 y2 = None 75 76 rates = np.zeros((len(steps), 4)) 77 # run 78 for i, step in enumerate(steps): 79 core.run(step - core.n_steps) 80 y[i] = core.products["dv/dlnr"].get()[0] 81 if return_nv: 82 (y2[i],) = core.products["N(v)"].get() 83 (rates[i, 0],) = core.products["cr"].get() 84 (rates[i, 1],) = core.products["crd"].get() 85 (rates[i, 2],) = core.products["cor"].get() 86 (rates[i, 3],) = core.products["br"].get() 87 88 x = (settings.radius_bins_edges[:-1] / si.micrometres,)[0] 89 90 return namedtuple("_", ("x", "y", "y2", "rates"))(x=x, y=y, y2=y2, rates=rates) 91 92 93def run_box_NObreakup(settings, steps=None, backend_class=CPU): 94 builder = Builder( 95 n_sd=settings.n_sd, 96 backend=backend_class(settings.formulae), 97 environment=Box(dv=settings.dv, dt=settings.dt), 98 ) 99 builder.particulator.environment["rhod"] = 1.0 100 attributes = {} 101 attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity( 102 settings.spectrum 103 ).sample(settings.n_sd) 104 coal = Coalescence( 105 collision_kernel=settings.kernel, 106 coalescence_efficiency=settings.coal_eff, 107 adaptive=settings.adaptive, 108 ) 109 builder.add_dynamic(coal) 110 products = ( 111 ParticleVolumeVersusRadiusLogarithmSpectrum( 112 radius_bins_edges=settings.radius_bins_edges, name="dv/dlnr" 113 ), 114 CollisionRatePerGridbox(name="cr"), 115 CollisionRateDeficitPerGridbox(name="crd"), 116 CoalescenceRatePerGridbox(name="cor"), 117 ) 118 core = builder.build(attributes, products) 119 120 if steps is None: 121 steps = settings.output_steps 122 y = np.ndarray((len(steps), len(settings.radius_bins_edges) - 1)) 123 rates = np.zeros((len(steps), 4)) 124 # run 125 for i, step in enumerate(steps): 126 core.run(step - core.n_steps) 127 (y[i],) = core.products["dv/dlnr"].get() 128 (rates[i, 0],) = core.products["cr"].get() 129 (rates[i, 1],) = core.products["crd"].get() 130 (rates[i, 2],) = core.products["cor"].get() 131 132 x = (settings.radius_bins_edges[:-1] / si.micrometres,)[0] 133 134 return (x, y, rates)
def
run_box_breakup( settings, steps=None, backend_class=<class 'PySDM.backends.numba.Numba'>, sample_in_radius=False, return_nv=False):
28def run_box_breakup( 29 settings, steps=None, backend_class=CPU, sample_in_radius=False, return_nv=False 30): 31 builder = Builder( 32 n_sd=settings.n_sd, 33 backend=backend_class(settings.formulae), 34 environment=Box(dv=settings.dv, dt=settings.dt), 35 ) 36 builder.particulator.environment["rhod"] = 1.0 37 attributes = {} 38 if sample_in_radius: 39 diams, attributes["multiplicity"] = Logarithmic(settings.spectrum).sample( 40 settings.n_sd 41 ) 42 radii = diams / 2 43 attributes["volume"] = Formulae().trivia.volume(radius=radii) 44 else: 45 attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity( 46 settings.spectrum 47 ).sample(settings.n_sd) 48 breakup = Collision( 49 collision_kernel=settings.kernel, 50 coalescence_efficiency=settings.coal_eff, 51 breakup_efficiency=settings.break_eff, 52 fragmentation_function=settings.fragmentation, 53 adaptive=settings.adaptive, 54 warn_overflows=settings.warn_overflows, 55 ) 56 builder.add_dynamic(breakup) 57 products = ( 58 ParticleVolumeVersusRadiusLogarithmSpectrum( 59 radius_bins_edges=settings.radius_bins_edges, name="dv/dlnr" 60 ), 61 NumberSizeSpectrum(radius_bins_edges=settings.radius_bins_edges, name="N(v)"), 62 CollisionRatePerGridbox(name="cr"), 63 CollisionRateDeficitPerGridbox(name="crd"), 64 CoalescenceRatePerGridbox(name="cor"), 65 BreakupRatePerGridbox(name="br"), 66 ) 67 core = builder.build(attributes, products) 68 69 if steps is None: 70 steps = settings.output_steps 71 y = np.ndarray((len(steps), len(settings.radius_bins_edges) - 1)) 72 if return_nv: 73 y2 = np.ndarray((len(steps), len(settings.radius_bins_edges) - 1)) 74 else: 75 y2 = None 76 77 rates = np.zeros((len(steps), 4)) 78 # run 79 for i, step in enumerate(steps): 80 core.run(step - core.n_steps) 81 y[i] = core.products["dv/dlnr"].get()[0] 82 if return_nv: 83 (y2[i],) = core.products["N(v)"].get() 84 (rates[i, 0],) = core.products["cr"].get() 85 (rates[i, 1],) = core.products["crd"].get() 86 (rates[i, 2],) = core.products["cor"].get() 87 (rates[i, 3],) = core.products["br"].get() 88 89 x = (settings.radius_bins_edges[:-1] / si.micrometres,)[0] 90 91 return namedtuple("_", ("x", "y", "y2", "rates"))(x=x, y=y, y2=y2, rates=rates)
94def run_box_NObreakup(settings, steps=None, backend_class=CPU): 95 builder = Builder( 96 n_sd=settings.n_sd, 97 backend=backend_class(settings.formulae), 98 environment=Box(dv=settings.dv, dt=settings.dt), 99 ) 100 builder.particulator.environment["rhod"] = 1.0 101 attributes = {} 102 attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity( 103 settings.spectrum 104 ).sample(settings.n_sd) 105 coal = Coalescence( 106 collision_kernel=settings.kernel, 107 coalescence_efficiency=settings.coal_eff, 108 adaptive=settings.adaptive, 109 ) 110 builder.add_dynamic(coal) 111 products = ( 112 ParticleVolumeVersusRadiusLogarithmSpectrum( 113 radius_bins_edges=settings.radius_bins_edges, name="dv/dlnr" 114 ), 115 CollisionRatePerGridbox(name="cr"), 116 CollisionRateDeficitPerGridbox(name="crd"), 117 CoalescenceRatePerGridbox(name="cor"), 118 ) 119 core = builder.build(attributes, products) 120 121 if steps is None: 122 steps = settings.output_steps 123 y = np.ndarray((len(steps), len(settings.radius_bins_edges) - 1)) 124 rates = np.zeros((len(steps), 4)) 125 # run 126 for i, step in enumerate(steps): 127 core.run(step - core.n_steps) 128 (y[i],) = core.products["dv/dlnr"].get() 129 (rates[i, 0],) = core.products["cr"].get() 130 (rates[i, 1],) = core.products["crd"].get() 131 (rates[i, 2],) = core.products["cor"].get() 132 133 x = (settings.radius_bins_edges[:-1] / si.micrometres,)[0] 134 135 return (x, y, rates)