Skip to content

Commit 329b916

Browse files
slayooslayoo
authored
Merge pull request #888 from claresinger/surface_tension_test
surface tension tests
2 parents dade5e4 + 04680f4 commit 329b916

File tree

6 files changed

+269
-62
lines changed

6 files changed

+269
-62
lines changed

‎.github/workflows/tests+pypi.yml

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@ jobs:
2121
platform: [ubuntu-latest, macos-12, windows-latest]
2222
python-version: ["3.7", "3.8", "3.9"]
2323
runs-on: ${{ matrix.platform }}
24+
timeout-minutes: 90
2425
steps:
2526
- uses: actions/checkout@v2
2627

@@ -57,6 +58,7 @@ jobs:
5758

5859
dist:
5960
runs-on: ubuntu-latest
61+
timeout-minutes: 90
6062
needs: [build]
6163
steps:
6264
- uses: actions/checkout@v2

‎PySDM/physics/surface_tension/compressed_film_ruehl.py

Lines changed: 41 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -13,9 +13,9 @@
1313
# pylint: disable=too-many-arguments
1414
@numba.njit(**{**jit_flags, "parallel": False})
1515
def minfun(f_surf, Cb_iso, RUEHL_C0, RUEHL_A0, A_iso, c):
16-
return Cb_iso * (1 - f_surf) / RUEHL_C0 - np.exp(
17-
c * (RUEHL_A0**2 - (A_iso / f_surf) ** 2)
18-
)
16+
lhs = Cb_iso * (1 - f_surf) / RUEHL_C0
17+
rhs = np.exp(c * (RUEHL_A0**2 - (A_iso / f_surf) ** 2))
18+
return lhs - rhs
1919

2020

2121
within_tolerance = numba.njit(
@@ -48,37 +48,44 @@ def sigma(const, T, v_wet, v_dry, f_org):
4848
# wet radius (m)
4949
r_wet = ((3 * v_wet) / (4 * const.PI)) ** (1 / 3)
5050

51-
# C_bulk is the concentration of the organic in the bulk phase
52-
# Cb_iso = C_bulk / (1-f_surf)
53-
Cb_iso = (f_org * v_dry / const.RUEHL_nu_org) / (v_wet / const.nu_w)
54-
55-
# A is the area one molecule of organic occupies at the droplet surface
56-
# A_iso = A*f_surf (m^2)
57-
A_iso = (4 * const.PI * r_wet**2) / (
58-
f_org * v_dry * const.N_A / const.RUEHL_nu_org
59-
)
60-
61-
# solve implicitly for fraction of organic at surface
62-
c = (const.RUEHL_m_sigma * const.N_A) / (2 * const.R_str * T)
63-
64-
args = (Cb_iso, const.RUEHL_C0, const.RUEHL_A0, A_iso, c)
65-
rtol = 1e-6
66-
max_iters = 1e2
67-
bracket = (1e-20, 1)
68-
f_surf, iters = toms748_solve(
69-
minfun,
70-
args,
71-
*bracket,
72-
minfun(bracket[0], *args),
73-
minfun(bracket[1], *args),
74-
rtol,
75-
max_iters,
76-
within_tolerance
77-
)
78-
assert iters != max_iters
79-
80-
# calculate surface tension
81-
sgm = const.sgm_w - (const.RUEHL_A0 - A_iso / f_surf) * const.RUEHL_m_sigma
51+
if f_org == 0:
52+
sgm = const.sgm_w
53+
elif f_org == 1:
54+
sgm = const.RUEHL_sgm_min
55+
else:
56+
# C_bulk is the concentration of the organic in the bulk phase
57+
# Cb_iso = C_bulk / (1-f_surf)
58+
Cb_iso = (f_org * v_dry / const.RUEHL_nu_org) / (v_wet / const.nu_w)
59+
60+
# A is the area one molecule of organic occupies at the droplet surface
61+
# A_iso = A*f_surf (m^2)
62+
A_iso = (4 * const.PI * r_wet**2) / (
63+
f_org * v_dry * const.N_A / const.RUEHL_nu_org
64+
)
65+
66+
# solve implicitly for fraction of organic at surface
67+
c = (const.RUEHL_m_sigma * const.N_A) / (2 * const.R_str * T)
68+
69+
args = (Cb_iso, const.RUEHL_C0, const.RUEHL_A0, A_iso, c)
70+
rtol = 1e-6
71+
max_iters = 1e2
72+
bracket = (rtol, 1)
73+
f_surf, iters = toms748_solve(
74+
minfun,
75+
args,
76+
*bracket,
77+
minfun(bracket[0], *args),
78+
minfun(bracket[1], *args),
79+
rtol,
80+
max_iters,
81+
within_tolerance
82+
)
83+
assert iters != max_iters
84+
85+
# calculate surface tension
86+
sgm = const.sgm_w - (const.RUEHL_A0 - A_iso / f_surf) * const.RUEHL_m_sigma
87+
88+
# surface tension bounded between sgm_min and sgm_w
8289
sgm = np.minimum(np.maximum(sgm, const.RUEHL_sgm_min), const.sgm_w)
8390
return sgm
8491

‎PySDM/physics/surface_tension/szyszkowski_langmuir.py

Lines changed: 33 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -25,32 +25,38 @@ def __init__(self, const):
2525

2626
@staticmethod
2727
def sigma(const, T, v_wet, v_dry, f_org):
28-
r_wet = ((3 * v_wet) / (4 * np.pi)) ** (1 / 3) # m - wet radius
29-
30-
# C_bulk is the concentration of the organic in the bulk phase
31-
# Cb_iso = C_bulk / (1-f_surf)
32-
Cb_iso = (f_org * v_dry / const.RUEHL_nu_org) / (v_wet / const.nu_w)
33-
34-
# A is the area that one molecule of organic occupies at the droplet surface
35-
# A_iso = A*f_surf (m^2)
36-
A_iso = (4 * np.pi * r_wet**2) / (
37-
f_org * v_dry * const.N_A / const.RUEHL_nu_org
38-
)
39-
40-
# fraction of organic at surface
41-
# quadratic formula, solve equation of state analytically
42-
a = -const.RUEHL_A0 / A_iso
43-
b = (
44-
const.RUEHL_A0 / A_iso
45-
+ (const.RUEHL_A0 / A_iso) * (const.RUEHL_C0 / Cb_iso)
46-
+ 1
47-
)
48-
c = -1
49-
f_surf = (-b + np.sqrt(b**2 - 4 * a * c)) / (2 * a)
50-
51-
# calculate surface tension
52-
sgm = const.sgm_w - ((const.R_str * T) / (const.RUEHL_A0 * const.N_A)) * np.log(
53-
1 + Cb_iso * (1 - f_surf) / const.RUEHL_C0
54-
)
28+
# wet radius (m)
29+
r_wet = ((3 * v_wet) / (4 * np.pi)) ** (1 / 3)
30+
31+
if f_org == 0:
32+
sgm = const.sgm_w
33+
else:
34+
# C_bulk is the concentration of the organic in the bulk phase
35+
# Cb_iso = C_bulk / (1-f_surf)
36+
Cb_iso = (f_org * v_dry / const.RUEHL_nu_org) / (v_wet / const.nu_w)
37+
38+
# A is the area that one molecule of organic occupies at the droplet surface
39+
# A_iso = A*f_surf (m^2)
40+
A_iso = (4 * np.pi * r_wet**2) / (
41+
f_org * v_dry * const.N_A / const.RUEHL_nu_org
42+
)
43+
44+
# fraction of organic at surface
45+
# quadratic formula, solve equation of state analytically
46+
a = -const.RUEHL_A0 / A_iso
47+
b = (
48+
const.RUEHL_A0 / A_iso
49+
+ (const.RUEHL_A0 / A_iso) * (const.RUEHL_C0 / Cb_iso)
50+
+ 1
51+
)
52+
c = -1
53+
f_surf = (-b + np.sqrt(b**2 - 4 * a * c)) / (2 * a)
54+
55+
# calculate surface tension
56+
sgm = const.sgm_w - (
57+
(const.R_str * T) / (const.RUEHL_A0 * const.N_A)
58+
) * np.log(1 + Cb_iso * (1 - f_surf) / const.RUEHL_C0)
59+
60+
# surface tension bounded between sgm_min and sgm_w
5561
sgm = np.minimum(np.maximum(sgm, const.RUEHL_sgm_min), const.sgm_w)
5662
return sgm

‎tests/smoke_tests/lowe_et_al_2019/test_zero_forg.py

Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,101 @@
1+
# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring
2+
import numpy as np
3+
from matplotlib import pyplot
4+
from PySDM_examples.Lowe_et_al_2019 import Settings, Simulation
5+
from PySDM_examples.Lowe_et_al_2019.aerosol import AerosolBoreal, AerosolMarine
6+
7+
from PySDM.initialisation.sampling import spectral_sampling as spec_sampling
8+
from PySDM.physics import si
9+
10+
11+
def test_zero_forg(plot=False):
12+
nRes = 3
13+
updraft_list = np.geomspace(0.1, 10, nRes)
14+
15+
subplot_list = ["a", "b", "c", "d"]
16+
models = ("Constant", "CompressedFilmOvadnevaite")
17+
18+
Acc = {"a": 30, "b": 134, "c": 160, "d": 540}
19+
20+
consts = {
21+
"delta_min": 0.1,
22+
"MAC": 1,
23+
"HAC": 1,
24+
"c_pd": 1006 * si.joule / si.kilogram / si.kelvin,
25+
"g_std": 9.81 * si.metre / si.second**2,
26+
"BDF": False,
27+
}
28+
29+
cdnc_compare = np.zeros((len(models), len(subplot_list), len(updraft_list)))
30+
for i, w in enumerate(updraft_list):
31+
for k, subplot in enumerate(subplot_list):
32+
for m, model in enumerate(models):
33+
settings = Settings(
34+
dz=1 * si.m,
35+
n_sd_per_mode=20,
36+
model=model,
37+
aerosol={
38+
"a": AerosolMarine(Forg=0, Acc_N2=Acc["a"]),
39+
"b": AerosolMarine(Forg=0, Acc_N2=Acc["b"]),
40+
"c": AerosolBoreal(Forg=0, Acc_N2=Acc["c"]),
41+
"d": AerosolBoreal(Forg=0, Acc_N2=Acc["d"]),
42+
}[subplot],
43+
w=w * si.m / si.s,
44+
spectral_sampling=spec_sampling.ConstantMultiplicity,
45+
**consts,
46+
)
47+
simulation = Simulation(settings)
48+
output = simulation.run()
49+
cdnc_compare[m, k, i] = np.array(output["n_c_cm3"])[-1]
50+
51+
mrkr = ["o", "s", "*", "v", "^", "D", "h", "x", "+", "8", "p", "<", ">", "d", "H"]
52+
_, axes = pyplot.subplots(
53+
1,
54+
len(subplot_list),
55+
figsize=(len(subplot_list) * 4, 4),
56+
constrained_layout=True,
57+
sharex=True,
58+
sharey=False,
59+
)
60+
61+
for k, subplot in enumerate(subplot_list):
62+
if len(subplot_list) > 1:
63+
ax = axes[k]
64+
else:
65+
ax = axes
66+
67+
for m, model in enumerate(models):
68+
for i, w in enumerate(updraft_list):
69+
if i == 0:
70+
ax.scatter(
71+
w * (1 + 0.1 * m),
72+
cdnc_compare[m, k, i],
73+
color="C" + str(m),
74+
marker=mrkr[m],
75+
label=model,
76+
)
77+
else:
78+
ax.scatter(
79+
w * (1 + 0.1 * m),
80+
cdnc_compare[m, k, i],
81+
color="C" + str(m),
82+
marker=mrkr[m],
83+
)
84+
ax.set_xscale("log")
85+
86+
ax.set_title(subplot, loc="left", weight="bold")
87+
ax.set_xlabel("Updraft velocity, w [m s$^{-1}$]")
88+
if k == 0:
89+
ax.set_ylabel("CDNC [cm$^{-3}$]")
90+
ax.legend()
91+
92+
pyplot.suptitle("zero organics")
93+
94+
if plot:
95+
pyplot.show()
96+
97+
np.testing.assert_allclose(
98+
cdnc_compare[0, :, :],
99+
cdnc_compare[1, :, :],
100+
rtol=1e-2,
101+
)

‎tests/unit_tests/backends/test_freezing_methods.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -124,7 +124,7 @@ def _plot_fit(fit_x, fit_y, low, hgh, total_time):
124124
pylab.legend()
125125
pylab.yscale("log")
126126
pylab.ylim(fit_y[-1], fit_y[0])
127-
pylab.xlim(0, total_time)
127+
pylab.xlim(None, total_time)
128128
pylab.xlabel("time")
129129
pylab.ylabel("unfrozen fraction")
130130
pylab.grid()

‎tests/unit_tests/physics/test_surface_tension.py

Lines changed: 91 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,91 @@
1+
# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring
2+
import warnings
3+
4+
import matplotlib.pyplot as plt
5+
import numpy as np
6+
from numba.core.errors import NumbaExperimentalFeatureWarning
7+
8+
from PySDM import Formulae
9+
from PySDM.physics import si
10+
11+
12+
def test_surface_tension(plot=False):
13+
ls = ["-", ":", "--", "-."]
14+
15+
st_type = (
16+
"Constant",
17+
"CompressedFilmOvadnevaite",
18+
"SzyszkowskiLangmuir",
19+
"CompressedFilmRuehl",
20+
)
21+
# Arrange
22+
TRIVIA = Formulae().trivia
23+
R_WET = np.logspace(np.log(100 * si.nm), np.log(1000 * si.nm), base=np.e, num=100)
24+
R_DRY = 50 * si.nm
25+
V_WET = TRIVIA.volume(R_WET)
26+
V_DRY = TRIVIA.volume(R_DRY)
27+
TEMPERATURE = 300 * si.K
28+
29+
################
30+
# test zero organic
31+
################
32+
33+
F_ORG = 0.0
34+
sgm = calc_sigma(st_type, TEMPERATURE, V_WET, V_DRY, F_ORG)
35+
36+
for i, st in enumerate(st_type):
37+
np.testing.assert_allclose(sgm[0, :], sgm[i, :], rtol=1e-6)
38+
plt.plot(R_WET / si.um, sgm[i, :], ls[i], label=st)
39+
40+
if plot:
41+
plt.title(f"Forg = {F_ORG}")
42+
plt.xlabel("wet radius (um)")
43+
plt.xscale("log")
44+
plt.ylabel("surface tension [N/m]")
45+
plt.legend()
46+
plt.show()
47+
48+
################
49+
# test all organic
50+
################
51+
52+
F_ORG = 1.0
53+
sgm = calc_sigma(st_type, TEMPERATURE, V_WET, V_DRY, F_ORG)
54+
55+
for i, st in enumerate(st_type):
56+
if i > 0:
57+
np.testing.assert_array_less(sgm[i, :], sgm[0, :])
58+
plt.plot(R_WET / si.um, sgm[i, :], ls[i], label=st)
59+
60+
if plot:
61+
plt.title(f"Forg = {F_ORG}")
62+
plt.xlabel("wet radius (um)")
63+
plt.xscale("log")
64+
plt.ylabel("surface tension [N/m]")
65+
plt.legend()
66+
plt.show()
67+
68+
69+
def calc_sigma(st_type, TEMPERATURE, V_WET, V_DRY, F_ORG):
70+
sgm = np.zeros((len(st_type), len(V_WET)))
71+
72+
with warnings.catch_warnings():
73+
warnings.simplefilter("ignore", category=NumbaExperimentalFeatureWarning)
74+
for i, st in enumerate(st_type):
75+
f = Formulae(
76+
surface_tension=st,
77+
constants={
78+
"sgm_org": 10 * si.mN / si.m,
79+
"delta_min": 1 * si.nm,
80+
"RUEHL_A0": 1e-17 * si.m * si.m,
81+
"RUEHL_C0": 1e-8,
82+
"RUEHL_m_sigma": 1e17 * si.J / si.m**2,
83+
"RUEHL_sgm_min": 10 * si.mN / si.m,
84+
"RUEHL_nu_org": 1e2 * si.cm**3 / si.mole,
85+
},
86+
)
87+
88+
for j, vw in enumerate(V_WET):
89+
sgm[i, j] = f.surface_tension.sigma(TEMPERATURE, vw, V_DRY, F_ORG)
90+
91+
return sgm

0 commit comments

Comments
 (0)