PySDM_examples.deJong_Mackay_et_al_2023.simulation_ss
1import pickle as pkl 2 3import numpy as np 4from PySDM_examples.deJong_Mackay_et_al_2023.settings_0D import Settings0D 5from PySDM_examples.deJong_Mackay_et_al_2023.simulation_0D import run_box_breakup 6 7from PySDM.dynamics.collisions.breakup_fragmentations import Straub2010Nf 8from PySDM.dynamics.collisions.coalescence_efficiencies import Straub2010Ec 9from PySDM.initialisation.spectra import Exponential 10from PySDM.physics import si 11 12 13def run_to_steady_state(parameterization, n_sd, steps, nruns=1, dt=1 * si.s): 14 rain_rate = 54 * si.mm / si.h 15 mp_scale = 4.1e3 * (rain_rate / si.mm * si.h) ** (-0.21) / si.m 16 mp_N0 = 8e6 / si.m**4 17 n_part = mp_N0 / mp_scale 18 19 nbins = 64 20 21 y_ensemble = np.zeros((nruns, len(steps), nbins - 1)) 22 y2_ensemble = np.zeros((nruns, len(steps), nbins - 1)) 23 irun = 0 24 25 while irun < nruns: 26 if parameterization == "Straub2010": 27 settings = Settings0D( 28 seed=7 ** (irun + 1), 29 fragmentation=Straub2010Nf( 30 vmin=(0.01 * si.mm) ** 3 * np.pi / 6, 31 nfmax=10000, 32 ), 33 ) 34 settings.coal_eff = Straub2010Ec() 35 else: 36 print("parameterization not recognized") 37 return 38 39 settings.dv = 1e6 * si.m**3 40 settings.dt = dt 41 settings.spectrum = Exponential( 42 norm_factor=n_part * settings.dv, scale=1 / mp_scale 43 ) 44 settings.n_sd = n_sd 45 settings.radius_bins_edges = np.linspace( 46 1e-3 * si.mm, 4 * si.mm, num=nbins, endpoint=True 47 ) 48 49 settings.warn_overflows = False 50 settings._steps = steps # pylint: disable=protected-access 51 try: 52 res = run_box_breakup(settings, sample_in_radius=True, return_nv=True) 53 x, y, y2, rates = res.x, res.y, res.y2, res.rates 54 y_ensemble[irun] = y 55 y2_ensemble[irun] = y2 56 print("Success with run #" + str(irun + 1)) 57 irun += 1 58 except: # pylint: disable=bare-except 59 if dt > 0.5 * si.s: 60 print( 61 "Error in steady state sim for " 62 + str(n_sd) 63 + " superdroplets, moving on with dt=" 64 + str(dt / 2) 65 ) 66 dt = dt / 2 67 else: 68 print( 69 "Error in steady state sim for " 70 + str(n_sd) 71 + " superdroplets, proceeding to next iteration" 72 ) 73 rates = np.nan 74 x = (settings.radius_bins_edges[:-1] / si.micrometres,)[0] 75 y_ensemble[irun] = np.ones((len(steps), nbins - 1)) * np.nan 76 irun += 1 77 dt = 1 * si.s 78 79 data_filename = "steadystate_" + parameterization + "_" + str(n_sd) + "sd.pkl" 80 81 with open(data_filename, "wb") as handle: 82 pkl.dump( 83 (x, y_ensemble, y2_ensemble, rates), handle, protocol=pkl.HIGHEST_PROTOCOL 84 ) 85 86 87def get_straub_fig10_init(): 88 rain_rate = 54 * si.mm / si.h 89 mp_scale = 4.1e3 * (rain_rate / si.mm * si.h) ** (-0.21) / si.m 90 mp_N0 = 8e6 / si.m**4 91 92 straub_x_init = np.linspace(1e-3, 4.0, 100) * si.mm 93 straub_y_init = mp_N0 * np.exp(-1.0 * mp_scale * (straub_x_init)) * si.mm 94 95 dx = np.concatenate( 96 [np.diff(straub_x_init), [straub_x_init[-1] - straub_x_init[-2]]] 97 ) 98 x_c = straub_x_init + dx / 2 99 m_c = np.pi / 6 * x_c**3 100 straub_dvdlnr_init = m_c * straub_y_init / si.mm * x_c 101 102 return (straub_x_init, straub_y_init, straub_dvdlnr_init) 103 104 105def get_straub_fig10_data(): 106 graph_x = ( 107 np.array( 108 [ 109 0.08988764, 110 0.086142322, 111 0.097378277, 112 0.108614232, 113 0.119850187, 114 0.142322097, 115 0.164794007, 116 0.194756554, 117 0.224719101, 118 0.262172285, 119 0.314606742, 120 0.36329588, 121 0.419475655, 122 0.479400749, 123 0.558052434, 124 0.68164794, 125 0.816479401, 126 0.943820225, 127 1.071161049, 128 1.213483146, 129 1.370786517, 130 1.617977528, 131 1.865168539, 132 2.074906367, 133 2.322097378, 134 2.546816479, 135 2.801498127, 136 3.018726592, 137 3.220973783, 138 3.378277154, 139 3.543071161, 140 3.651685393, 141 ] 142 ) 143 * si.mm 144 ) 145 graph_log_y = np.array( 146 [ 147 1.055803251, 148 1.003199334, 149 1.14357294, 150 1.316140369, 151 1.509917836, 152 1.811447329, 153 2.159352957, 154 2.607176908, 155 3.048912888, 156 3.453402358, 157 3.849578422, 158 3.955529874, 159 3.849578422, 160 3.634485506, 161 3.328848104, 162 2.897005075, 163 2.525213782, 164 2.294463222, 165 2.125139584, 166 2.045222871, 167 2.02571776, 168 2.032198707, 169 2, 170 1.930947254, 171 1.811447329, 172 1.685826681, 173 1.531778159, 174 1.389585281, 175 1.262606891, 176 1.169430766, 177 1.069379699, 178 1.0064089029687935, 179 ] 180 ) 181 182 dx = np.concatenate([np.diff(graph_x), [graph_x[-1] - graph_x[-2]]]) 183 lnr = np.log(graph_x / 2) 184 dlnr = np.concatenate([np.diff(lnr), [lnr[-1] - lnr[-2]]]) 185 x_c = graph_x + dx / 2 186 m_c = np.pi / 6 * x_c**3 187 n = np.power(10.0, graph_log_y) * dx 188 straub_dvdlnr_ss = n * m_c / dlnr 189 190 return (graph_x, graph_log_y, straub_dvdlnr_ss)
def
run_to_steady_state(parameterization, n_sd, steps, nruns=1, dt=1.0):
14def run_to_steady_state(parameterization, n_sd, steps, nruns=1, dt=1 * si.s): 15 rain_rate = 54 * si.mm / si.h 16 mp_scale = 4.1e3 * (rain_rate / si.mm * si.h) ** (-0.21) / si.m 17 mp_N0 = 8e6 / si.m**4 18 n_part = mp_N0 / mp_scale 19 20 nbins = 64 21 22 y_ensemble = np.zeros((nruns, len(steps), nbins - 1)) 23 y2_ensemble = np.zeros((nruns, len(steps), nbins - 1)) 24 irun = 0 25 26 while irun < nruns: 27 if parameterization == "Straub2010": 28 settings = Settings0D( 29 seed=7 ** (irun + 1), 30 fragmentation=Straub2010Nf( 31 vmin=(0.01 * si.mm) ** 3 * np.pi / 6, 32 nfmax=10000, 33 ), 34 ) 35 settings.coal_eff = Straub2010Ec() 36 else: 37 print("parameterization not recognized") 38 return 39 40 settings.dv = 1e6 * si.m**3 41 settings.dt = dt 42 settings.spectrum = Exponential( 43 norm_factor=n_part * settings.dv, scale=1 / mp_scale 44 ) 45 settings.n_sd = n_sd 46 settings.radius_bins_edges = np.linspace( 47 1e-3 * si.mm, 4 * si.mm, num=nbins, endpoint=True 48 ) 49 50 settings.warn_overflows = False 51 settings._steps = steps # pylint: disable=protected-access 52 try: 53 res = run_box_breakup(settings, sample_in_radius=True, return_nv=True) 54 x, y, y2, rates = res.x, res.y, res.y2, res.rates 55 y_ensemble[irun] = y 56 y2_ensemble[irun] = y2 57 print("Success with run #" + str(irun + 1)) 58 irun += 1 59 except: # pylint: disable=bare-except 60 if dt > 0.5 * si.s: 61 print( 62 "Error in steady state sim for " 63 + str(n_sd) 64 + " superdroplets, moving on with dt=" 65 + str(dt / 2) 66 ) 67 dt = dt / 2 68 else: 69 print( 70 "Error in steady state sim for " 71 + str(n_sd) 72 + " superdroplets, proceeding to next iteration" 73 ) 74 rates = np.nan 75 x = (settings.radius_bins_edges[:-1] / si.micrometres,)[0] 76 y_ensemble[irun] = np.ones((len(steps), nbins - 1)) * np.nan 77 irun += 1 78 dt = 1 * si.s 79 80 data_filename = "steadystate_" + parameterization + "_" + str(n_sd) + "sd.pkl" 81 82 with open(data_filename, "wb") as handle: 83 pkl.dump( 84 (x, y_ensemble, y2_ensemble, rates), handle, protocol=pkl.HIGHEST_PROTOCOL 85 )
def
get_straub_fig10_init():
88def get_straub_fig10_init(): 89 rain_rate = 54 * si.mm / si.h 90 mp_scale = 4.1e3 * (rain_rate / si.mm * si.h) ** (-0.21) / si.m 91 mp_N0 = 8e6 / si.m**4 92 93 straub_x_init = np.linspace(1e-3, 4.0, 100) * si.mm 94 straub_y_init = mp_N0 * np.exp(-1.0 * mp_scale * (straub_x_init)) * si.mm 95 96 dx = np.concatenate( 97 [np.diff(straub_x_init), [straub_x_init[-1] - straub_x_init[-2]]] 98 ) 99 x_c = straub_x_init + dx / 2 100 m_c = np.pi / 6 * x_c**3 101 straub_dvdlnr_init = m_c * straub_y_init / si.mm * x_c 102 103 return (straub_x_init, straub_y_init, straub_dvdlnr_init)
def
get_straub_fig10_data():
106def get_straub_fig10_data(): 107 graph_x = ( 108 np.array( 109 [ 110 0.08988764, 111 0.086142322, 112 0.097378277, 113 0.108614232, 114 0.119850187, 115 0.142322097, 116 0.164794007, 117 0.194756554, 118 0.224719101, 119 0.262172285, 120 0.314606742, 121 0.36329588, 122 0.419475655, 123 0.479400749, 124 0.558052434, 125 0.68164794, 126 0.816479401, 127 0.943820225, 128 1.071161049, 129 1.213483146, 130 1.370786517, 131 1.617977528, 132 1.865168539, 133 2.074906367, 134 2.322097378, 135 2.546816479, 136 2.801498127, 137 3.018726592, 138 3.220973783, 139 3.378277154, 140 3.543071161, 141 3.651685393, 142 ] 143 ) 144 * si.mm 145 ) 146 graph_log_y = np.array( 147 [ 148 1.055803251, 149 1.003199334, 150 1.14357294, 151 1.316140369, 152 1.509917836, 153 1.811447329, 154 2.159352957, 155 2.607176908, 156 3.048912888, 157 3.453402358, 158 3.849578422, 159 3.955529874, 160 3.849578422, 161 3.634485506, 162 3.328848104, 163 2.897005075, 164 2.525213782, 165 2.294463222, 166 2.125139584, 167 2.045222871, 168 2.02571776, 169 2.032198707, 170 2, 171 1.930947254, 172 1.811447329, 173 1.685826681, 174 1.531778159, 175 1.389585281, 176 1.262606891, 177 1.169430766, 178 1.069379699, 179 1.0064089029687935, 180 ] 181 ) 182 183 dx = np.concatenate([np.diff(graph_x), [graph_x[-1] - graph_x[-2]]]) 184 lnr = np.log(graph_x / 2) 185 dlnr = np.concatenate([np.diff(lnr), [lnr[-1] - lnr[-2]]]) 186 x_c = graph_x + dx / 2 187 m_c = np.pi / 6 * x_c**3 188 n = np.power(10.0, graph_log_y) * dx 189 straub_dvdlnr_ss = n * m_c / dlnr 190 191 return (graph_x, graph_log_y, straub_dvdlnr_ss)