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