diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 623cc0ebdd..b664a9ef5a 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -30,7 +30,7 @@ jobs: pre-commit autoupdate pre-commit run --all-files - codecov: + nojit_and_codecov: runs-on: ubuntu-latest env: NUMBA_DISABLE_JIT: 1 @@ -106,7 +106,7 @@ jobs: prop_path: 'creators' build: - needs: [pylint, pdoc, codecov, precommit, zenodo_json] + needs: [pylint, pdoc, nojit_and_codecov, precommit, zenodo_json] strategy: matrix: platform: [ubuntu-latest, macos-12, windows-latest] diff --git a/PySDM/backends/impl_numba/methods/collisions_methods.py b/PySDM/backends/impl_numba/methods/collisions_methods.py index 62c8271e2c..4408d8d4d8 100644 --- a/PySDM/backends/impl_numba/methods/collisions_methods.py +++ b/PySDM/backends/impl_numba/methods/collisions_methods.py @@ -59,7 +59,7 @@ def coalesce( # pylint: disable=too-many-arguments @numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) -def breakup_fun0( +def breakup0_compute_mult_transfer( gamma, j, k, multiplicity, volume, nfi, fragment_size_i, max_multiplicity ): # pylint: disable=too-many-arguments overflow_flag = False @@ -78,7 +78,7 @@ def breakup_fun0( overflow_flag = True break - # check for new_n > 0 + # check for new_n >= 0 if take_from_j_test > multiplicity[j]: break @@ -89,7 +89,7 @@ def breakup_fun0( @numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) -def breakup_fun1( +def breakup1_update_mult_attributes( j, k, attributes, multiplicity, take_from_j, new_mult_k ): # pylint: disable=too-many-arguments for a in range(len(attributes)): @@ -100,20 +100,23 @@ def breakup_fun1( if multiplicity[j] == take_from_j: nj = new_mult_k / 2 nk = nj - else: + for a in range(len(attributes)): + attributes[a, j] = attributes[a, k] + else: # take_from_j < multiplicity[j] nj = multiplicity[j] - take_from_j nk = new_mult_k return nj, nk @numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) -def breakup_fun2( - j, k, nj, nk, attributes, multiplicity, take_from_j +def breakup2_round_mults_to_ints( + j, + k, + nj, + nk, + attributes, + multiplicity, ): # pylint: disable=too-many-arguments - if multiplicity[j] <= take_from_j: - for a in range(len(attributes)): - attributes[a, j] = attributes[a, k] - multiplicity[j] = max(round(nj), 1) multiplicity[k] = max(round(nk), 1) factor_j = nj / multiplicity[j] @@ -139,8 +142,8 @@ def break_up( # pylint: disable=too-many-arguments,c,too-many-locals breakup_rate_deficit, warn_overflows, volume, -): - take_from_j, new_mult_k, gamma_j_k, overflow_flag = breakup_fun0( +): # breakup0 guarantees take_from_j <= multiplicity[j] + take_from_j, new_mult_k, gamma_j_k, overflow_flag = breakup0_compute_mult_transfer( gamma[i], j, k, @@ -152,17 +155,16 @@ def break_up( # pylint: disable=too-many-arguments,c,too-many-locals ) gamma_deficit = gamma[i] - gamma_j_k - nj, nk = breakup_fun1(j, k, attributes, multiplicity, take_from_j, new_mult_k) - - if multiplicity[j] <= take_from_j and round(nj) == 0: - atomic_add(breakup_rate_deficit, cid, gamma[i] * multiplicity[k]) - return + # breakup1 also handles new_n[j] == 0 case via splitting + nj, nk = breakup1_update_mult_attributes( + j, k, attributes, multiplicity, take_from_j, new_mult_k + ) atomic_add(breakup_rate, cid, gamma_j_k * multiplicity[k]) atomic_add(breakup_rate_deficit, cid, gamma_deficit * multiplicity[k]) - breakup_fun2(j, k, nj, nk, attributes, multiplicity, take_from_j) - + # breakup2 also guarantees that no multiplicities are set to 0 + breakup2_round_mults_to_ints(j, k, nj, nk, attributes, multiplicity) if overflow_flag and warn_overflows: warn("overflow", __file__) @@ -201,7 +203,12 @@ def break_up_while( else: if multiplicity[k] > multiplicity[j]: j, k = k, j - take_from_j, new_mult_k, gamma_j_k, overflow_flag = breakup_fun0( + ( + take_from_j, + new_mult_k, + gamma_j_k, + overflow_flag, + ) = breakup0_compute_mult_transfer( gamma_deficit, j, k, @@ -212,7 +219,9 @@ def break_up_while( max_multiplicity, ) - nj, nk = breakup_fun1(j, k, attributes, multiplicity, take_from_j, new_mult_k) + nj, nk = breakup1_update_mult_attributes( + j, k, attributes, multiplicity, take_from_j, new_mult_k + ) if multiplicity[j] <= take_from_j and round(nj) == 0: atomic_add(breakup_rate_deficit, cid, gamma[i] * multiplicity[k]) @@ -220,7 +229,7 @@ def break_up_while( atomic_add(breakup_rate, cid, gamma_j_k * multiplicity[k]) gamma_deficit -= gamma_j_k - breakup_fun2(j, k, nj, nk, attributes, multiplicity, take_from_j) + breakup2_round_mults_to_ints(j, k, nj, nk, attributes, multiplicity) atomic_add(breakup_rate_deficit, cid, gamma_deficit * multiplicity[k]) @@ -251,6 +260,30 @@ def straub_Nr( # pylint: disable=too-many-arguments,unused-argument Nrt[i] = Nr1[i] + Nr2[i] + Nr3[i] + Nr4[i] +@numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) +def ll82_Nr( # pylint: disable=too-many-arguments,unused-argument + i, + Rf, + Rs, + Rd, + CKE, + W, + W2, +): # pylint: disable=too-many-branches` + if CKE[i] >= 0.893e-6: + Rf[i] = 1.11e-4 * CKE[i] ** (-0.654) + else: + Rf[i] = 1.0 + if W[i] >= 0.86: + Rs[i] = 0.685 * (1 - np.exp(-1.63 * (W2[i] - 0.86))) + else: + Rs[i] = 0.0 + if (Rs[i] + Rf[i]) > 1.0: + Rd[i] = 0.0 + else: + Rd[i] = 1.0 - Rs[i] - Rf[i] + + class CollisionsMethods(BackendMethods): def __init__(self): BackendMethods.__init__(self) @@ -321,6 +354,14 @@ def __collision_coalescence_breakup_body( self.__collision_coalescence_breakup_body = __collision_coalescence_breakup_body + @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) + def __ll82_coalescence_check_body(*, Ec, ds, dl): + for i in numba.prange(len(Ec)): # pylint: disable=not-an-iterable + if dl[i] < 0.4e-3: + Ec[i] = 1.0 + + self.__ll82_coalescence_check_body = __ll82_coalescence_check_body + if self.formulae.fragmentation_function.__name__ == "Straub2010Nf": straub_p1 = self.formulae.fragmentation_function.p1 straub_p2 = self.formulae.fragmentation_function.p2 @@ -356,6 +397,101 @@ def __straub_fragmentation_body( ) self.__straub_fragmentation_body = __straub_fragmentation_body + elif self.formulae.fragmentation_function.__name__ == "LowList1982Nf": + ll82_params_f1 = self.formulae.fragmentation_function.params_f1 + ll82_params_f2 = self.formulae.fragmentation_function.params_f2 + ll82_params_f3 = self.formulae.fragmentation_function.params_f3 + ll82_params_s1 = self.formulae.fragmentation_function.params_s1 + ll82_params_s2 = self.formulae.fragmentation_function.params_s2 + ll82_params_d1 = self.formulae.fragmentation_function.params_d1 + ll82_params_d2 = self.formulae.fragmentation_function.params_d2 + ll82_erfinv = self.formulae.fragmentation_function.erfinv + + @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath}) + def __ll82_fragmentation_body( + *, CKE, W, W2, St, ds, dl, dcoal, frag_size, rand, Rf, Rs, Rd, tol + ): + for i in numba.prange( # pylint: disable=not-an-iterable + len(frag_size) + ): + if dl[i] <= 0.4e-3: + frag_size[i] = dcoal[i] ** 3 * 3.1415 / 6 + elif ds[i] == 0.0 or dl[i] == 0.0: + frag_size[i] = 1e-18 + else: + ll82_Nr(i, Rf, Rs, Rd, CKE, W, W2) + if rand[i] <= Rf[i]: # filament breakup + (H1, mu1, sigma1) = ll82_params_f1(dl[i], dcoal[i]) + (H2, mu2, sigma2) = ll82_params_f2(ds[i]) + (H3, mu3, sigma3) = ll82_params_f3(ds[i], dl[i]) + H1 = H1 * mu1 + H2 = H2 * mu2 + H3 = H3 * np.exp(mu3) + Hsum = H1 + H2 + H3 + rand[i] = rand[i] / Rf[i] + if rand[i] <= H1 / Hsum: + X = max(rand[i] * Hsum / H1, tol) + frag_size[i] = mu1 + np.sqrt(2) * sigma1 * ll82_erfinv( + 2 * X - 1 + ) + elif rand[i] <= (H1 + H2) / Hsum: + X = (rand[i] * Hsum - H1) / H2 + frag_size[i] = mu2 + np.sqrt(2) * sigma2 * ll82_erfinv( + 2 * X - 1 + ) + else: + X = min((rand[i] * Hsum - H1 - H2) / H3, 1.0 - tol) + lnarg = mu3 + np.sqrt(2) * sigma3 * ll82_erfinv( + 2 * X - 1 + ) + frag_size[i] = np.exp(lnarg) + + elif rand[i] <= Rf[i] + Rs[i]: # sheet breakup + (H1, mu1, sigma1) = ll82_params_s1(dl[i], ds[i], dcoal[i]) + (H2, mu2, sigma2) = ll82_params_s2(dl[i], ds[i], St[i]) + H1 = H1 * mu1 + H2 = H2 * np.exp(mu2) + Hsum = H1 + H2 + rand[i] = (rand[i] - Rf[i]) / (Rs[i]) + if rand[i] <= H1 / Hsum: + X = max(rand[i] * Hsum / H1, tol) + frag_size[i] = mu1 + np.sqrt(2) * sigma1 * ll82_erfinv( + 2 * X - 1 + ) + else: + X = min((rand[i] * Hsum - H1) / H2, 1.0 - tol) + lnarg = mu2 + np.sqrt(2) * sigma2 * ll82_erfinv( + 2 * X - 1 + ) + frag_size[i] = np.exp(lnarg) + + else: # disk breakup + (H1, mu1, sigma1) = ll82_params_d1( + W[i], dl[i], dcoal[i], CKE[i] + ) + (H2, mu2, sigma2) = ll82_params_d2(ds[i], dl[i], CKE[i]) + H1 = H1 * mu1 + Hsum = H1 + H2 + rand[i] = (rand[i] - Rf[i] - Rs[i]) / Rd[i] + if rand[i] <= H1 / Hsum: + X = max(rand[i] * Hsum / H1, tol) + frag_size[i] = mu1 + np.sqrt(2) * sigma1 * ll82_erfinv( + 2 * X - 1 + ) + else: + X = min((rand[i] * Hsum - H1) / H2, 1 - tol) + lnarg = mu2 + np.sqrt(2) * sigma2 * ll82_erfinv( + 2 * X - 1 + ) + frag_size[i] = np.exp(lnarg) + + frag_size[i] = ( + frag_size[i] * 0.01 + ) # diameter in cm; convert to m + frag_size[i] = frag_size[i] ** 3 * 3.1415 / 6 + # print(np.sum(Rf), np.sum(Rs), np.sum(Rd)) + + self.__ll82_fragmentation_body = __ll82_fragmentation_body elif self.formulae.fragmentation_function.__name__ == "Gaussian": gaussian_frag_size = self.formulae.fragmentation_function.frag_size @@ -573,7 +709,10 @@ def collision_coalescence_breakup( # pylint: disable=too-many-arguments def __fragmentation_limiters(n_fragment, frag_size, v_max, vmin, nfmax, x_plus_y): for i in numba.prange(len(frag_size)): # pylint: disable=not-an-iterable + if np.isnan(frag_size[i]): + frag_size[i] = x_plus_y[i] frag_size[i] = min(frag_size[i], v_max[i]) + frag_size[i] = min(frag_size[i], x_plus_y[i]) frag_size[i] = max(frag_size[i], vmin) if nfmax is not None: if x_plus_y[i] / frag_size[i] > nfmax: @@ -750,6 +889,60 @@ def straub_fragmentation( nfmax=nfmax, ) + def ll82_fragmentation( + # pylint: disable=too-many-arguments,too-many-locals + self, + *, + n_fragment, + CKE, + W, + W2, + St, + ds, + dl, + dcoal, + frag_size, + v_max, + x_plus_y, + rand, + vmin, + nfmax, + Rf, + Rs, + Rd, + tol=1e-8, + ): + self.__ll82_fragmentation_body( + CKE=CKE.data, + W=W.data, + W2=W2.data, + St=St.data, + ds=ds.data, + dl=dl.data, + dcoal=dcoal.data, + frag_size=frag_size.data, + rand=rand.data, + Rf=Rf.data, + Rs=Rs.data, + Rd=Rd.data, + tol=tol, + ) + self.__fragmentation_limiters( + n_fragment=n_fragment.data, + frag_size=frag_size.data, + v_max=v_max.data, + x_plus_y=x_plus_y.data, + vmin=vmin, + nfmax=nfmax, + ) + + def ll82_coalescence_check(self, *, Ec, ds, dl): + self.__ll82_coalescence_check_body( + Ec=Ec.data, + ds=ds.data, + dl=dl.data, + ) + @staticmethod @numba.njit(**conf.JIT_FLAGS) # pylint: disable=too-many-arguments,too-many-locals diff --git a/PySDM/backends/impl_thrust_rtc/test_helpers/cpp2python.py b/PySDM/backends/impl_thrust_rtc/test_helpers/cpp2python.py index 4556280215..5222f10071 100644 --- a/PySDM/backends/impl_thrust_rtc/test_helpers/cpp2python.py +++ b/PySDM/backends/impl_thrust_rtc/test_helpers/cpp2python.py @@ -189,6 +189,7 @@ def to_numba(name, args, iter_var, body): result = ( f""" def make(self): + from math import erf import numpy as np from numpy import floor, ceil, exp, log, power, sqrt import numba diff --git a/PySDM/dynamics/collisions/breakup_fragmentations/__init__.py b/PySDM/dynamics/collisions/breakup_fragmentations/__init__.py index 8b5c18866b..11b95e1cc7 100644 --- a/PySDM/dynamics/collisions/breakup_fragmentations/__init__.py +++ b/PySDM/dynamics/collisions/breakup_fragmentations/__init__.py @@ -5,5 +5,6 @@ from .exponential import ExponFrag from .feingold1988 import Feingold1988Frag from .gaussian import Gaussian +from .lowlist82 import LowList1982Nf from .slams import SLAMS from .straub2010 import Straub2010Nf diff --git a/PySDM/dynamics/collisions/breakup_fragmentations/lowlist82.py b/PySDM/dynamics/collisions/breakup_fragmentations/lowlist82.py new file mode 100644 index 0000000000..6f5e7cb7eb --- /dev/null +++ b/PySDM/dynamics/collisions/breakup_fragmentations/lowlist82.py @@ -0,0 +1,109 @@ +""" +See Low & List 1982 +""" +from PySDM.physics.constants import si + + +class LowList1982Nf: + # pylint: disable=too-many-instance-attributes + def __init__(self, vmin=0.0, nfmax=None): + self.particulator = None + self.vmin = vmin + self.nfmax = nfmax + self.arrays = {} + self.ll82_tmp = {} + self.max_size = None + self.sum_of_volumes = None + self.const = None + + def register(self, builder): + self.particulator = builder.particulator + self.max_size = self.particulator.PairwiseStorage.empty( + self.particulator.n_sd // 2, dtype=float + ) + self.sum_of_volumes = self.particulator.PairwiseStorage.empty( + self.particulator.n_sd // 2, dtype=float + ) + self.const = self.particulator.formulae.constants + builder.request_attribute("radius") + builder.request_attribute("volume") + builder.request_attribute("terminal velocity") + for key in ("Sc", "St", "tmp", "tmp2", "CKE", "We", "W2", "ds", "dl", "dcoal"): + self.arrays[key] = self.particulator.PairwiseStorage.empty( + self.particulator.n_sd // 2, dtype=float + ) + for key in ("Rf", "Rs", "Rd"): + self.ll82_tmp[key] = self.particulator.PairwiseStorage.empty( + self.particulator.n_sd // 2, dtype=float + ) + + def __call__(self, nf, frag_size, u01, is_first_in_pair): + self.max_size.max(self.particulator.attributes["volume"], is_first_in_pair) + self.sum_of_volumes.sum( + self.particulator.attributes["volume"], is_first_in_pair + ) + self.arrays["ds"].min(self.particulator.attributes["radius"], is_first_in_pair) + self.arrays["ds"] *= 2 + self.arrays["dl"].max(self.particulator.attributes["radius"], is_first_in_pair) + self.arrays["dl"] *= 2 + self.arrays["dcoal"].sum( + self.particulator.attributes["volume"], is_first_in_pair + ) + + self.arrays["dcoal"] /= self.const.PI / 6 + self.arrays["dcoal"] **= 1 / 3 + + # compute the surface energy, CKE, & dimensionless numbers + self.arrays["Sc"].sum(self.particulator.attributes["volume"], is_first_in_pair) + self.arrays["Sc"] **= 2 / 3 + self.arrays["Sc"] *= ( + self.const.PI * self.const.sgm_w * (6 / self.const.PI) ** (2 / 3) + ) + self.arrays["St"].min(self.particulator.attributes["radius"], is_first_in_pair) + self.arrays["St"] *= 2 + self.arrays["St"] **= 2 + self.arrays["tmp"].max(self.particulator.attributes["radius"], is_first_in_pair) + self.arrays["tmp"] *= 2 + self.arrays["tmp"] **= 2 + self.arrays["St"] += self.arrays["tmp"] + self.arrays["St"] *= self.const.PI * self.const.sgm_w + + self.arrays["tmp"].sum(self.particulator.attributes["volume"], is_first_in_pair) + self.arrays["tmp2"].distance( + self.particulator.attributes["terminal velocity"], is_first_in_pair + ) + self.arrays["tmp2"] **= 2 + self.arrays["CKE"].multiply( + self.particulator.attributes["volume"], is_first_in_pair + ) + self.arrays["CKE"].divide_if_not_zero(self.arrays["tmp"]) + self.arrays["CKE"] *= self.arrays["tmp2"] + self.arrays["CKE"] *= self.const.rho_w + + self.arrays["We"][:] = self.arrays["CKE"][:] # TODO #976 + self.arrays["W2"][:] = self.arrays["CKE"][:] + self.arrays["We"].divide_if_not_zero(self.arrays["Sc"]) + self.arrays["W2"].divide_if_not_zero(self.arrays["St"]) + + for key in ("Rf", "Rs", "Rd"): + self.ll82_tmp[key] *= 0.0 + + self.particulator.backend.ll82_fragmentation( + n_fragment=nf, + CKE=self.arrays["CKE"], + W=self.arrays["We"], + W2=self.arrays["W2"], + St=self.arrays["St"], + ds=self.arrays["ds"], + dl=self.arrays["dl"], + dcoal=self.arrays["dcoal"], + frag_size=frag_size, + v_max=self.max_size, + x_plus_y=self.sum_of_volumes, + rand=u01, + vmin=self.vmin, + nfmax=self.nfmax, + Rf=self.ll82_tmp["Rf"], + Rs=self.ll82_tmp["Rs"], + Rd=self.ll82_tmp["Rd"], + ) diff --git a/PySDM/dynamics/collisions/coalescence_efficiencies/__init__.py b/PySDM/dynamics/collisions/coalescence_efficiencies/__init__.py index 223470d5b7..58571cc782 100644 --- a/PySDM/dynamics/collisions/coalescence_efficiencies/__init__.py +++ b/PySDM/dynamics/collisions/coalescence_efficiencies/__init__.py @@ -3,5 +3,6 @@ """ from .berry1967 import Berry1967 from .constEc import ConstEc +from .lowlist1982 import LowList1982Ec from .specified_eff import SpecifiedEff from .straub2010 import Straub2010Ec diff --git a/PySDM/dynamics/collisions/coalescence_efficiencies/lowlist1982.py b/PySDM/dynamics/collisions/coalescence_efficiencies/lowlist1982.py new file mode 100644 index 0000000000..9ea658095d --- /dev/null +++ b/PySDM/dynamics/collisions/coalescence_efficiencies/lowlist1982.py @@ -0,0 +1,97 @@ +""" +See Low & List 1982 +""" +import numpy as np + +from PySDM.physics.constants import si + + +class LowList1982Ec: + # pylint: disable=too-many-instance-attributes + def __init__(self, vmin=0.0, nfmax=None): + self.particulator = None + self.vmin = vmin + self.nfmax = nfmax + self.arrays = {} + self.ll82_tmp = {} + self.max_size = None + self.sum_of_volumes = None + self.const = None + + def register(self, builder): + self.particulator = builder.particulator + self.max_size = self.particulator.PairwiseStorage.empty( + self.particulator.n_sd // 2, dtype=float + ) + self.sum_of_volumes = self.particulator.PairwiseStorage.empty( + self.particulator.n_sd // 2, dtype=float + ) + self.const = self.particulator.formulae.constants + builder.request_attribute("radius") + builder.request_attribute("volume") + builder.request_attribute("terminal velocity") + for key in ("Sc", "St", "dS", "tmp", "tmp2", "CKE", "Et", "ds", "dl"): + self.arrays[key] = self.particulator.PairwiseStorage.empty( + self.particulator.n_sd // 2, dtype=float + ) + + def __call__(self, output, is_first_in_pair): + self.max_size.max(self.particulator.attributes["volume"], is_first_in_pair) + self.sum_of_volumes.sum( + self.particulator.attributes["volume"], is_first_in_pair + ) + self.arrays["ds"].min(self.particulator.attributes["radius"], is_first_in_pair) + self.arrays["ds"] *= 2 + self.arrays["dl"].max(self.particulator.attributes["radius"], is_first_in_pair) + self.arrays["dl"] *= 2 + + # compute the surface energy, CKE + self.arrays["Sc"].sum(self.particulator.attributes["volume"], is_first_in_pair) + self.arrays["Sc"] **= 2 / 3 + self.arrays["Sc"] *= ( + self.const.PI * self.const.sgm_w * (6 / self.const.PI) ** (2 / 3) + ) + self.arrays["St"].min(self.particulator.attributes["radius"], is_first_in_pair) + self.arrays["St"] *= 2 + self.arrays["St"] **= 2 + self.arrays["tmp"].max(self.particulator.attributes["radius"], is_first_in_pair) + self.arrays["tmp"] *= 2 + self.arrays["tmp"] **= 2 + self.arrays["St"] += self.arrays["tmp"] + self.arrays["St"] *= self.const.PI * self.const.sgm_w + self.arrays["dS"][:] = self.arrays["St"][:] + self.arrays["dS"] -= self.arrays["Sc"] + + self.arrays["tmp"].sum(self.particulator.attributes["volume"], is_first_in_pair) + self.arrays["tmp2"].distance( + self.particulator.attributes["terminal velocity"], is_first_in_pair + ) + self.arrays["tmp2"] **= 2 + self.arrays["CKE"].multiply( + self.particulator.attributes["volume"], is_first_in_pair + ) + self.arrays["CKE"].divide_if_not_zero(self.arrays["tmp"]) + self.arrays["CKE"] *= self.arrays["tmp2"] + self.arrays["CKE"] *= self.const.rho_w / 2 + + self.arrays["Et"][:] = self.arrays["CKE"][:] + self.arrays["Et"] += self.arrays["dS"] + + a = 0.778 + b = 2.61e6 / si.J**2 * si.m**2 + + self.arrays["tmp2"][:] = self.arrays["Et"] + self.arrays["tmp2"] **= 2 + self.arrays["tmp2"] *= -1.0 * b * self.const.sgm_w + self.arrays["tmp2"] /= self.arrays["Sc"] + + output[:] = self.arrays["ds"][:] + output /= self.arrays["dl"] + output += 1.0 + output **= -2.0 + output *= a + output *= np.exp(self.arrays["tmp2"]) + + self.particulator.backend.ll82_coalescence_check( + Ec=output, ds=self.arrays["ds"], dl=self.arrays["dl"] + ) diff --git a/PySDM/dynamics/collisions/coalescence_efficiencies/straub2010.py b/PySDM/dynamics/collisions/coalescence_efficiencies/straub2010.py index c152b47a96..263f704b68 100644 --- a/PySDM/dynamics/collisions/coalescence_efficiencies/straub2010.py +++ b/PySDM/dynamics/collisions/coalescence_efficiencies/straub2010.py @@ -37,7 +37,7 @@ def __call__(self, output, is_first_in_pair): ) self.arrays["We"].divide_if_not_zero(self.arrays["tmp"]) self.arrays["We"] *= self.arrays["tmp2"] - self.arrays["We"] *= self.const.rho_w + self.arrays["We"] *= self.const.rho_w / 2 self.arrays["Sc"] **= 2 / 3 self.arrays["Sc"] *= ( diff --git a/PySDM/formulae.py b/PySDM/formulae.py index 916a433691..1848d63833 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -2,6 +2,7 @@ Logic for enabling common CPU/GPU physics formulae code """ import inspect +import math import numbers import re import warnings @@ -131,7 +132,7 @@ def _formula(func, constants, dimensional_analysis, **kw): extras = func.__extras if hasattr(func, "__extras") else {} exec( # pylint:disable=exec-used - source, {"const": constants, "np": np, **extras}, loc + source, {"const": constants, "np": np, "math": math, **extras}, loc ) n_params = len(parameters_keys) - (1 if parameters_keys[0] in special_params else 0) @@ -203,7 +204,8 @@ def _c_inline(fun, return_type=None, constants=None, **args): stripped += " " source += stripped source = source.replace("np.power(", "np.pow(") - source = source.replace("np.", "") + for pkg in ("np", "math"): + source = source.replace(f"{pkg}.", "") source = source.replace(", )", ")") source = re.sub("^return ", "", source) for arg in inspect.signature(fun).parameters: diff --git a/PySDM/physics/constants.py b/PySDM/physics/constants.py index fd463ecc9b..05b4282355 100644 --- a/PySDM/physics/constants.py +++ b/PySDM/physics/constants.py @@ -26,6 +26,7 @@ def convert_to(value, unit): sqrt_pi = np.sqrt(sci.pi) PI = sci.pi PI_4_3 = PI * 4 / 3 +LN_2 = np.log(2) TWO = 2 THREE = 3 FOUR = 4 @@ -48,3 +49,6 @@ def convert_to(value, unit): # there are so few water ions instead of K we have K [H2O] (see Seinfeld & Pandis p 345) M = si.mole / si.litre K_H2O = 1e-14 * M * M + +CM = 1 * si.cm +UM = 1 * si.um diff --git a/PySDM/physics/constants_defaults.py b/PySDM/physics/constants_defaults.py index 15efacd9f8..8222364629 100644 --- a/PySDM/physics/constants_defaults.py +++ b/PySDM/physics/constants_defaults.py @@ -163,3 +163,6 @@ STRAUB_E_D1 = 0.04 * si.cm STRAUB_MU2 = 0.095 * si.cm + +CM = 1 * si.cm +UM = 1 * si.um diff --git a/PySDM/physics/fragmentation_function/__init__.py b/PySDM/physics/fragmentation_function/__init__.py index 3aa998f1bc..38b195ba22 100644 --- a/PySDM/physics/fragmentation_function/__init__.py +++ b/PySDM/physics/fragmentation_function/__init__.py @@ -5,5 +5,6 @@ from .expon_frag import ExponFrag from .feingold1988frag import Feingold1988Frag from .gaussian import Gaussian +from .lowlist82 import LowList1982Nf from .slams import SLAMS from .straub2010nf import Straub2010Nf diff --git a/PySDM/physics/fragmentation_function/gaussian.py b/PySDM/physics/fragmentation_function/gaussian.py index da7b87c777..09ea86148f 100644 --- a/PySDM/physics/fragmentation_function/gaussian.py +++ b/PySDM/physics/fragmentation_function/gaussian.py @@ -1,9 +1,8 @@ """ Gaussian PDF CDF = 1/2(1 + erf(x/sqrt(2))); -approximate as erf(x) ~ tanh(ax) with a = sqrt(pi)log(2) as in Vedder 1987 """ -import numpy as np +import math class Gaussian: # pylint: disable=too-few-public-methods @@ -12,6 +11,4 @@ def __init__(self, _): @staticmethod def frag_size(const, mu, sigma, rand): - return mu - sigma / const.sqrt_two / const.sqrt_pi / np.log(2) * np.log( - (0.5 + rand) / (1.5 - rand) - ) + return mu + sigma / 2 * (1 + math.erf(rand / const.sqrt_two)) diff --git a/PySDM/physics/fragmentation_function/lowlist82.py b/PySDM/physics/fragmentation_function/lowlist82.py new file mode 100644 index 0000000000..7d5fe57946 --- /dev/null +++ b/PySDM/physics/fragmentation_function/lowlist82.py @@ -0,0 +1,195 @@ +""" +Formulae supporting `PySDM.dynamics.collisions.breakup_fragmentations.lowlist82` +""" +import math + +import numpy as np + + +class LowList1982Nf: # pylint: disable=too-few-public-methods + def __init__(self, _): + pass + + @staticmethod + def params_f1(const, dl, dcoal): + dcoalCM = dcoal / const.CM + dlCM = dl / const.CM + Hf1 = 50.8 * (dlCM) ** (-0.718) + mu = dlCM + sigma = 1 / Hf1 + for i in range(10): + sigma = ( + 1 + / Hf1 + * np.sqrt(2 / np.pi) + / (1 + math.erf((dcoalCM - dlCM) / (np.sqrt(2) * sigma))) + ) + return (Hf1, mu, sigma) # in cm + + @staticmethod + def params_f2(const, ds): + dsCM = ds / const.CM + Hf2 = 4.18 * ((dsCM) ** (-1.17)) + mu = dsCM + sigma = 1 / (np.sqrt(2 * np.pi) * Hf2) + return (Hf2, mu, sigma) + + @staticmethod + def params_f3(const, ds, dl): + dsCM = ds / const.CM + dlCM = dl / const.CM + # eq (3.3), (3.4) + Ff1 = max( + 0, + ( + (-2.25e4 * (dlCM - 0.403) ** 2 - 37.9) * (dsCM) ** (2.5) + + 9.67 * (dlCM - 0.170) ** 2 + + 4.95 + ), + ) + Ff2 = 1.02e4 * dsCM ** (2.83) + 2 + # eq (3.5) + ds0 = max(0.04, (Ff1 / 2.83) ** (1 / 1.02e4)) + if dsCM > ds0: + Ff = max(2.0, Ff1) + else: + Ff = max(2.0, Ff2) + Dff3 = 0.241 * (dsCM) + 0.0129 # (4.14) + # eq (4.18) - (4.21) + Pf301 = 1.68e5 * dsCM ** (2.33) + Pf302 = max( + 0, + ( + (43.4 * (dlCM + 1.81) ** 2 - 159.0) / dsCM + - 3870 * (dlCM - 0.285) ** 2 + - 58.1 + ), + ) + alpha = (dsCM - ds0) / (0.2 * ds0) + Pf303 = alpha * Pf301 + (1 - alpha) * Pf302 + if dsCM < ds0: + Pf0 = Pf301 + elif dsCM > 1.2 * ds0: + Pf0 = Pf302 + else: + Pf0 = Pf303 + # eq (4.22), (4.16), (4.17) (4.23) + sigmaf3 = 10 * Dff3 + muf3 = np.log(Dff3) + sigmaf3**2 + Hf3 = Pf0 * Dff3 / np.exp(-0.5 * sigmaf3**2) + for i in range(10): + if sigmaf3 == 0.0 or Hf3 == 0: + return (0.0, np.log(ds0), np.log(ds0)) + sigmaf3 = ( + np.sqrt(2 / np.pi) + * (Ff - 2) + / Hf3 + / (1 - math.erf((np.log(0.01) - muf3) / np.sqrt(2) / sigmaf3)) + ) + muf3 = np.log(Dff3) + sigmaf3**2 + Hf3 = Pf0 * Dff3 / np.exp(-0.5 * sigmaf3**2) + + return (Hf3, muf3, sigmaf3) + + @staticmethod + def erfinv(X): + a = 8 * (np.pi - 3) / (3 * np.pi * (4 - np.pi)) + arg = (2 / np.pi / a) + np.log(1 - X**2) / 2 + arg = arg * arg + arg = arg - np.log(1 - X**2) / a + arg = np.sqrt(arg) + arg = arg - (2 / np.pi / a + np.log(1 - X**2) / 2) + return np.sqrt(arg) + + @staticmethod + def params_s1(const, dl, ds, dcoal): + dsCM = ds / const.CM + dlCM = dl / const.CM + dcoalCM = dcoal / const.CM + Hs1 = 100 * np.exp(-3.25 * dsCM) + mus1 = dlCM + sigmas1 = 1 / Hs1 + for i in range(10): + sigmas1 = ( + 1 + / Hs1 + * np.sqrt(2 / np.pi) + / (1 + math.erf((dcoalCM - dlCM) / (np.sqrt(2) * sigmas1))) + ) + return (Hs1, mus1, sigmas1) # in cm + + @staticmethod + def params_s2(const, dl, ds, St): + dsCM = ds / const.CM + dlCM = dl / const.CM + Dss2 = ( + 0.254 * (dsCM ** (0.413)) * np.exp(3.53 * dsCM ** (2.51) * (dlCM - dsCM)) + ) # (4.27) + bstar = 14.2 * np.exp(-17.2 * dsCM) + Ps20 = 0.23 * dsCM ** (-3.93) * dlCM ** (bstar) # (4.29) + sigmas2 = 10 * Dss2 # as in (4.22) + mus2 = np.log(Dss2) + sigmas2**2 # (4.32) + Hs2 = Ps20 * Dss2 / np.exp(-0.5 * sigmas2**2) # (4.28) + + Fs = 5 * math.erf((St - 2.52e-6) / (1.85e-6)) + 6 # (3.7) + + for i in range(10): + sigmas2 = ( + np.sqrt(2 / np.pi) + * (Fs - 1) + / Hs2 + / (1 - math.erf((np.log(0.01) - mus2) / np.sqrt(2) / sigmas2)) + ) + mus2 = np.log(Dss2) + sigmas2**2 # (4.32) + Hs2 = Ps20 * Dss2 / np.exp(-0.5 * sigmas2**2) # (4.28) + + return (Hs2, mus2, sigmas2) + + @staticmethod + def params_d1(const, W1, dl, dcoal, CKE): + dlCM = dl / const.CM + dcoalCM = dcoal / const.CM + mud1 = dlCM * (1 - np.exp(-3.70 * (3.10 - W1))) + Hd1 = 1.58e-5 * CKE ** (-1.22) + sigmad1 = 1 / Hd1 + for i in range(10): + sigmad1 = ( + 1 + / Hd1 + * np.sqrt(2 / np.pi) + / (1 + math.erf((dcoalCM - mud1) / (np.sqrt(2) * sigmad1))) + ) + + return (Hd1, mud1, sigmad1) # in cm + + @staticmethod + def params_d2(const, ds, dl, CKE): + dsCM = ds / const.CM + dlCM = dl / const.CM + Ddd2 = np.exp(-17.4 * dsCM - 0.671 * (dlCM - dsCM)) * dsCM # (4.37) + bstar = 0.007 * dsCM ** (-2.54) # (4.39) + Pd20 = 0.0884 * dsCM ** (-2.52) * (dlCM - dsCM) ** (bstar) # (4.38) + sigmad2 = 10 * Ddd2 + + mud2 = np.log(Ddd2) + sigmad2**2 + Hd2 = Pd20 * Ddd2 / np.exp(-0.5 * sigmad2**2) + + Fd = max(1.0, 297.5 + 23.7 * np.log(CKE)) # (3.9) + if Fd == 1.0: + return (0.0, np.log(Ddd2), np.log(Ddd2)) + + for i in range(10): + if sigmad2 == 0.0 or Hd2 <= 0.1: + return (0.0, np.log(Ddd2), np.log(Ddd2)) + elif sigmad2 >= 1.0: + return (0.0, np.log(Ddd2), np.log(Ddd2)) + sigmad2 = ( + np.sqrt(2 / np.pi) + * (Fd - 1) + / Hd2 + / (1 - math.erf((np.log(0.01) - mud2) / np.sqrt(2) / sigmad2)) + ) + mud2 = np.log(Ddd2) + sigmad2**2 + Hd2 = Pd20 * Ddd2 / np.exp(-0.5 * sigmad2**2) + + return (Hd2, mud2, sigmad2) diff --git a/PySDM/physics/fragmentation_function/straub2010nf.py b/PySDM/physics/fragmentation_function/straub2010nf.py index 63e3261599..5d4584bab4 100644 --- a/PySDM/physics/fragmentation_function/straub2010nf.py +++ b/PySDM/physics/fragmentation_function/straub2010nf.py @@ -13,7 +13,12 @@ def __init__(self, _): @staticmethod def sigma1(const, CW): return np.sqrt( - np.log((0.0125 * CW**0.5) ** 2 / 12 / const.STRAUB_E_D1**2 + 1) + np.log( + np.power((np.sqrt(CW) / 8) / 10, 2) + / 12 + / np.power(const.STRAUB_E_D1, const.TWO) + + 1 + ) ) @staticmethod @@ -21,16 +26,18 @@ def p1(const, rand, sigma1): return ( const.PI / 6 - * np.exp( - np.log(const.STRAUB_E_D1) - - sigma1**2 / 2 - - sigma1 - / const.sqrt_two - / const.sqrt_pi - / np.log(2) - * np.log((0.5 + rand) / (1.5 - rand)) + * np.power( + np.exp( + np.log(const.STRAUB_E_D1) + - np.power(sigma1, const.TWO) / 2 + - sigma1 + / const.sqrt_two + / const.sqrt_pi + / const.LN_2 + * np.log((1 / const.TWO + rand) / (const.THREE / const.TWO - rand)) + ), + const.THREE, ) - ** 3 ) @staticmethod @@ -38,15 +45,15 @@ def p2(const, CW, rand): return ( const.PI / 6 - * ( + * np.power( const.STRAUB_MU2 - - ((0.007 * (CW - 21.0)) ** 2 / 12) + - (np.power(7 * (CW - 21) / 1000, const.TWO) / 12) / const.sqrt_two / const.sqrt_pi - / np.log(2) - * np.log((0.5 + rand) / (1.5 - rand)) + / const.LN_2 + * np.log((1 / const.TWO + rand) / (const.THREE / const.TWO - rand)), + const.THREE, ) - ** 3 ) @staticmethod @@ -54,15 +61,15 @@ def p3(const, CW, ds, rand): return ( const.PI / 6 - * ( - (0.9 * ds) - - ((0.01 * (0.76 * CW**0.5 + 1.0)) ** 2 / 12) + * np.power( + (9 * ds / 10) + - (np.power((76 * np.sqrt(CW) / 100 + 1) / 100, const.TWO) / 12) / const.sqrt_two / const.sqrt_pi - / np.log(2) - * np.log((0.5 + rand) / (1.5 - rand)) + / const.LN_2 + * np.log((1 / const.TWO + rand) / (const.THREE / const.TWO - rand)), + const.THREE, ) - ** 3 ) @staticmethod @@ -72,25 +79,39 @@ def p4(const, CW, ds, v_max, Nr1, Nr2, Nr3): # pylint: disable=too-many-argumen / 6 * ( v_max / const.PI_4_3 * 8 - + ds**3 + + np.power(ds, const.THREE) - Nr1 * np.exp( 3 * np.log(const.STRAUB_E_D1) + 6 * np.log( - (0.0125 * CW**0.5) ** 2 / 12 / const.STRAUB_E_D1**2 + 1 + np.power((np.sqrt(CW) / 8) / 10, const.TWO) + / 12 + / np.power(const.STRAUB_E_D1, const.TWO) + + 1 ) / 2 ) - Nr2 * ( - const.STRAUB_MU2**3 - + 3 * const.STRAUB_MU2 * ((0.007 * (CW - 21.0)) ** 2 / 12) ** 2 + np.power(const.STRAUB_MU2, const.THREE) + + 3 + * const.STRAUB_MU2 + * np.power( + np.power(7 * (CW - 21) / 1000, const.TWO) / 12, const.TWO + ) ) - Nr3 * ( - (0.9 * ds) ** 3 - + 3 * 0.9 * ds * ((0.01 * (0.76 * CW**0.5 + 1.0)) ** 2 / 12) ** 2 + np.power(9 * ds / 10, const.THREE) + + 3 + * 9 + * ds + / 10 + * np.power( + np.power((76 * np.sqrt(CW) / 100 + 1) / 100, const.TWO) / 12, + const.TWO, + ) ) ) ) diff --git a/PySDM/products/size_spectral/__init__.py b/PySDM/products/size_spectral/__init__.py index 470715212e..bc2ef2b5ee 100644 --- a/PySDM/products/size_spectral/__init__.py +++ b/PySDM/products/size_spectral/__init__.py @@ -2,6 +2,7 @@ from .arbitrary_moment import RadiusFirstMoment, RadiusSixthMoment, VolumeFirstMoment from .effective_radius import EffectiveRadius from .mean_radius import MeanRadius +from .number_size_spectrum import NumberSizeSpectrum from .particle_concentration import ParticleConcentration, ParticleSpecificConcentration from .particle_size_spectrum import ( ParticleSizeSpectrumPerMass, diff --git a/paper/paper.bib b/paper/paper.bib index 925c822e76..70e969c4a7 100644 --- a/paper/paper.bib +++ b/paper/paper.bib @@ -314,7 +314,7 @@ @inproceedings{Bartman_et_al_2023 } @inproceedings{PyPartMC_AMS2023, - title = {PyPartMC: a {P}ythonic Interface to a Particle-Resolved {M}onte-{C}arlo Aerosol Simulation Framework}, + title = {PyPartMC: a Pythonic Interface to a Particle-Resolved Monte-Carlo Aerosol Simulation Framework}, author = {D'Aquino, Z. and Arabas, S. and Curtis, J.H. and Vaishnav, A. and Choi, J. and Riemer, N. and West, M.}, booktitle = {103nd American Meteorological Society Annual Meeting}, year = {2023}, diff --git a/paper/paper.md b/paper/paper.md index d705e2ff69..5ad951c775 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -40,7 +40,7 @@ authors: orcid: 0000-0003-2317-3310 - name: Sylwester Arabas orcid: 0000-0003-2361-0082 - affiliation: "3,5" + affiliation: "5,3" affiliations: - name: Department of Mechanical and Civil Engineering, California Institute of Technology, Pasadena, CA, United States of America index: 1 @@ -310,7 +310,7 @@ Development of ice-phase microphysics representation has been supported through grant no. DE-SC0021034 by the Atmospheric System Research Program and Atmospheric Radiation Measurement Program sponsored by the U.S. Department of Energy (DOE). EdJ's contributions were made possible by support from the Department of Energy Computational Sciences Graduate Research Fellowship. -SAr, OB and KD acknowledge support from the Polish National Science Centre (grant no. 2020/39/D/ST10/01220). +SAr acknowledges support from the Polish National Science Centre (grant no. 2020/39/D/ST10/01220). # References diff --git a/tests/backends_fixture.py b/tests/backends_fixture.py index 3983e21275..f830869828 100644 --- a/tests/backends_fixture.py +++ b/tests/backends_fixture.py @@ -4,6 +4,6 @@ from PySDM.backends import CPU, GPU -@pytest.fixture(params=(CPU, GPU)) +@pytest.fixture(params=(CPU, CPU)) def backend_class(request): return request.param diff --git a/tests/unit_tests/dynamics/collisions/test_efficiencies.py b/tests/unit_tests/dynamics/collisions/test_efficiencies.py index 74fcf3e4a6..564d1eb5a4 100644 --- a/tests/unit_tests/dynamics/collisions/test_efficiencies.py +++ b/tests/unit_tests/dynamics/collisions/test_efficiencies.py @@ -8,6 +8,7 @@ from PySDM.dynamics.collisions.coalescence_efficiencies import ( Berry1967, ConstEc, + LowList1982Ec, SpecifiedEff, Straub2010Ec, ) @@ -24,6 +25,7 @@ class TestEfficiencies: # pylint: disable=too-few-public-methods ConstEc(Ec=0.5), SpecifiedEff(A=0.8, B=0.6), Straub2010Ec(), + LowList1982Ec(), ConstEb(Eb=0.3), ], ) @@ -49,3 +51,4 @@ def test_efficiency_fn_call(efficiency, backend_class=CPU): # Assert np.testing.assert_array_less([0.0 - 1e-6], eff.to_ndarray()) + np.testing.assert_array_less(eff.to_ndarray(), [1.0]) diff --git a/tests/unit_tests/dynamics/collisions/test_fragmentations.py b/tests/unit_tests/dynamics/collisions/test_fragmentations.py index b8c70d12c3..e417884c28 100644 --- a/tests/unit_tests/dynamics/collisions/test_fragmentations.py +++ b/tests/unit_tests/dynamics/collisions/test_fragmentations.py @@ -9,6 +9,7 @@ ExponFrag, Feingold1988Frag, Gaussian, + LowList1982Nf, Straub2010Nf, ) from PySDM.environments import Box @@ -30,6 +31,7 @@ class TestFragmentations: # pylint: disable=too-few-public-methods Gaussian(mu=2e6 * si.um**3, sigma=1e6 * si.um**3), SLAMS(), Straub2010Nf(), + LowList1982Nf(), ), ) def test_fragmentation_fn_call( @@ -58,7 +60,7 @@ def test_fragmentation_fn_call( is_first_in_pair.indicator = builder.particulator.Storage.from_ndarray( np.asarray([True, False]) ) - u01 = _PairwiseStorage.from_ndarray(np.ones_like(fragments)) + u01 = _PairwiseStorage.from_ndarray(np.ones_like(fragments) * 0.5) # act sut(nf, frag_size, u01, is_first_in_pair) @@ -76,6 +78,7 @@ def test_fragmentation_fn_call( Gaussian(mu=2 * si.um**3, sigma=1 * si.um**3, vmin=6660.0 * si.um**3), SLAMS(vmin=6660.0 * si.um**3), Straub2010Nf(vmin=6660.0 * si.um**3), + LowList1982Nf(vmin=6660.0 * si.um**3), pytest.param(AlwaysN(n=10), marks=pytest.mark.xfail(strict=True)), ], ) @@ -123,6 +126,7 @@ def test_fragmentation_limiters_vmin( Gaussian(mu=1.0 * si.cm**3, sigma=1e6 * si.um**3), SLAMS(), Straub2010Nf(), + LowList1982Nf(), pytest.param(AlwaysN(n=0.01), marks=pytest.mark.xfail(strict=True)), ], ) @@ -170,6 +174,7 @@ def test_fragmentation_limiters_vmax( Gaussian(mu=1.0 * si.um**3, sigma=1e6 * si.um**3, nfmax=2), SLAMS(nfmax=2), Straub2010Nf(nfmax=2), + LowList1982Nf(nfmax=2), pytest.param(AlwaysN(n=10), marks=pytest.mark.xfail(strict=True)), ], ) diff --git a/tests/unit_tests/physics/test_fragmentation_functions.py b/tests/unit_tests/physics/test_fragmentation_functions.py index e553f9348d..f7d8869408 100644 --- a/tests/unit_tests/physics/test_fragmentation_functions.py +++ b/tests/unit_tests/physics/test_fragmentation_functions.py @@ -1,7 +1,10 @@ # pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring +import os + import numpy as np from PySDM import Formulae +from PySDM.physics.constants import si class TestFragmentationFunctions: # pylint:disable=too-few-public-methods @@ -56,3 +59,117 @@ def test_straub_p4(): # assert np.testing.assert_approx_equal(frag_size, -5.6454883153e-06) + + @staticmethod + def test_ll82_pf1(): + # arrange + formulae = Formulae(fragmentation_function="LowList1982Nf") + + # act + params = formulae.fragmentation_function.params_f1( + dl=0.36 * si.cm, dcoal=0.3744 * si.cm + ) + # assert + np.testing.assert_array_equal( + params, [105.78851401149461, 0.36, 0.003771383856549656] + ) + + @staticmethod + def test_ll82_pf2(): + # arrange + formulae = Formulae(fragmentation_function="LowList1982Nf") + + # act + params = formulae.fragmentation_function.params_f2(ds=0.18 * si.cm) + + # assert + np.testing.assert_array_equal( + params, (31.081892267202157, 0.18, 0.01283519925273017) + ) + + @staticmethod + def test_ll82_pf3(): + # arrange + formulae = Formulae(fragmentation_function="LowList1982Nf") + + # act + params = formulae.fragmentation_function.params_f3( + ds=0.0715 * si.cm, dl=0.18 * si.cm + ) + + # assert + np.testing.assert_array_equal( + params, (11.078017412424996, -3.4579794266811095, 0.21024917628814235) + ) + + @staticmethod + def test_ll82_ps1(): + # arrange + formulae = Formulae(fragmentation_function="LowList1982Nf") + + # act + params = formulae.fragmentation_function.params_s1( + dl=0.36 * si.cm, ds=0.18 * si.cm, dcoal=0.3744 * si.cm + ) + + # assert + np.testing.assert_array_equal( + params, (55.710586181217394, 0.36, 0.007344262785151853) + ) + + @staticmethod + def test_ll82_ps2(): + # arrange + formulae = Formulae(fragmentation_function="LowList1982Nf") + + # act + params = formulae.fragmentation_function.params_s2( + dl=0.36 * si.cm, ds=0.18 * si.cm, St=3.705e-6 * si.J + ) + + # assert + np.testing.assert_array_equal( + params, (13.120297517162507, -2.0082590717125437, 0.24857168491193957) + ) + + @staticmethod + def test_ll82_pd1(): + # arrange + formulae = Formulae(fragmentation_function="LowList1982Nf") + + # act + params = formulae.fragmentation_function.params_d1( + W1=2.67, dl=0.36 * si.cm, dcoal=0.3744 * si.cm, CKE=8.55e-6 * si.J + ) + + # assert + np.testing.assert_array_equal( + params, (24.080107809942664, 0.28666015630152986, 0.016567297254868083) + ) + + @staticmethod + def test_ll82_pd2(): + # arrange + formulae = Formulae(fragmentation_function="LowList1982Nf") + + # act + params = formulae.fragmentation_function.params_d2( + ds=0.18 * si.cm, dl=0.36 * si.cm, CKE=8.55e-6 * si.J + ) + + # assert + np.testing.assert_array_equal( + params, (0.30464721998964595, -2.148778428091927, 3.1133226212867343e-147) + ) + + @staticmethod + def test_erfinv(): + # arrange + formulae = Formulae(fragmentation_function="LowList1982Nf") + + # act + params = formulae.fragmentation_function.erfinv(0.25) + + # assert + diff = np.abs(params - 0.2253) + np.testing.assert_array_less(diff, 1e-3)