PyMPDATA_examples.Olesik_et_al_2022.settings
1import numpy as np 2import pint 3from PyMPDATA_examples.Olesik_et_al_2022 import ( 4 East_and_Marshall_1954, 5 equilibrium_drop_growth, 6) 7from pystrict import strict 8from scipy import integrate, optimize 9 10default_nr = 64 11default_GC_max = 0.5 12default_mixing_ratios_g_kg = np.array([1, 2, 4, 6, 8, 10]) 13default_opt_set = { 14 "a": {"n_iters": 1}, 15 "b": {"n_iters": 2}, 16 "c": {"n_iters": 2, "infinite_gauge": True}, 17 "d": {"n_iters": 2, "infinite_gauge": True, "nonoscillatory": True}, 18 "e": {"n_iters": 2, "DPDC": True, "infinite_gauge": True, "nonoscillatory": True}, 19 "f": {"n_iters": 3, "third_order_terms": True}, 20 "g": {"n_iters": 3}, 21 "h": { 22 "n_iters": 3, 23 "third_order_terms": True, 24 "infinite_gauge": True, 25 "nonoscillatory": True, 26 }, 27} 28colors = ["red", "blue", "crimson", "orange", "olive", "navy", "green", "blueviolet"] 29colors = {key: colors.pop(0) for key in default_opt_set} 30 31 32def option_string(opts): 33 str_repl = [ 34 ["'n_iters': 1", "upwind"], 35 ["'n_iters': 2", "MPDATA 2 iterations"], 36 ["'n_iters': 3", "MPDATA 3 iterations"], 37 ["'", ""], 38 [": True", ""], 39 ["_", " "], 40 ["{", ""], 41 ["}", ""], 42 [",", ""], 43 ["nonoscillatory", "non-oscillatory"], 44 ] 45 for repl in str_repl: 46 opts = opts.replace(repl[0], repl[1]) 47 return opts 48 49 50si = pint.UnitRegistry() 51ksi_1 = 100 * si.micrometre**2 / si.second 52 53 54# based on Fig. 3 from East 1957 55@strict 56class Settings: 57 def __init__( 58 self, 59 nr: int = default_nr, 60 mixing_ratios_g_kg: np.ndarray = default_mixing_ratios_g_kg, 61 ): 62 self.si = si 63 self.nr = nr 64 self.r_min = 1 * si.micrometre 65 self.r_max = 25 * si.micrometre 66 self.rho_w = 1 * si.kilogram / si.decimetre**3 67 self.rho_a = 1 * si.kilogram / si.metre**3 68 self.mixing_ratios = mixing_ratios_g_kg * si.gram / si.kilogram 69 S = 1.00075 70 self.drdt = equilibrium_drop_growth.DrDt(ksi_1, S) 71 self.size_distribution = East_and_Marshall_1954.SizeDistribution(si) 72 73 self.C = (1 * si.gram / si.kilogram) / self.mixing_ratio( 74 self.size_distribution.pdf 75 ) 76 np.testing.assert_approx_equal(self.C.to(si.dimensionless).magnitude, 1, 4) 77 self.out_times = self.find_out_steps() 78 79 def find_out_steps(self): 80 out_steps = [] 81 for mr in self.mixing_ratios: 82 t_unit = self.si.second 83 84 def findroot(ti, mr=mr, t_unit=t_unit): 85 return ( 86 mr 87 - self.mixing_ratio( 88 equilibrium_drop_growth.PdfEvolver( 89 self.pdf, self.drdt, ti * t_unit 90 ) 91 ) 92 ).magnitude 93 94 t = optimize.brentq(findroot, 0, (1 * self.si.hour).to(t_unit).magnitude) 95 out_steps.append(t) 96 return out_steps 97 98 def mixing_ratio(self, pdf): 99 xunit = self.si.micrometre 100 yunit = 1 / self.si.micrometre / self.si.centimetre**3 101 102 r_min = 0.1 * self.si.um 103 while not np.isfinite(pdf(r_min).magnitude): 104 r_min *= 1.01 105 106 def pdfarg(r_nounit): 107 r = r_nounit * xunit 108 result = pdf(r) * r**3 109 return result.to(yunit * xunit**3).magnitude 110 111 I = ( 112 integrate.quad(pdfarg, r_min.to(xunit).magnitude, np.inf)[0] 113 * yunit 114 * xunit**4 115 ) 116 return (I * 4 / 3 * np.pi * self.rho_w / self.rho_a).to( 117 self.si.gram / self.si.kilogram 118 ) 119 120 def pdf(self, r): 121 return self.C * self.size_distribution.pdf(r)
default_nr =
64
default_GC_max =
0.5
default_mixing_ratios_g_kg =
array([ 1, 2, 4, 6, 8, 10])
default_opt_set =
{'a': {'n_iters': 1}, 'b': {'n_iters': 2}, 'c': {'n_iters': 2, 'infinite_gauge': True}, 'd': {'n_iters': 2, 'infinite_gauge': True, 'nonoscillatory': True}, 'e': {'n_iters': 2, 'DPDC': True, 'infinite_gauge': True, 'nonoscillatory': True}, 'f': {'n_iters': 3, 'third_order_terms': True}, 'g': {'n_iters': 3}, 'h': {'n_iters': 3, 'third_order_terms': True, 'infinite_gauge': True, 'nonoscillatory': True}}
colors =
{'a': 'red', 'b': 'blue', 'c': 'crimson', 'd': 'orange', 'e': 'olive', 'f': 'navy', 'g': 'green', 'h': 'blueviolet'}
def
option_string(opts):
33def option_string(opts): 34 str_repl = [ 35 ["'n_iters': 1", "upwind"], 36 ["'n_iters': 2", "MPDATA 2 iterations"], 37 ["'n_iters': 3", "MPDATA 3 iterations"], 38 ["'", ""], 39 [": True", ""], 40 ["_", " "], 41 ["{", ""], 42 ["}", ""], 43 [",", ""], 44 ["nonoscillatory", "non-oscillatory"], 45 ] 46 for repl in str_repl: 47 opts = opts.replace(repl[0], repl[1]) 48 return opts
si =
<pint.registry.UnitRegistry object>
ksi_1 =
<Quantity(100.0, 'micrometer ** 2 / second')>
@strict
class
Settings:
56@strict 57class Settings: 58 def __init__( 59 self, 60 nr: int = default_nr, 61 mixing_ratios_g_kg: np.ndarray = default_mixing_ratios_g_kg, 62 ): 63 self.si = si 64 self.nr = nr 65 self.r_min = 1 * si.micrometre 66 self.r_max = 25 * si.micrometre 67 self.rho_w = 1 * si.kilogram / si.decimetre**3 68 self.rho_a = 1 * si.kilogram / si.metre**3 69 self.mixing_ratios = mixing_ratios_g_kg * si.gram / si.kilogram 70 S = 1.00075 71 self.drdt = equilibrium_drop_growth.DrDt(ksi_1, S) 72 self.size_distribution = East_and_Marshall_1954.SizeDistribution(si) 73 74 self.C = (1 * si.gram / si.kilogram) / self.mixing_ratio( 75 self.size_distribution.pdf 76 ) 77 np.testing.assert_approx_equal(self.C.to(si.dimensionless).magnitude, 1, 4) 78 self.out_times = self.find_out_steps() 79 80 def find_out_steps(self): 81 out_steps = [] 82 for mr in self.mixing_ratios: 83 t_unit = self.si.second 84 85 def findroot(ti, mr=mr, t_unit=t_unit): 86 return ( 87 mr 88 - self.mixing_ratio( 89 equilibrium_drop_growth.PdfEvolver( 90 self.pdf, self.drdt, ti * t_unit 91 ) 92 ) 93 ).magnitude 94 95 t = optimize.brentq(findroot, 0, (1 * self.si.hour).to(t_unit).magnitude) 96 out_steps.append(t) 97 return out_steps 98 99 def mixing_ratio(self, pdf): 100 xunit = self.si.micrometre 101 yunit = 1 / self.si.micrometre / self.si.centimetre**3 102 103 r_min = 0.1 * self.si.um 104 while not np.isfinite(pdf(r_min).magnitude): 105 r_min *= 1.01 106 107 def pdfarg(r_nounit): 108 r = r_nounit * xunit 109 result = pdf(r) * r**3 110 return result.to(yunit * xunit**3).magnitude 111 112 I = ( 113 integrate.quad(pdfarg, r_min.to(xunit).magnitude, np.inf)[0] 114 * yunit 115 * xunit**4 116 ) 117 return (I * 4 / 3 * np.pi * self.rho_w / self.rho_a).to( 118 self.si.gram / self.si.kilogram 119 ) 120 121 def pdf(self, r): 122 return self.C * self.size_distribution.pdf(r)
Settings( nr: int = 64, mixing_ratios_g_kg: numpy.ndarray = array([ 1, 2, 4, 6, 8, 10]))
58 def __init__( 59 self, 60 nr: int = default_nr, 61 mixing_ratios_g_kg: np.ndarray = default_mixing_ratios_g_kg, 62 ): 63 self.si = si 64 self.nr = nr 65 self.r_min = 1 * si.micrometre 66 self.r_max = 25 * si.micrometre 67 self.rho_w = 1 * si.kilogram / si.decimetre**3 68 self.rho_a = 1 * si.kilogram / si.metre**3 69 self.mixing_ratios = mixing_ratios_g_kg * si.gram / si.kilogram 70 S = 1.00075 71 self.drdt = equilibrium_drop_growth.DrDt(ksi_1, S) 72 self.size_distribution = East_and_Marshall_1954.SizeDistribution(si) 73 74 self.C = (1 * si.gram / si.kilogram) / self.mixing_ratio( 75 self.size_distribution.pdf 76 ) 77 np.testing.assert_approx_equal(self.C.to(si.dimensionless).magnitude, 1, 4) 78 self.out_times = self.find_out_steps()
def
find_out_steps(self):
80 def find_out_steps(self): 81 out_steps = [] 82 for mr in self.mixing_ratios: 83 t_unit = self.si.second 84 85 def findroot(ti, mr=mr, t_unit=t_unit): 86 return ( 87 mr 88 - self.mixing_ratio( 89 equilibrium_drop_growth.PdfEvolver( 90 self.pdf, self.drdt, ti * t_unit 91 ) 92 ) 93 ).magnitude 94 95 t = optimize.brentq(findroot, 0, (1 * self.si.hour).to(t_unit).magnitude) 96 out_steps.append(t) 97 return out_steps
def
mixing_ratio(self, pdf):
99 def mixing_ratio(self, pdf): 100 xunit = self.si.micrometre 101 yunit = 1 / self.si.micrometre / self.si.centimetre**3 102 103 r_min = 0.1 * self.si.um 104 while not np.isfinite(pdf(r_min).magnitude): 105 r_min *= 1.01 106 107 def pdfarg(r_nounit): 108 r = r_nounit * xunit 109 result = pdf(r) * r**3 110 return result.to(yunit * xunit**3).magnitude 111 112 I = ( 113 integrate.quad(pdfarg, r_min.to(xunit).magnitude, np.inf)[0] 114 * yunit 115 * xunit**4 116 ) 117 return (I * 4 / 3 * np.pi * self.rho_w / self.rho_a).to( 118 self.si.gram / self.si.kilogram 119 )