From a978eaa3a60408e1e1f44fa4e619fa2573ce9114 Mon Sep 17 00:00:00 2001 From: Oleksii Bulenok Date: Fri, 8 Jul 2022 11:59:40 +0200 Subject: [PATCH 01/23] Change thd_tol to RH_tol --- .../methods/condensation_methods.py | 43 +++++++++++-------- PySDM/backends/impl_numba/test_helpers/bdf.py | 4 +- .../methods/condensation_methods.py | 2 +- PySDM/dynamics/condensation.py | 18 +++++--- PySDM/particulator.py | 5 ++- .../test_single_supersaturation_peak.py | 38 ++++++++++------ 6 files changed, 67 insertions(+), 43 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/condensation_methods.py b/PySDM/backends/impl_numba/methods/condensation_methods.py index 94c406d75d..b83cd03d66 100644 --- a/PySDM/backends/impl_numba/methods/condensation_methods.py +++ b/PySDM/backends/impl_numba/methods/condensation_methods.py @@ -29,6 +29,7 @@ def condensation( rhod, thd, qv, + RH, dv, prhod, pthd, @@ -36,7 +37,7 @@ def condensation( kappa, f_org, rtol_x, - rtol_thd, + rtol_RH, timestep, counters, cell_order, @@ -58,6 +59,7 @@ def condensation( rhod=rhod.data, thd=thd.data, qv=qv.data, + RH=RH.data, dv_mean=dv, prhod=prhod.data, pthd=pthd.data, @@ -65,7 +67,7 @@ def condensation( kappa=kappa.data, f_org=f_org.data, rtol_x=rtol_x, - rtol_thd=rtol_thd, + rtol_RH=rtol_RH, timestep=timestep, counter_n_substeps=counters["n_substeps"].data, counter_n_activating=counters["n_activating"].data, @@ -92,6 +94,7 @@ def _condensation( rhod, thd, qv, + RH, dv_mean, prhod, pthd, @@ -99,7 +102,7 @@ def _condensation( kappa, f_org, rtol_x, - rtol_thd, + rtol_RH, timestep, counter_n_substeps, counter_n_activating, @@ -143,12 +146,13 @@ def _condensation( f_org, thd[cell_id], qv[cell_id], + RH[cell_id], dthd_dt, dqv_dt, md, rhod_mean, rtol_x, - rtol_thd, + rtol_RH, timestep, counter_n_substeps[cell_id], ) @@ -175,7 +179,7 @@ def make_adapt_substeps( n_substeps_min = math.ceil(timestep / dt_range[1]) @numba.njit(**jit_flags) - def adapt_substeps(args, n_substeps, thd, rtol_thd): + def adapt_substeps(args, n_substeps, thd, RH, rtol_RH): n_substeps = np.maximum(n_substeps_min, n_substeps // multiplier) success = False for burnout in range(fuse + 1): @@ -189,24 +193,26 @@ def adapt_substeps(args, n_substeps, thd, rtol_thd): ), return_value=(0, False), ) - thd_new_long, success = step_fake(args, timestep, n_substeps) + thd_new_long, RH_new_long, success = step_fake(args, timestep, n_substeps) if success: break n_substeps *= multiplier for burnout in range(fuse + 1): if burnout == fuse: return warn("burnout (short)", __file__, return_value=(0, False)) - thd_new_short, success = step_fake( + thd_new_short, RH_new_short, success = step_fake( args, timestep, n_substeps * multiplier ) if not success: return warn("short failed", __file__, return_value=(0, False)) - dthd_long = thd_new_long - thd - dthd_short = thd_new_short - thd - error_estimate = np.abs(dthd_long - multiplier * dthd_short) - thd_new_long = thd_new_short - if within_tolerance(error_estimate, thd, rtol_thd): + + dRH_long = RH_new_long - RH + dRH_short = RH_new_short - RH + error_estimate = np.abs(dRH_long - multiplier * dRH_short) + RH_new_long = RH_new_short + if within_tolerance(error_estimate, RH, rtol_RH): break + n_substeps *= multiplier if n_substeps > n_substeps_max: break @@ -219,8 +225,8 @@ def make_step_fake(jit_flags, step_impl): @numba.njit(**jit_flags) def step_fake(args, dt, n_substeps): dt /= n_substeps - _, thd_new, _, _, _, _, success = step_impl(*args, dt, 1, True) - return thd_new, success + _, thd_new, RH_new, _, _, _, _, success = step_impl(*args, dt, 1, True) + return thd_new, RH_new, success return step_fake @@ -328,6 +334,7 @@ def step_impl( # pylint: disable=too-many-arguments return ( qv, thd, + RH, count_activating, count_deactivating, count_ripening, @@ -412,7 +419,7 @@ def calculate_ml_new( # pylint: disable=too-many-arguments lambdaK = phys_lambdaK(T, p) lambdaD = phys_lambdaD(DTp, T) for drop in cell_idx: - if v[drop] < 0: + if v[drop] < 0: # TODO: use _unfrozen from freezing_methods.py continue x_old = x(v[drop]) r_old = radius(v[drop]) @@ -659,12 +666,13 @@ def solve( # pylint: disable=too-many-arguments f_org, thd, qv, + RH, dthd_dt, dqv_dt, m_d, rhod_mean, rtol_x, - rtol_thd, + rtol_RH, timestep, n_substeps, ): @@ -686,11 +694,12 @@ def solve( # pylint: disable=too-many-arguments ) success = True if adaptive: - n_substeps, success = adapt_substeps(args, n_substeps, thd, rtol_thd) + n_substeps, success = adapt_substeps(args, n_substeps, thd, RH, rtol_RH) if success: ( qv, thd, + _, n_activating, n_deactivating, n_ripening, diff --git a/PySDM/backends/impl_numba/test_helpers/bdf.py b/PySDM/backends/impl_numba/test_helpers/bdf.py index 65d61fc203..319c9e6a34 100644 --- a/PySDM/backends/impl_numba/test_helpers/bdf.py +++ b/PySDM/backends/impl_numba/test_helpers/bdf.py @@ -25,7 +25,7 @@ def patch_particulator(particulator): def _bdf_condensation( - particulator, *, rtol_x, rtol_thd, counters, RH_max, success, cell_order + particulator, *, rtol_x, rtol_RH, counters, RH_max, success, cell_order ): func = Numba._condensation if not numba.config.DISABLE_JIT: # pylint: disable=no-member @@ -50,7 +50,7 @@ def _bdf_condensation( kappa=particulator.attributes["kappa"].data, f_org=particulator.attributes["dry volume organic fraction"].data, rtol_x=rtol_x, - rtol_thd=rtol_thd, + rtol_RH=rtol_RH, timestep=particulator.dt, counter_n_substeps=counters["n_substeps"], counter_n_activating=counters["n_activating"], diff --git a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py index 75086382b7..635ace699a 100644 --- a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py +++ b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py @@ -325,7 +325,7 @@ def condensation( kappa, f_org, rtol_x, - rtol_thd, + rtol_RH, timestep, counters, cell_order, diff --git a/PySDM/dynamics/condensation.py b/PySDM/dynamics/condensation.py index 27aed8ed68..7c66017553 100644 --- a/PySDM/dynamics/condensation.py +++ b/PySDM/dynamics/condensation.py @@ -8,9 +8,14 @@ from ..physics import si -DEFAULTS = namedtuple("_", ("rtol_x", "rtol_thd", "cond_range", "schedule"))( +DEFAULTS = namedtuple("_", ( + "rtol_x", + "rtol_RH", + "cond_range", + "schedule" +))( rtol_x=1e-6, - rtol_thd=1e-6, + rtol_RH=1e-4, cond_range=(1e-4 * si.second, 1 * si.second), schedule="dynamic", ) @@ -21,7 +26,7 @@ def __init__( self, *, rtol_x=DEFAULTS.rtol_x, - rtol_thd=DEFAULTS.rtol_thd, + rtol_RH=DEFAULTS.rtol_RH, substeps: int = 1, adaptive: bool = True, dt_cond_range: tuple = DEFAULTS.cond_range, @@ -29,12 +34,11 @@ def __init__( max_iters: int = 16, update_thd: bool = True, ): - self.particulator = None self.enable = True self.rtol_x = rtol_x - self.rtol_thd = rtol_thd + self.rtol_RH = rtol_RH self.rh_max = None self.success = None @@ -58,7 +62,7 @@ def register(self, builder): adaptive=self.adaptive, fuse=32, multiplier=2, - RH_rtol=1e-7, + RH_rtol=1e-7, # TODO: use rtol_RH? max_iters=self.max_iters, ) builder.request_attribute("critical volume") @@ -95,7 +99,7 @@ def __call__(self): self.particulator.condensation( rtol_x=self.rtol_x, - rtol_thd=self.rtol_thd, + rtol_RH=self.rtol_RH, counters=self.counters, RH_max=self.rh_max, success=self.success, diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 72896d1235..31eddd74e3 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -101,7 +101,7 @@ def update_TpRH(self): RH=self.environment.get_predicted("RH"), ) - def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_order): + def condensation(self, *, rtol_x, rtol_RH, counters, RH_max, success, cell_order): self.backend.condensation( solver=self.condensation_solver, n_cell=self.mesh.n_cell, @@ -113,6 +113,7 @@ def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_orde rhod=self.environment["rhod"], thd=self.environment["thd"], qv=self.environment["qv"], + RH=self.environment["RH"], dv=self.environment.dv, prhod=self.environment.get_predicted("rhod"), pthd=self.environment.get_predicted("thd"), @@ -120,7 +121,7 @@ def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_orde kappa=self.attributes["kappa"], f_org=self.attributes["dry volume organic fraction"], rtol_x=rtol_x, - rtol_thd=rtol_thd, + rtol_RH=rtol_RH, v_cr=self.attributes["critical volume"], timestep=self.dt, counters=counters, diff --git a/tests/smoke_tests/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py b/tests/smoke_tests/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py index 55910b1309..8c0ff7df6d 100644 --- a/tests/smoke_tests/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py +++ b/tests/smoke_tests/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py @@ -17,18 +17,22 @@ @pytest.mark.parametrize( - "rtol_thd", + "rtol_RH", ( - pytest.param(1e-6, marks=pytest.mark.xfail(strict=True)), - pytest.param(1e-7, marks=pytest.mark.xfail(strict=True)), - 1e-8, - 1e-9, + pytest.param(1e-2, marks=pytest.mark.xfail(strict=True)), + pytest.param(1e-3, marks=pytest.mark.xfail(strict=True)), + 1e-4, + 1e-5, ), ) @pytest.mark.parametrize("rtol_x", (1e-7,)) -@pytest.mark.parametrize("adaptive", (True,)) @pytest.mark.parametrize("scheme", ("PySDM",)) -def test_single_supersaturation_peak(adaptive, scheme, rtol_x, rtol_thd, plot=False): +def test_single_supersaturation_peak( + scheme, + rtol_x, + rtol_RH, + plot=False +): # arrange products = ( PySDM_products.WaterMixingRatio(unit="g/kg", name="ql"), @@ -36,15 +40,18 @@ def test_single_supersaturation_peak(adaptive, scheme, rtol_x, rtol_thd, plot=Fa PySDM_products.AmbientRelativeHumidity(name="RH"), PySDM_products.ParcelDisplacement(name="z"), ) + dt = 2 * si.s + w = 0.5 * si.m / si.s env = Parcel( - dt=2 * si.s, + dt=dt, mass_of_dry_air=1e3 * si.kg, p0=1000 * si.hPa, q0=22.76 * si.g / si.kg, - w=0.5 * si.m / si.s, + w=w, T0=300 * si.K, ) - n_steps = 70 + z_max = 70 * si.m + n_steps = int(z_max / (w * dt)) n_sd = 2 kappa = 0.4 spectrum = Lognormal(norm_factor=5000 / si.cm**3, m_mode=50.0 * si.nm, s_geom=2.0) @@ -53,9 +60,9 @@ def test_single_supersaturation_peak(adaptive, scheme, rtol_x, rtol_thd, plot=Fa builder.add_dynamic(AmbientThermodynamics()) builder.add_dynamic( Condensation( - adaptive=adaptive, + adaptive=True, rtol_x=rtol_x, - rtol_thd=rtol_thd, + rtol_RH=rtol_RH ) ) @@ -108,9 +115,12 @@ def test_single_supersaturation_peak(adaptive, scheme, rtol_x, rtol_thd, plot=Fa twin.set_xlim(-0.001, 0.0015) pyplot.legend(loc="lower right") pyplot.grid() - pyplot.title(f"rtol_thd={rtol_thd}; rtol_x={rtol_x}") + pyplot.title(f"rtol_RH={rtol_RH}; rtol_x={rtol_x}") if plot: pyplot.show() # assert - assert signal.argrelextrema(np.asarray(output["RH"]), np.greater)[0].shape[0] == 1 + n_max = signal.argrelextrema(np.asarray(output["RH"]), np.greater)[0].shape[0] + n_min = signal.argrelextrema(np.asarray(output["RH"]), np.less)[0].shape[0] + assert n_max == 1 and n_min == 0 + From df1d9fb318859d2036aec0ef1c3daa428d6e4b88 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sat, 9 Jul 2022 00:15:55 +0200 Subject: [PATCH 02/23] making pylint & black happy; bumping test-time-req for examples --- .../impl_numba/methods/condensation_methods.py | 4 ++-- PySDM/dynamics/condensation.py | 7 +------ test-time-requirements.txt | 2 +- .../test_single_supersaturation_peak.py | 16 ++-------------- 4 files changed, 6 insertions(+), 23 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/condensation_methods.py b/PySDM/backends/impl_numba/methods/condensation_methods.py index b83cd03d66..a3d1928551 100644 --- a/PySDM/backends/impl_numba/methods/condensation_methods.py +++ b/PySDM/backends/impl_numba/methods/condensation_methods.py @@ -193,14 +193,14 @@ def adapt_substeps(args, n_substeps, thd, RH, rtol_RH): ), return_value=(0, False), ) - thd_new_long, RH_new_long, success = step_fake(args, timestep, n_substeps) + _, RH_new_long, success = step_fake(args, timestep, n_substeps) if success: break n_substeps *= multiplier for burnout in range(fuse + 1): if burnout == fuse: return warn("burnout (short)", __file__, return_value=(0, False)) - thd_new_short, RH_new_short, success = step_fake( + _, RH_new_short, success = step_fake( args, timestep, n_substeps * multiplier ) if not success: diff --git a/PySDM/dynamics/condensation.py b/PySDM/dynamics/condensation.py index 7c66017553..faef14ca82 100644 --- a/PySDM/dynamics/condensation.py +++ b/PySDM/dynamics/condensation.py @@ -8,12 +8,7 @@ from ..physics import si -DEFAULTS = namedtuple("_", ( - "rtol_x", - "rtol_RH", - "cond_range", - "schedule" -))( +DEFAULTS = namedtuple("_", ("rtol_x", "rtol_RH", "cond_range", "schedule"))( rtol_x=1e-6, rtol_RH=1e-4, cond_range=(1e-4 * si.second, 1 * si.second), diff --git a/test-time-requirements.txt b/test-time-requirements.txt index b72e22fb21..2a410f5439 100644 --- a/test-time-requirements.txt +++ b/test-time-requirements.txt @@ -4,6 +4,6 @@ ghapi pytest # note: if cloning both PySDM and PySDM examples, consider "pip install -e" -PySDM-examples @ git+https://github.com/atmos-cloud-sim-uj/PySDM-examples@4b1fed8#egg=PySDM-examples +PySDM-examples @ git+https://github.com/slayoo/PySDM-examples@30601b6#egg=PySDM-examples PyMPDATA @ git+https://github.com/atmos-cloud-sim-uj/PyMPDATA@e7b73a7#egg=PyMPDATA PyPartMC==0.0.9; platform_machine == 'x86_64' diff --git a/tests/smoke_tests/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py b/tests/smoke_tests/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py index 8c0ff7df6d..7ca64ccab4 100644 --- a/tests/smoke_tests/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py +++ b/tests/smoke_tests/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py @@ -27,12 +27,7 @@ ) @pytest.mark.parametrize("rtol_x", (1e-7,)) @pytest.mark.parametrize("scheme", ("PySDM",)) -def test_single_supersaturation_peak( - scheme, - rtol_x, - rtol_RH, - plot=False -): +def test_single_supersaturation_peak(scheme, rtol_x, rtol_RH, plot=False): # arrange products = ( PySDM_products.WaterMixingRatio(unit="g/kg", name="ql"), @@ -58,13 +53,7 @@ def test_single_supersaturation_peak( builder = Builder(backend=CPU(), n_sd=n_sd) builder.set_environment(env) builder.add_dynamic(AmbientThermodynamics()) - builder.add_dynamic( - Condensation( - adaptive=True, - rtol_x=rtol_x, - rtol_RH=rtol_RH - ) - ) + builder.add_dynamic(Condensation(adaptive=True, rtol_x=rtol_x, rtol_RH=rtol_RH)) r_dry, concentration = ConstantMultiplicity(spectrum).sample(n_sd) v_dry = builder.formulae.trivia.volume(radius=r_dry) @@ -123,4 +112,3 @@ def test_single_supersaturation_peak( n_max = signal.argrelextrema(np.asarray(output["RH"]), np.greater)[0].shape[0] n_min = signal.argrelextrema(np.asarray(output["RH"]), np.less)[0].shape[0] assert n_max == 1 and n_min == 0 - From dde3ab55f7b857fd975dfb8c94a41fe4365d88ba Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sat, 9 Jul 2022 00:34:11 +0200 Subject: [PATCH 03/23] fix GPU backend condensation signature --- PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py | 1 + 1 file changed, 1 insertion(+) diff --git a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py index 635ace699a..8044841bad 100644 --- a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py +++ b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py @@ -318,6 +318,7 @@ def condensation( rhod, thd, qv, + RH, dv, prhod, pthd, From ba9fb9a6f56f43c0003f18b4b2099a7d243b03b4 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sat, 9 Jul 2022 01:08:04 +0200 Subject: [PATCH 04/23] fix BDF solver API --- PySDM/backends/impl_numba/test_helpers/bdf.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/PySDM/backends/impl_numba/test_helpers/bdf.py b/PySDM/backends/impl_numba/test_helpers/bdf.py index 319c9e6a34..c74aa2a1e8 100644 --- a/PySDM/backends/impl_numba/test_helpers/bdf.py +++ b/PySDM/backends/impl_numba/test_helpers/bdf.py @@ -43,6 +43,7 @@ def _bdf_condensation( rhod=particulator.environment["rhod"].data, thd=particulator.environment["thd"].data, qv=particulator.environment["qv"].data, + RH=particulator.environment["RH"].data, dv_mean=particulator.environment.dv, prhod=particulator.environment.get_predicted("rhod").data, pthd=particulator.environment.get_predicted("thd").data, @@ -177,14 +178,15 @@ def solve( # pylint: disable=too-many-arguments f_org, thd, qv, + __, dthd_dt, dqv_dt, m_d_mean, rhod_mean, - __, ___, - dt, ____, + dt, + _____, ): n_sd_in_cell = len(cell_idx) y0 = np.empty(n_sd_in_cell + idx_x) From 2023c904c71193599b7423efbc8d54ae13f048a1 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sat, 9 Jul 2022 21:00:27 +0200 Subject: [PATCH 05/23] default rtol_RH=1e-6; bump examples req; reduce significant digit number in arabas_and_shima_2017/test_conservation --- PySDM/dynamics/condensation.py | 2 +- test-time-requirements.txt | 2 +- tests/smoke_tests/arabas_and_shima_2017/test_conservation.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/PySDM/dynamics/condensation.py b/PySDM/dynamics/condensation.py index faef14ca82..d60c4bedb5 100644 --- a/PySDM/dynamics/condensation.py +++ b/PySDM/dynamics/condensation.py @@ -10,7 +10,7 @@ DEFAULTS = namedtuple("_", ("rtol_x", "rtol_RH", "cond_range", "schedule"))( rtol_x=1e-6, - rtol_RH=1e-4, + rtol_RH=1e-6, cond_range=(1e-4 * si.second, 1 * si.second), schedule="dynamic", ) diff --git a/test-time-requirements.txt b/test-time-requirements.txt index 2a410f5439..625c4dc05f 100644 --- a/test-time-requirements.txt +++ b/test-time-requirements.txt @@ -4,6 +4,6 @@ ghapi pytest # note: if cloning both PySDM and PySDM examples, consider "pip install -e" -PySDM-examples @ git+https://github.com/slayoo/PySDM-examples@30601b6#egg=PySDM-examples +PySDM-examples @ git+https://github.com/slayoo/PySDM-examples@9c8d1cd#egg=PySDM-examples PyMPDATA @ git+https://github.com/atmos-cloud-sim-uj/PyMPDATA@e7b73a7#egg=PyMPDATA PyPartMC==0.0.9; platform_machine == 'x86_64' diff --git a/tests/smoke_tests/arabas_and_shima_2017/test_conservation.py b/tests/smoke_tests/arabas_and_shima_2017/test_conservation.py index b4f41be0ed..197603cfb1 100644 --- a/tests/smoke_tests/arabas_and_shima_2017/test_conservation.py +++ b/tests/smoke_tests/arabas_and_shima_2017/test_conservation.py @@ -47,7 +47,7 @@ def test_water_mass_conservation(settings_idx, mass_of_dry_air, scheme, coord): # Assert ql = liquid_water_mixing_ratio(simulation) qt = simulation.particulator.environment["qv"].to_ndarray() + ql - significant = 6 if scheme == "GPU" else 14 # TODO #540 + significant = 6 if scheme == "GPU" else 12 # TODO #540 np.testing.assert_approx_equal(qt, qt0, significant) if scheme != "BDF": assert simulation.particulator.products["S_max"].get() >= output["S"][-1] From 51e7138aa9803bb7e3d8b48455e74ec37309db6d Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Mon, 8 May 2023 13:54:46 +0200 Subject: [PATCH 06/23] point to correct branch in ttr --- test-time-requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test-time-requirements.txt b/test-time-requirements.txt index 08dc5e6a48..b55d0f4a77 100644 --- a/test-time-requirements.txt +++ b/test-time-requirements.txt @@ -11,7 +11,7 @@ pytest pytest-timeout # note: if cloning both PySDM and PySDM examples, consider "pip install -e" -PySDM-examples @ git+https://github.com/edejong-caltech/PySDM-examples@4bf5122#egg=PySDM-examples +PySDM-examples @ git+https://github.com/slayoo/PySDM-examples@13a5fe3#egg=PySDM-examples PyMPDATA @ git+https://github.com/open-atmos/PyMPDATA@63f1c11#egg=PyMPDATA PyPartMC==0.0.32 From d2771d269b7e888673eaf3bb85e59102726937cf Mon Sep 17 00:00:00 2001 From: Oleksii Bulenok Date: Fri, 8 Jul 2022 11:59:40 +0200 Subject: [PATCH 07/23] Change thd_tol to RH_tol --- .../methods/condensation_methods.py | 43 +++++++++++-------- .../scipy_ode_condensation_solver.py | 4 +- .../methods/condensation_methods.py | 2 +- PySDM/dynamics/condensation.py | 17 +++++--- PySDM/particulator.py | 5 ++- .../test_single_supersaturation_peak.py | 38 +++++++++------- 6 files changed, 66 insertions(+), 43 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/condensation_methods.py b/PySDM/backends/impl_numba/methods/condensation_methods.py index 90b2e4f6fb..febfe708cc 100644 --- a/PySDM/backends/impl_numba/methods/condensation_methods.py +++ b/PySDM/backends/impl_numba/methods/condensation_methods.py @@ -29,6 +29,7 @@ def condensation( rhod, thd, qv, + RH, dv, prhod, pthd, @@ -36,7 +37,7 @@ def condensation( kappa, f_org, rtol_x, - rtol_thd, + rtol_RH, timestep, counters, cell_order, @@ -58,6 +59,7 @@ def condensation( rhod=rhod.data, thd=thd.data, qv=qv.data, + RH=RH.data, dv_mean=dv, prhod=prhod.data, pthd=pthd.data, @@ -65,7 +67,7 @@ def condensation( kappa=kappa.data, f_org=f_org.data, rtol_x=rtol_x, - rtol_thd=rtol_thd, + rtol_RH=rtol_RH, timestep=timestep, counter_n_substeps=counters["n_substeps"].data, counter_n_activating=counters["n_activating"].data, @@ -92,6 +94,7 @@ def _condensation( rhod, thd, qv, + RH, dv_mean, prhod, pthd, @@ -99,7 +102,7 @@ def _condensation( kappa, f_org, rtol_x, - rtol_thd, + rtol_RH, timestep, counter_n_substeps, counter_n_activating, @@ -144,12 +147,13 @@ def _condensation( f_org, thd[cell_id], qv[cell_id], + RH[cell_id], dthd_dt, dqv_dt, md, rhod_mean, rtol_x, - rtol_thd, + rtol_RH, timestep, counter_n_substeps[cell_id], ) @@ -176,7 +180,7 @@ def make_adapt_substeps( n_substeps_min = math.ceil(timestep / dt_range[1]) @numba.njit(**jit_flags) - def adapt_substeps(args, n_substeps, thd, rtol_thd): + def adapt_substeps(args, n_substeps, thd, RH, rtol_RH): n_substeps = np.maximum(n_substeps_min, n_substeps // multiplier) success = False for burnout in range(fuse + 1): @@ -190,24 +194,26 @@ def adapt_substeps(args, n_substeps, thd, rtol_thd): ), return_value=(0, False), ) - thd_new_long, success = step_fake(args, timestep, n_substeps) + thd_new_long, RH_new_long, success = step_fake(args, timestep, n_substeps) if success: break n_substeps *= multiplier for burnout in range(fuse + 1): if burnout == fuse: return warn("burnout (short)", __file__, return_value=(0, False)) - thd_new_short, success = step_fake( + thd_new_short, RH_new_short, success = step_fake( args, timestep, n_substeps * multiplier ) if not success: return warn("short failed", __file__, return_value=(0, False)) - dthd_long = thd_new_long - thd - dthd_short = thd_new_short - thd - error_estimate = np.abs(dthd_long - multiplier * dthd_short) - thd_new_long = thd_new_short - if within_tolerance(error_estimate, thd, rtol_thd): + + dRH_long = RH_new_long - RH + dRH_short = RH_new_short - RH + error_estimate = np.abs(dRH_long - multiplier * dRH_short) + RH_new_long = RH_new_short + if within_tolerance(error_estimate, RH, rtol_RH): break + n_substeps *= multiplier if n_substeps > n_substeps_max: break @@ -220,8 +226,8 @@ def make_step_fake(jit_flags, step_impl): @numba.njit(**jit_flags) def step_fake(args, dt, n_substeps): dt /= n_substeps - _, thd_new, _, _, _, _, success = step_impl(*args, dt, 1, True) - return thd_new, success + _, thd_new, RH_new, _, _, _, _, success = step_impl(*args, dt, 1, True) + return thd_new, RH_new, success return step_fake @@ -329,6 +335,7 @@ def step_impl( # pylint: disable=too-many-arguments,too-many-locals return ( qv, thd, + RH, count_activating, count_deactivating, count_ripening, @@ -413,7 +420,7 @@ def calculate_ml_new( # pylint: disable=too-many-arguments,too-many-statements, lambdaK = phys_lambdaK(T, p) lambdaD = phys_lambdaD(DTp, T) for drop in cell_idx: - if v[drop] < 0: + if v[drop] < 0: # TODO: use _unfrozen from freezing_methods.py continue x_old = x(v[drop]) r_old = radius(v[drop]) @@ -661,12 +668,13 @@ def solve( # pylint: disable=too-many-arguments f_org, thd, qv, + RH, dthd_dt, dqv_dt, m_d, rhod_mean, rtol_x, - rtol_thd, + rtol_RH, timestep, n_substeps, ): @@ -688,11 +696,12 @@ def solve( # pylint: disable=too-many-arguments ) success = True if adaptive: - n_substeps, success = adapt_substeps(args, n_substeps, thd, rtol_thd) + n_substeps, success = adapt_substeps(args, n_substeps, thd, RH, rtol_RH) if success: ( qv, thd, + _, n_activating, n_deactivating, n_ripening, diff --git a/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py b/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py index 80bf1df5d8..a9e9c62009 100644 --- a/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py +++ b/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py @@ -24,7 +24,7 @@ def patch_particulator(particulator): def _condensation( - particulator, *, rtol_x, rtol_thd, counters, RH_max, success, cell_order + particulator, *, rtol_x, rtol_RH, counters, RH_max, success, cell_order ): func = Numba._condensation if not numba.config.DISABLE_JIT: # pylint: disable=no-member @@ -49,7 +49,7 @@ def _condensation( kappa=particulator.attributes["kappa"].data, f_org=particulator.attributes["dry volume organic fraction"].data, rtol_x=rtol_x, - rtol_thd=rtol_thd, + rtol_RH=rtol_RH, timestep=particulator.dt, counter_n_substeps=counters["n_substeps"], counter_n_activating=counters["n_activating"], diff --git a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py index 63db0e87bc..7f48825b29 100644 --- a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py +++ b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py @@ -327,7 +327,7 @@ def condensation( kappa, f_org, rtol_x, - rtol_thd, + rtol_RH, timestep, counters, cell_order, diff --git a/PySDM/dynamics/condensation.py b/PySDM/dynamics/condensation.py index fbebef962d..eba777875a 100644 --- a/PySDM/dynamics/condensation.py +++ b/PySDM/dynamics/condensation.py @@ -8,9 +8,14 @@ from ..physics import si -DEFAULTS = namedtuple("_", ("rtol_x", "rtol_thd", "cond_range", "schedule"))( +DEFAULTS = namedtuple("_", ( + "rtol_x", + "rtol_RH", + "cond_range", + "schedule" +))( rtol_x=1e-6, - rtol_thd=1e-6, + rtol_RH=1e-4, cond_range=(1e-4 * si.second, 1 * si.second), schedule="dynamic", ) @@ -21,7 +26,7 @@ def __init__( self, *, rtol_x=DEFAULTS.rtol_x, - rtol_thd=DEFAULTS.rtol_thd, + rtol_RH=DEFAULTS.rtol_RH, substeps: int = 1, adaptive: bool = True, dt_cond_range: tuple = DEFAULTS.cond_range, @@ -33,7 +38,7 @@ def __init__( self.enable = True self.rtol_x = rtol_x - self.rtol_thd = rtol_thd + self.rtol_RH = rtol_RH self.rh_max = None self.success = None @@ -57,7 +62,7 @@ def register(self, builder): adaptive=self.adaptive, fuse=32, multiplier=2, - RH_rtol=1e-7, + RH_rtol=1e-7, # TODO: use rtol_RH? max_iters=self.max_iters, ) builder.request_attribute("critical volume") @@ -94,7 +99,7 @@ def __call__(self): self.particulator.condensation( rtol_x=self.rtol_x, - rtol_thd=self.rtol_thd, + rtol_RH=self.rtol_RH, counters=self.counters, RH_max=self.rh_max, success=self.success, diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 2e1055d478..93e1f33f1d 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -101,7 +101,7 @@ def update_TpRH(self): RH=self.environment.get_predicted("RH"), ) - def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_order): + def condensation(self, *, rtol_x, rtol_RH, counters, RH_max, success, cell_order): self.backend.condensation( solver=self.condensation_solver, n_cell=self.mesh.n_cell, @@ -113,6 +113,7 @@ def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_orde rhod=self.environment["rhod"], thd=self.environment["thd"], qv=self.environment["qv"], + RH=self.environment["RH"], dv=self.environment.dv, prhod=self.environment.get_predicted("rhod"), pthd=self.environment.get_predicted("thd"), @@ -120,7 +121,7 @@ def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_orde kappa=self.attributes["kappa"], f_org=self.attributes["dry volume organic fraction"], rtol_x=rtol_x, - rtol_thd=rtol_thd, + rtol_RH=rtol_RH, v_cr=self.attributes["critical volume"], timestep=self.dt, counters=counters, diff --git a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py index 62c9296efd..dfb3befebb 100644 --- a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py +++ b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py @@ -17,20 +17,22 @@ @pytest.mark.parametrize( - "rtol_thd", + "rtol_RH", ( - pytest.param(1e-6, marks=pytest.mark.xfail(strict=True)), - pytest.param(1e-7, marks=pytest.mark.xfail(strict=True)), - 1e-8, - 1e-9, + pytest.param(1e-2, marks=pytest.mark.xfail(strict=True)), + pytest.param(1e-3, marks=pytest.mark.xfail(strict=True)), + 1e-4, + 1e-5, ), ) @pytest.mark.parametrize("rtol_x", (1e-7,)) -@pytest.mark.parametrize("adaptive", (True,)) @pytest.mark.parametrize("scheme", ("PySDM",)) def test_single_supersaturation_peak( - adaptive, scheme, rtol_x, rtol_thd, plot=False -): # pylint: disable=too-many-locals + scheme, + rtol_x, + rtol_RH, + plot=False +): # arrange products = ( PySDM_products.WaterMixingRatio(unit="g/kg", name="ql"), @@ -38,15 +40,18 @@ def test_single_supersaturation_peak( PySDM_products.AmbientRelativeHumidity(name="RH"), PySDM_products.ParcelDisplacement(name="z"), ) + dt = 2 * si.s + w = 0.5 * si.m / si.s env = Parcel( - dt=2 * si.s, + dt=dt, mass_of_dry_air=1e3 * si.kg, p0=1000 * si.hPa, q0=22.76 * si.g / si.kg, - w=0.5 * si.m / si.s, + w=w, T0=300 * si.K, ) - n_steps = 70 + z_max = 70 * si.m + n_steps = int(z_max / (w * dt)) n_sd = 2 kappa = 0.4 spectrum = Lognormal(norm_factor=5000 / si.cm**3, m_mode=50.0 * si.nm, s_geom=2.0) @@ -55,9 +60,9 @@ def test_single_supersaturation_peak( builder.add_dynamic(AmbientThermodynamics()) builder.add_dynamic( Condensation( - adaptive=adaptive, + adaptive=True, rtol_x=rtol_x, - rtol_thd=rtol_thd, + rtol_RH=rtol_RH ) ) @@ -110,9 +115,12 @@ def test_single_supersaturation_peak( twin.set_xlim(-0.001, 0.0015) pyplot.legend(loc="lower right") pyplot.grid() - pyplot.title(f"rtol_thd={rtol_thd}; rtol_x={rtol_x}") + pyplot.title(f"rtol_RH={rtol_RH}; rtol_x={rtol_x}") if plot: pyplot.show() # assert - assert signal.argrelextrema(np.asarray(output["RH"]), np.greater)[0].shape[0] == 1 + n_max = signal.argrelextrema(np.asarray(output["RH"]), np.greater)[0].shape[0] + n_min = signal.argrelextrema(np.asarray(output["RH"]), np.less)[0].shape[0] + assert n_max == 1 and n_min == 0 + From 88ac2a511ff9b16b6c942b2414d6a4ff11bfdd8d Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sat, 9 Jul 2022 00:15:55 +0200 Subject: [PATCH 08/23] making pylint & black happy; bumping test-time-req for examples --- .../impl_numba/methods/condensation_methods.py | 4 ++-- PySDM/dynamics/condensation.py | 7 +------ .../test_single_supersaturation_peak.py | 16 ++-------------- 3 files changed, 5 insertions(+), 22 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/condensation_methods.py b/PySDM/backends/impl_numba/methods/condensation_methods.py index febfe708cc..6d65240053 100644 --- a/PySDM/backends/impl_numba/methods/condensation_methods.py +++ b/PySDM/backends/impl_numba/methods/condensation_methods.py @@ -194,14 +194,14 @@ def adapt_substeps(args, n_substeps, thd, RH, rtol_RH): ), return_value=(0, False), ) - thd_new_long, RH_new_long, success = step_fake(args, timestep, n_substeps) + _, RH_new_long, success = step_fake(args, timestep, n_substeps) if success: break n_substeps *= multiplier for burnout in range(fuse + 1): if burnout == fuse: return warn("burnout (short)", __file__, return_value=(0, False)) - thd_new_short, RH_new_short, success = step_fake( + _, RH_new_short, success = step_fake( args, timestep, n_substeps * multiplier ) if not success: diff --git a/PySDM/dynamics/condensation.py b/PySDM/dynamics/condensation.py index eba777875a..5a2f35587f 100644 --- a/PySDM/dynamics/condensation.py +++ b/PySDM/dynamics/condensation.py @@ -8,12 +8,7 @@ from ..physics import si -DEFAULTS = namedtuple("_", ( - "rtol_x", - "rtol_RH", - "cond_range", - "schedule" -))( +DEFAULTS = namedtuple("_", ("rtol_x", "rtol_RH", "cond_range", "schedule"))( rtol_x=1e-6, rtol_RH=1e-4, cond_range=(1e-4 * si.second, 1 * si.second), diff --git a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py index dfb3befebb..7f270e31d8 100644 --- a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py +++ b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py @@ -27,12 +27,7 @@ ) @pytest.mark.parametrize("rtol_x", (1e-7,)) @pytest.mark.parametrize("scheme", ("PySDM",)) -def test_single_supersaturation_peak( - scheme, - rtol_x, - rtol_RH, - plot=False -): +def test_single_supersaturation_peak(scheme, rtol_x, rtol_RH, plot=False): # arrange products = ( PySDM_products.WaterMixingRatio(unit="g/kg", name="ql"), @@ -58,13 +53,7 @@ def test_single_supersaturation_peak( builder = Builder(backend=CPU(), n_sd=n_sd) builder.set_environment(env) builder.add_dynamic(AmbientThermodynamics()) - builder.add_dynamic( - Condensation( - adaptive=True, - rtol_x=rtol_x, - rtol_RH=rtol_RH - ) - ) + builder.add_dynamic(Condensation(adaptive=True, rtol_x=rtol_x, rtol_RH=rtol_RH)) r_dry, concentration = ConstantMultiplicity(spectrum).sample(n_sd) v_dry = builder.formulae.trivia.volume(radius=r_dry) @@ -123,4 +112,3 @@ def test_single_supersaturation_peak( n_max = signal.argrelextrema(np.asarray(output["RH"]), np.greater)[0].shape[0] n_min = signal.argrelextrema(np.asarray(output["RH"]), np.less)[0].shape[0] assert n_max == 1 and n_min == 0 - From 9b732919990fe18d7e88d613043dc4c4164c234f Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sat, 9 Jul 2022 00:34:11 +0200 Subject: [PATCH 09/23] fix GPU backend condensation signature --- PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py | 1 + 1 file changed, 1 insertion(+) diff --git a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py index 7f48825b29..de1a721d8f 100644 --- a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py +++ b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py @@ -320,6 +320,7 @@ def condensation( rhod, thd, qv, + RH, dv, prhod, pthd, From 441a97a8cb0f8a9e02d720dca4ed5a9a21bc1bfb Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sat, 9 Jul 2022 01:08:04 +0200 Subject: [PATCH 10/23] fix BDF solver API --- .../test_helpers/scipy_ode_condensation_solver.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py b/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py index a9e9c62009..3b9b17261d 100644 --- a/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py +++ b/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py @@ -42,6 +42,7 @@ def _condensation( rhod=particulator.environment["rhod"].data, thd=particulator.environment["thd"].data, qv=particulator.environment["qv"].data, + RH=particulator.environment["RH"].data, dv_mean=particulator.environment.dv, prhod=particulator.environment.get_predicted("rhod").data, pthd=particulator.environment.get_predicted("thd").data, @@ -176,14 +177,15 @@ def solve( # pylint: disable=too-many-arguments,too-many-locals f_org, thd, qv, + __, dthd_dt, dqv_dt, m_d_mean, rhod_mean, - __, ___, - dt, ____, + dt, + _____, ): n_sd_in_cell = len(cell_idx) y0 = np.empty(n_sd_in_cell + idx_x) From 38694f7506d12faa2f34b59f47b05d7be09f817b Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sat, 9 Jul 2022 21:00:27 +0200 Subject: [PATCH 11/23] default rtol_RH=1e-6; bump examples req; reduce significant digit number in arabas_and_shima_2017/test_conservation --- PySDM/dynamics/condensation.py | 2 +- .../parcel/arabas_and_shima_2017/test_conservation.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/PySDM/dynamics/condensation.py b/PySDM/dynamics/condensation.py index 5a2f35587f..90069ab91f 100644 --- a/PySDM/dynamics/condensation.py +++ b/PySDM/dynamics/condensation.py @@ -10,7 +10,7 @@ DEFAULTS = namedtuple("_", ("rtol_x", "rtol_RH", "cond_range", "schedule"))( rtol_x=1e-6, - rtol_RH=1e-4, + rtol_RH=1e-6, cond_range=(1e-4 * si.second, 1 * si.second), schedule="dynamic", ) diff --git a/tests/smoke_tests/parcel/arabas_and_shima_2017/test_conservation.py b/tests/smoke_tests/parcel/arabas_and_shima_2017/test_conservation.py index 462b1f9a8c..e5d40d4254 100644 --- a/tests/smoke_tests/parcel/arabas_and_shima_2017/test_conservation.py +++ b/tests/smoke_tests/parcel/arabas_and_shima_2017/test_conservation.py @@ -47,7 +47,7 @@ def test_water_mass_conservation(settings_idx, mass_of_dry_air, scheme, coord): # Assert ql = liquid_water_mixing_ratio(simulation) qt = simulation.particulator.environment["qv"].to_ndarray() + ql - significant = 6 if scheme == "GPU" else 14 # TODO #540 + significant = 6 if scheme == "GPU" else 12 # TODO #540 np.testing.assert_approx_equal(qt, qt0, significant) if scheme != "SciPy": assert simulation.particulator.products["S_max"].get() >= output["S"][-1] From 9ad06a204a02a58543e62181625ca656e829c196 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 7 Jun 2023 22:20:27 +0200 Subject: [PATCH 12/23] updating examples code --- .../Arabas_and_Shima_2017/settings.py | 2 +- .../Arabas_and_Shima_2017/simulation.py | 2 +- .../fig_5_SCIPY_VS_ADAPTIVE.py | 2 +- .../Bartman_et_al_2021/demo.ipynb | 2 +- .../Bartman_et_al_2021/demo_fig2.ipynb | 10 +++++----- .../Bartman_et_al_2021/demo_fig3.ipynb | 2 +- .../Morrison_and_Grabowski_2007/common.py | 2 +- examples/PySDM_examples/Pyrcel/simulation.py | 4 ++-- .../Shipway_and_Hill_2012/settings.py | 2 +- .../Shipway_and_Hill_2012/simulation.py | 2 +- .../Szumowski_et_al_1998/gui_settings.py | 18 +++++++++--------- .../Szumowski_et_al_1998/simulation.py | 2 +- .../PySDM_examples/Yang_et_al_2018/fig_2.ipynb | 10 +++++----- .../PySDM_examples/Yang_et_al_2018/settings.py | 2 +- .../Yang_et_al_2018/simulation.py | 2 +- tests/devops_tests | 2 +- 16 files changed, 33 insertions(+), 33 deletions(-) diff --git a/examples/PySDM_examples/Arabas_and_Shima_2017/settings.py b/examples/PySDM_examples/Arabas_and_Shima_2017/settings.py index e2ef220425..a4fa31a0e9 100644 --- a/examples/PySDM_examples/Arabas_and_Shima_2017/settings.py +++ b/examples/PySDM_examples/Arabas_and_Shima_2017/settings.py @@ -37,7 +37,7 @@ def __init__( self.n_output = 500 self.rtol_x = condensation.DEFAULTS.rtol_x - self.rtol_thd = condensation.DEFAULTS.rtol_thd + self.rtol_RH = condensation.DEFAULTS.rtol_RH self.coord = "volume logarithm" self.dt_cond_range = condensation.DEFAULTS.cond_range diff --git a/examples/PySDM_examples/Arabas_and_Shima_2017/simulation.py b/examples/PySDM_examples/Arabas_and_Shima_2017/simulation.py index 52f7db4b5e..3e69187862 100644 --- a/examples/PySDM_examples/Arabas_and_Shima_2017/simulation.py +++ b/examples/PySDM_examples/Arabas_and_Shima_2017/simulation.py @@ -34,7 +34,7 @@ def __init__(self, settings, backend=CPU): builder.add_dynamic( Condensation( rtol_x=settings.rtol_x, - rtol_thd=settings.rtol_thd, + rtol_RH=settings.rtol_RH, dt_cond_range=settings.dt_cond_range, ) ) diff --git a/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_5_SCIPY_VS_ADAPTIVE.py b/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_5_SCIPY_VS_ADAPTIVE.py index f87e1cf2df..1436eac3dc 100644 --- a/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_5_SCIPY_VS_ADAPTIVE.py +++ b/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_5_SCIPY_VS_ADAPTIVE.py @@ -34,7 +34,7 @@ def data(n_output, rtols, schemes, setups_num): for settings_idx in range(setups_num): settings = setups[settings_idx] settings.rtol_x = rtol - settings.rtol_thd = rtol + settings.rtol_RH = rtol settings.n_output = n_output simulation = Simulation( settings, backend=CPU if scheme == "CPU" else GPU diff --git a/examples/PySDM_examples/Bartman_et_al_2021/demo.ipynb b/examples/PySDM_examples/Bartman_et_al_2021/demo.ipynb index b5bb39d47b..9f7cf9f1e0 100644 --- a/examples/PySDM_examples/Bartman_et_al_2021/demo.ipynb +++ b/examples/PySDM_examples/Bartman_et_al_2021/demo.ipynb @@ -60,7 +60,7 @@ "settings.simulation_time = .175 * settings.spin_up_time\n", "settings.output_interval = 1 * si.minute\n", "settings.condensation_rtol_x = 1e-6\n", - "settings.condensation_rtol_thd = 5e-7\n", + "settings.condensation_rtol_RH = 5e-7\n", "\n", "settings.condensation_dt_cond_range = (.25*si.s, settings.dt)\n", "settings.coalescence_dt_coal_range = settings.condensation_dt_cond_range\n", diff --git a/examples/PySDM_examples/Bartman_et_al_2021/demo_fig2.ipynb b/examples/PySDM_examples/Bartman_et_al_2021/demo_fig2.ipynb index feefdbd56a..f27f1bbc99 100644 --- a/examples/PySDM_examples/Bartman_et_al_2021/demo_fig2.ipynb +++ b/examples/PySDM_examples/Bartman_et_al_2021/demo_fig2.ipynb @@ -45,10 +45,10 @@ "tol = .5e-6\n", "m = 2\n", "runs = (\n", - " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_thd': tol}},\n", - " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_thd': tol*m}},\n", - " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_thd': tol*m*m}},\n", - " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_thd': tol*m*m*m}},\n", + " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_RH': tol}},\n", + " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_RH': tol*m}},\n", + " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_RH': tol*m*m}},\n", + " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_RH': tol*m*m*m}},\n", " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': False, 'condensation_substeps': 128}},\n", " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': False, 'condensation_substeps': 32}},\n", " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': False, 'condensation_substeps': 8}},\n", @@ -87,7 +87,7 @@ " settings.simulation_time = settings.spin_up_time * (1 if 'CI' not in os.environ else .1)\n", " settings.output_interval = settings.dt\n", " settings.condensation_rtol_x = 1e-6\n", - " settings.condensation_rtol_thd = -1\n", + " settings.condensation_rtol_RH = -1\n", " settings.condensation_schedule = 'dynamic'\n", " settings.kappa = .8\n", " \n", diff --git a/examples/PySDM_examples/Bartman_et_al_2021/demo_fig3.ipynb b/examples/PySDM_examples/Bartman_et_al_2021/demo_fig3.ipynb index 75d155730c..169b9b3be3 100644 --- a/examples/PySDM_examples/Bartman_et_al_2021/demo_fig3.ipynb +++ b/examples/PySDM_examples/Bartman_et_al_2021/demo_fig3.ipynb @@ -88,7 +88,7 @@ " settings.output_interval = settings.dt\n", " settings.condensation_adaptive = True\n", " settings.condensation_rtol_x = 1e-6\n", - " settings.condensation_rtol_thd = 2e-5/7/7\n", + " settings.condensation_rtol_RH = 2e-5/7/7\n", " settings.condensation_schedule = 'dynamic'\n", " settings.kappa = .8\n", " \n", diff --git a/examples/PySDM_examples/Morrison_and_Grabowski_2007/common.py b/examples/PySDM_examples/Morrison_and_Grabowski_2007/common.py index 7cda18e28f..faa8d2040a 100644 --- a/examples/PySDM_examples/Morrison_and_Grabowski_2007/common.py +++ b/examples/PySDM_examples/Morrison_and_Grabowski_2007/common.py @@ -21,7 +21,7 @@ def __init__(self, formulae: Formulae): const = formulae.constants self.condensation_rtol_x = condensation.DEFAULTS.rtol_x - self.condensation_rtol_thd = condensation.DEFAULTS.rtol_thd + self.condensation_rtol_RH = condensation.DEFAULTS.rtol_RH self.condensation_adaptive = True self.condensation_substeps = -1 self.condensation_dt_cond_range = condensation.DEFAULTS.cond_range diff --git a/examples/PySDM_examples/Pyrcel/simulation.py b/examples/PySDM_examples/Pyrcel/simulation.py index 330936cb26..dd966db8e0 100644 --- a/examples/PySDM_examples/Pyrcel/simulation.py +++ b/examples/PySDM_examples/Pyrcel/simulation.py @@ -13,7 +13,7 @@ class Simulation(BasicSimulation): def __init__( - self, settings, products=None, scipy_solver=False, rtol_thd=1e-10, rtol_x=1e-10 + self, settings, products=None, scipy_solver=False, rtol_RH=1e-6, rtol_x=1e-10 ): env = Parcel( dt=settings.timestep, @@ -27,7 +27,7 @@ def __init__( builder = Builder(n_sd=n_sd, backend=CPU(formulae=settings.formulae)) builder.set_environment(env) builder.add_dynamic(AmbientThermodynamics()) - builder.add_dynamic(Condensation(rtol_thd=rtol_thd, rtol_x=rtol_x)) + builder.add_dynamic(Condensation(rtol_RH=rtol_RH, rtol_x=rtol_x)) volume = env.mass_of_dry_air / settings.initial_air_density attributes = { diff --git a/examples/PySDM_examples/Shipway_and_Hill_2012/settings.py b/examples/PySDM_examples/Shipway_and_Hill_2012/settings.py index 5773dcf62e..e6e881290b 100644 --- a/examples/PySDM_examples/Shipway_and_Hill_2012/settings.py +++ b/examples/PySDM_examples/Shipway_and_Hill_2012/settings.py @@ -127,7 +127,7 @@ def drhod_dz(z_above_reservoir, rhod): self.mpdata_settings = {"n_iters": 3, "iga": True, "fct": True, "tot": True} self.condensation_rtol_x = condensation.DEFAULTS.rtol_x - self.condensation_rtol_thd = condensation.DEFAULTS.rtol_thd + self.condensation_rtol_RH = condensation.DEFAULTS.rtol_RH self.condensation_adaptive = True self.condensation_update_thd = False self.coalescence_adaptive = True diff --git a/examples/PySDM_examples/Shipway_and_Hill_2012/simulation.py b/examples/PySDM_examples/Shipway_and_Hill_2012/simulation.py index 1712e202ca..9d0d70c5a7 100644 --- a/examples/PySDM_examples/Shipway_and_Hill_2012/simulation.py +++ b/examples/PySDM_examples/Shipway_and_Hill_2012/simulation.py @@ -69,7 +69,7 @@ def zZ_to_z_above_reservoir(zZ): self.builder.add_dynamic( Condensation( adaptive=settings.condensation_adaptive, - rtol_thd=settings.condensation_rtol_thd, + rtol_RH=settings.condensation_rtol_RH, rtol_x=settings.condensation_rtol_x, update_thd=settings.condensation_update_thd, ) diff --git a/examples/PySDM_examples/Szumowski_et_al_1998/gui_settings.py b/examples/PySDM_examples/Szumowski_et_al_1998/gui_settings.py index 9487692ac4..1140224812 100644 --- a/examples/PySDM_examples/Szumowski_et_al_1998/gui_settings.py +++ b/examples/PySDM_examples/Szumowski_et_al_1998/gui_settings.py @@ -80,16 +80,16 @@ def __init__(self, settings): description="simulation time $[s]$", ) self.ui_condensation_rtol_x = IntSlider( - value=np.log10(settings.condensation_rtol_thd), - min=-9, - max=-3, + value=np.log10(settings.condensation_rtol_RH), + min=-7, + max=-1, description="log$_{10}$(rtol$_x$)", ) - self.ui_condensation_rtol_thd = IntSlider( - value=np.log10(settings.condensation_rtol_thd), + self.ui_condensation_rtol_RH = IntSlider( + value=np.log10(settings.condensation_rtol_RH), min=-9, max=-3, - description="log$_{10}$(rtol$_\\theta$)", + description="log$_{10}$(rtol$_{RH}$)", ) self.ui_condensation_adaptive = Checkbox( value=settings.condensation_adaptive, @@ -274,8 +274,8 @@ def condensation_rtol_x(self): return 10**self.ui_condensation_rtol_x.value @property - def condensation_rtol_thd(self): - return 10**self.ui_condensation_rtol_thd.value + def condensation_rtol_RH(self): + return 10**self.ui_condensation_rtol_RH.value @property def condensation_adaptive(self): @@ -371,7 +371,7 @@ def box(self): self.ui_dt, self.ui_simulation_time, self.ui_condensation_rtol_x, - self.ui_condensation_rtol_thd, + self.ui_condensation_rtol_RH, self.ui_condensation_adaptive, self.ui_coalescence_adaptive, self.ui_displacement_rtol, diff --git a/examples/PySDM_examples/Szumowski_et_al_1998/simulation.py b/examples/PySDM_examples/Szumowski_et_al_1998/simulation.py index 9769147ee5..1602562fab 100644 --- a/examples/PySDM_examples/Szumowski_et_al_1998/simulation.py +++ b/examples/PySDM_examples/Szumowski_et_al_1998/simulation.py @@ -55,7 +55,7 @@ def reinit(self, products=None): if self.settings.processes["condensation"]: condensation = Condensation( rtol_x=self.settings.condensation_rtol_x, - rtol_thd=self.settings.condensation_rtol_thd, + rtol_RH=self.settings.condensation_rtol_RH, adaptive=self.settings.condensation_adaptive, substeps=self.settings.condensation_substeps, dt_cond_range=self.settings.condensation_dt_cond_range, diff --git a/examples/PySDM_examples/Yang_et_al_2018/fig_2.ipynb b/examples/PySDM_examples/Yang_et_al_2018/fig_2.ipynb index c44be0ac50..6d34a1a516 100644 --- a/examples/PySDM_examples/Yang_et_al_2018/fig_2.ipynb +++ b/examples/PySDM_examples/Yang_et_al_2018/fig_2.ipynb @@ -97,7 +97,7 @@ " result['T'] = np.array(output[\"T\"])\n", " result['n'] = settings.n / (settings.mass_of_dry_air * si.kilogram)\n", " result['dt_max'] = settings.dt_max\n", - " result['rtol_thd'] = settings.rtol_thd\n", + " result['rtol_RH'] = settings.rtol_RH\n", " result['rtol_x'] = settings.rtol_x\n", " result['dt_cond_max'] = output['dt_cond_max']\n", " result['dt_cond_min'] = output['dt_cond_min']\n", @@ -134,8 +134,8 @@ "for i, output in enumerate(outputs):\n", " dt_max = output['dt_max']\n", " rtol_x = output['rtol_x']\n", - " rtol_thd = output['rtol_thd']\n", - " tols = f'rtol_x = {rtol_x}, rtol_thd = {rtol_thd}'\n", + " rtol_RH = output['rtol_RH']\n", + " tols = f'rtol_x = {rtol_x}, rtol_RH = {rtol_RH}'\n", " \n", " ax[i,1].set_title(f'dt_max = {dt_max}, '+tols)\n", " \n", @@ -195,8 +195,8 @@ " tols = f'tolerance = {rtol}'\n", " else:\n", " rtol_x = output['rtol_x']\n", - " rtol_thd = output['rtol_thd']\n", - " tols = f'rtol_x = {rtol_x}, rtol_thd = {rtol_thd}'\n", + " rtol_RH = output['rtol_RH']\n", + " tols = f'rtol_x = {rtol_x}, rtol_RH = {rtol_RH}'\n", " \n", " hist = output['r_bins_values']\n", " hist = (hist[:,0:-1] + hist[:,1:])/2\n", diff --git a/examples/PySDM_examples/Yang_et_al_2018/settings.py b/examples/PySDM_examples/Yang_et_al_2018/settings.py index 98d54216a4..72ea0df5d7 100644 --- a/examples/PySDM_examples/Yang_et_al_2018/settings.py +++ b/examples/PySDM_examples/Yang_et_al_2018/settings.py @@ -44,7 +44,7 @@ def __init__( self.coord = "VolumeLogarithm" self.adaptive = True self.rtol_x = condensation.DEFAULTS.rtol_x - self.rtol_thd = condensation.DEFAULTS.rtol_thd + self.rtol_RH = condensation.DEFAULTS.rtol_RH self.dt_cond_range = condensation.DEFAULTS.cond_range self.T0 = 284.3 * si.kelvin diff --git a/examples/PySDM_examples/Yang_et_al_2018/simulation.py b/examples/PySDM_examples/Yang_et_al_2018/simulation.py index 8c9d77bcb2..be7ec1b7b2 100644 --- a/examples/PySDM_examples/Yang_et_al_2018/simulation.py +++ b/examples/PySDM_examples/Yang_et_al_2018/simulation.py @@ -36,7 +36,7 @@ def __init__(self, settings, backend=CPU): condensation = Condensation( adaptive=settings.adaptive, rtol_x=settings.rtol_x, - rtol_thd=settings.rtol_thd, + rtol_RH=settings.rtol_RH, dt_cond_range=settings.dt_cond_range, ) builder.add_dynamic(condensation) diff --git a/tests/devops_tests b/tests/devops_tests index 7eaf1552e0..3641b80ba3 160000 --- a/tests/devops_tests +++ b/tests/devops_tests @@ -1 +1 @@ -Subproject commit 7eaf1552e0667ab7eee0fca18242071df26edca2 +Subproject commit 3641b80ba396d2bb500623aa09138b33b70b5cd7 From 2fd9f5ba9750aa139ed396fa37f7f4cdcb429cc3 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 7 Jun 2023 22:30:16 +0200 Subject: [PATCH 13/23] restore devops tests to version from main --- tests/devops_tests | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/devops_tests b/tests/devops_tests index 3641b80ba3..7eaf1552e0 160000 --- a/tests/devops_tests +++ b/tests/devops_tests @@ -1 +1 @@ -Subproject commit 3641b80ba396d2bb500623aa09138b33b70b5cd7 +Subproject commit 7eaf1552e0667ab7eee0fca18242071df26edca2 From abcde66f92be215b43b0c94f7663261e23f8779f Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 7 Jun 2023 22:42:52 +0200 Subject: [PATCH 14/23] todo labels, pylint --- PySDM/backends/impl_numba/methods/condensation_methods.py | 2 +- PySDM/dynamics/condensation.py | 2 +- .../abdul_razzak_ghan_2000/test_single_supersaturation_peak.py | 1 + 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/condensation_methods.py b/PySDM/backends/impl_numba/methods/condensation_methods.py index 6d65240053..85a8279969 100644 --- a/PySDM/backends/impl_numba/methods/condensation_methods.py +++ b/PySDM/backends/impl_numba/methods/condensation_methods.py @@ -420,7 +420,7 @@ def calculate_ml_new( # pylint: disable=too-many-arguments,too-many-statements, lambdaK = phys_lambdaK(T, p) lambdaD = phys_lambdaD(DTp, T) for drop in cell_idx: - if v[drop] < 0: # TODO: use _unfrozen from freezing_methods.py + if v[drop] < 0: # TODO #1086: use _unfrozen from freezing_methods.py continue x_old = x(v[drop]) r_old = radius(v[drop]) diff --git a/PySDM/dynamics/condensation.py b/PySDM/dynamics/condensation.py index 90069ab91f..e07d738666 100644 --- a/PySDM/dynamics/condensation.py +++ b/PySDM/dynamics/condensation.py @@ -57,7 +57,7 @@ def register(self, builder): adaptive=self.adaptive, fuse=32, multiplier=2, - RH_rtol=1e-7, # TODO: use rtol_RH? + RH_rtol=1e-7, # TODO #868: use rtol_RH? max_iters=self.max_iters, ) builder.request_attribute("critical volume") diff --git a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py index 7f270e31d8..6fa6d8a174 100644 --- a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py +++ b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py @@ -28,6 +28,7 @@ @pytest.mark.parametrize("rtol_x", (1e-7,)) @pytest.mark.parametrize("scheme", ("PySDM",)) def test_single_supersaturation_peak(scheme, rtol_x, rtol_RH, plot=False): + # pylint: disable=too-many-locals # arrange products = ( PySDM_products.WaterMixingRatio(unit="g/kg", name="ql"), From d1f15a2fd66a6a1a43cdbf7996a5e8b02d3c00c2 Mon Sep 17 00:00:00 2001 From: Oleksii Bulenok Date: Fri, 8 Jul 2022 11:59:40 +0200 Subject: [PATCH 15/23] Change thd_tol to RH_tol --- .../methods/condensation_methods.py | 43 +++++++++++-------- .../scipy_ode_condensation_solver.py | 4 +- .../methods/condensation_methods.py | 2 +- PySDM/dynamics/condensation.py | 17 +++++--- PySDM/particulator.py | 5 ++- .../test_single_supersaturation_peak.py | 30 ++++++++----- 6 files changed, 62 insertions(+), 39 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/condensation_methods.py b/PySDM/backends/impl_numba/methods/condensation_methods.py index 67c8e1dba1..96a427a5bb 100644 --- a/PySDM/backends/impl_numba/methods/condensation_methods.py +++ b/PySDM/backends/impl_numba/methods/condensation_methods.py @@ -29,6 +29,7 @@ def condensation( rhod, thd, qv, + RH, dv, prhod, pthd, @@ -36,7 +37,7 @@ def condensation( kappa, f_org, rtol_x, - rtol_thd, + rtol_RH, timestep, counters, cell_order, @@ -58,6 +59,7 @@ def condensation( rhod=rhod.data, thd=thd.data, qv=qv.data, + RH=RH.data, dv_mean=dv, prhod=prhod.data, pthd=pthd.data, @@ -65,7 +67,7 @@ def condensation( kappa=kappa.data, f_org=f_org.data, rtol_x=rtol_x, - rtol_thd=rtol_thd, + rtol_RH=rtol_RH, timestep=timestep, counter_n_substeps=counters["n_substeps"].data, counter_n_activating=counters["n_activating"].data, @@ -92,6 +94,7 @@ def _condensation( rhod, thd, qv, + RH, dv_mean, prhod, pthd, @@ -99,7 +102,7 @@ def _condensation( kappa, f_org, rtol_x, - rtol_thd, + rtol_RH, timestep, counter_n_substeps, counter_n_activating, @@ -145,12 +148,13 @@ def _condensation( thd[cell_id], qv[cell_id], rhod[cell_id], + RH[cell_id], dthd_dt, dqv_dt, drhod_dt, md, rtol_x, - rtol_thd, + rtol_RH, timestep, counter_n_substeps[cell_id], ) @@ -177,7 +181,7 @@ def make_adapt_substeps( n_substeps_min = math.ceil(timestep / dt_range[1]) @numba.njit(**jit_flags) - def adapt_substeps(args, n_substeps, thd, rtol_thd): + def adapt_substeps(args, n_substeps, thd, RH, rtol_RH): n_substeps = np.maximum(n_substeps_min, n_substeps // multiplier) success = False for burnout in range(fuse + 1): @@ -191,24 +195,26 @@ def adapt_substeps(args, n_substeps, thd, rtol_thd): ), return_value=(0, False), ) - thd_new_long, success = step_fake(args, timestep, n_substeps) + thd_new_long, RH_new_long, success = step_fake(args, timestep, n_substeps) if success: break n_substeps *= multiplier for burnout in range(fuse + 1): if burnout == fuse: return warn("burnout (short)", __file__, return_value=(0, False)) - thd_new_short, success = step_fake( + thd_new_short, RH_new_short, success = step_fake( args, timestep, n_substeps * multiplier ) if not success: return warn("short failed", __file__, return_value=(0, False)) - dthd_long = thd_new_long - thd - dthd_short = thd_new_short - thd - error_estimate = np.abs(dthd_long - multiplier * dthd_short) - thd_new_long = thd_new_short - if within_tolerance(error_estimate, thd, rtol_thd): + + dRH_long = RH_new_long - RH + dRH_short = RH_new_short - RH + error_estimate = np.abs(dRH_long - multiplier * dRH_short) + RH_new_long = RH_new_short + if within_tolerance(error_estimate, RH, rtol_RH): break + n_substeps *= multiplier if n_substeps > n_substeps_max: break @@ -221,8 +227,8 @@ def make_step_fake(jit_flags, step_impl): @numba.njit(**jit_flags) def step_fake(args, dt, n_substeps): dt /= n_substeps - _, thd_new, _, _, _, _, success = step_impl(*args, dt, 1, True) - return thd_new, success + _, thd_new, RH_new, _, _, _, _, success = step_impl(*args, dt, 1, True) + return thd_new, RH_new, success return step_fake @@ -333,6 +339,7 @@ def step_impl( # pylint: disable=too-many-arguments,too-many-locals return ( qv, thd, + RH, count_activating, count_deactivating, count_ripening, @@ -417,7 +424,7 @@ def calculate_ml_new( # pylint: disable=too-many-arguments,too-many-statements, lambdaK = phys_lambdaK(T, p) lambdaD = phys_lambdaD(DTp, T) for drop in cell_idx: - if v[drop] < 0: + if v[drop] < 0: # TODO: use _unfrozen from freezing_methods.py continue x_old = x(v[drop]) r_old = radius(v[drop]) @@ -666,12 +673,13 @@ def solve( # pylint: disable=too-many-arguments thd, qv, rhod, + RH, dthd_dt, dqv_dt, drhod_dt, m_d, rtol_x, - rtol_thd, + rtol_RH, timestep, n_substeps, ): @@ -694,11 +702,12 @@ def solve( # pylint: disable=too-many-arguments ) success = True if adaptive: - n_substeps, success = adapt_substeps(args, n_substeps, thd, rtol_thd) + n_substeps, success = adapt_substeps(args, n_substeps, thd, RH, rtol_RH) if success: ( qv, thd, + _, n_activating, n_deactivating, n_ripening, diff --git a/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py b/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py index 5da22468cb..547fb9b9f2 100644 --- a/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py +++ b/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py @@ -24,7 +24,7 @@ def patch_particulator(particulator): def _condensation( - particulator, *, rtol_x, rtol_thd, counters, RH_max, success, cell_order + particulator, *, rtol_x, rtol_RH, counters, RH_max, success, cell_order ): func = Numba._condensation if not numba.config.DISABLE_JIT: # pylint: disable=no-member @@ -49,7 +49,7 @@ def _condensation( kappa=particulator.attributes["kappa"].data, f_org=particulator.attributes["dry volume organic fraction"].data, rtol_x=rtol_x, - rtol_thd=rtol_thd, + rtol_RH=rtol_RH, timestep=particulator.dt, counter_n_substeps=counters["n_substeps"], counter_n_activating=counters["n_activating"], diff --git a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py index 6a49945810..ff8c96525a 100644 --- a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py +++ b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py @@ -339,7 +339,7 @@ def condensation( kappa, f_org, rtol_x, - rtol_thd, + rtol_RH, timestep, counters, cell_order, diff --git a/PySDM/dynamics/condensation.py b/PySDM/dynamics/condensation.py index 2d3dae6460..6631800dc5 100644 --- a/PySDM/dynamics/condensation.py +++ b/PySDM/dynamics/condensation.py @@ -8,9 +8,14 @@ from ..physics import si -DEFAULTS = namedtuple("_", ("rtol_x", "rtol_thd", "cond_range", "schedule"))( +DEFAULTS = namedtuple("_", ( + "rtol_x", + "rtol_RH", + "cond_range", + "schedule" +))( rtol_x=1e-6, - rtol_thd=1e-6, + rtol_RH=1e-4, cond_range=(1e-4 * si.second, 1 * si.second), schedule="dynamic", ) @@ -21,7 +26,7 @@ def __init__( self, *, rtol_x=DEFAULTS.rtol_x, - rtol_thd=DEFAULTS.rtol_thd, + rtol_RH=DEFAULTS.rtol_RH, substeps: int = 1, adaptive: bool = True, dt_cond_range: tuple = DEFAULTS.cond_range, @@ -38,7 +43,7 @@ def __init__( self.enable = True self.rtol_x = rtol_x - self.rtol_thd = rtol_thd + self.rtol_RH = rtol_RH self.rh_max = None self.success = None @@ -62,7 +67,7 @@ def register(self, builder): adaptive=self.adaptive, fuse=32, multiplier=2, - RH_rtol=1e-7, + RH_rtol=1e-7, # TODO: use rtol_RH? max_iters=self.max_iters, ) builder.request_attribute("critical volume") @@ -99,7 +104,7 @@ def __call__(self): self.particulator.condensation( rtol_x=self.rtol_x, - rtol_thd=self.rtol_thd, + rtol_RH=self.rtol_RH, counters=self.counters, RH_max=self.rh_max, success=self.success, diff --git a/PySDM/particulator.py b/PySDM/particulator.py index b7241610d5..521c2104b9 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -101,7 +101,7 @@ def update_TpRH(self): RH=self.environment.get_predicted("RH"), ) - def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_order): + def condensation(self, *, rtol_x, rtol_RH, counters, RH_max, success, cell_order): """Updates droplet volumes by simulating condensation driven by prior changes in environment thermodynamic state, updates the environment state. In the case of parcel environment, condensation is driven solely by changes in @@ -121,6 +121,7 @@ def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_orde rhod=self.environment["rhod"], thd=self.environment["thd"], qv=self.environment["qv"], + RH=self.environment["RH"], dv=self.environment.dv, prhod=self.environment.get_predicted("rhod"), pthd=self.environment.get_predicted("thd"), @@ -128,7 +129,7 @@ def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_orde kappa=self.attributes["kappa"], f_org=self.attributes["dry volume organic fraction"], rtol_x=rtol_x, - rtol_thd=rtol_thd, + rtol_RH=rtol_RH, v_cr=self.attributes["critical volume"], timestep=self.dt, counters=counters, diff --git a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py index f31073ccfa..19355e57e5 100644 --- a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py +++ b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py @@ -17,7 +17,7 @@ @pytest.mark.parametrize( - "rtol_thd", + "rtol_RH", ( pytest.param(1e-6, marks=pytest.mark.xfail(strict=True)), 1e-7, @@ -26,11 +26,13 @@ ), ) @pytest.mark.parametrize("rtol_x", (1e-7,)) -@pytest.mark.parametrize("adaptive", (True,)) @pytest.mark.parametrize("scheme", ("PySDM",)) def test_single_supersaturation_peak( - adaptive, scheme, rtol_x, rtol_thd, plot=False -): # pylint: disable=too-many-locals + scheme, + rtol_x, + rtol_RH, + plot=False +): # arrange products = ( PySDM_products.WaterMixingRatio(unit="g/kg", name="ql"), @@ -38,15 +40,18 @@ def test_single_supersaturation_peak( PySDM_products.AmbientRelativeHumidity(name="RH"), PySDM_products.ParcelDisplacement(name="z"), ) + dt = 2 * si.s + w = 0.5 * si.m / si.s env = Parcel( - dt=2 * si.s, + dt=dt, mass_of_dry_air=1e3 * si.kg, p0=1000 * si.hPa, q0=22.76 * si.g / si.kg, - w=0.5 * si.m / si.s, + w=w, T0=300 * si.K, ) - n_steps = 70 + z_max = 70 * si.m + n_steps = int(z_max / (w * dt)) n_sd = 2 kappa = 0.4 spectrum = Lognormal(norm_factor=5000 / si.cm**3, m_mode=50.0 * si.nm, s_geom=2.0) @@ -55,9 +60,9 @@ def test_single_supersaturation_peak( builder.add_dynamic(AmbientThermodynamics()) builder.add_dynamic( Condensation( - adaptive=adaptive, + adaptive=True, rtol_x=rtol_x, - rtol_thd=rtol_thd, + rtol_RH=rtol_RH ) ) @@ -110,9 +115,12 @@ def test_single_supersaturation_peak( twin.set_xlim(-0.001, 0.0015) pyplot.legend(loc="lower right") pyplot.grid() - pyplot.title(f"rtol_thd={rtol_thd}; rtol_x={rtol_x}") + pyplot.title(f"rtol_RH={rtol_RH}; rtol_x={rtol_x}") if plot: pyplot.show() # assert - assert signal.argrelextrema(np.asarray(output["RH"]), np.greater)[0].shape[0] == 1 + n_max = signal.argrelextrema(np.asarray(output["RH"]), np.greater)[0].shape[0] + n_min = signal.argrelextrema(np.asarray(output["RH"]), np.less)[0].shape[0] + assert n_max == 1 and n_min == 0 + From 8f04c691873cb72470754fdbf740290cfaa00146 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sat, 9 Jul 2022 00:15:55 +0200 Subject: [PATCH 16/23] making pylint & black happy; bumping test-time-req for examples --- .../impl_numba/methods/condensation_methods.py | 4 ++-- PySDM/dynamics/condensation.py | 7 +------ .../test_single_supersaturation_peak.py | 16 ++-------------- 3 files changed, 5 insertions(+), 22 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/condensation_methods.py b/PySDM/backends/impl_numba/methods/condensation_methods.py index 96a427a5bb..c997ff167c 100644 --- a/PySDM/backends/impl_numba/methods/condensation_methods.py +++ b/PySDM/backends/impl_numba/methods/condensation_methods.py @@ -195,14 +195,14 @@ def adapt_substeps(args, n_substeps, thd, RH, rtol_RH): ), return_value=(0, False), ) - thd_new_long, RH_new_long, success = step_fake(args, timestep, n_substeps) + _, RH_new_long, success = step_fake(args, timestep, n_substeps) if success: break n_substeps *= multiplier for burnout in range(fuse + 1): if burnout == fuse: return warn("burnout (short)", __file__, return_value=(0, False)) - thd_new_short, RH_new_short, success = step_fake( + _, RH_new_short, success = step_fake( args, timestep, n_substeps * multiplier ) if not success: diff --git a/PySDM/dynamics/condensation.py b/PySDM/dynamics/condensation.py index 6631800dc5..e10973a968 100644 --- a/PySDM/dynamics/condensation.py +++ b/PySDM/dynamics/condensation.py @@ -8,12 +8,7 @@ from ..physics import si -DEFAULTS = namedtuple("_", ( - "rtol_x", - "rtol_RH", - "cond_range", - "schedule" -))( +DEFAULTS = namedtuple("_", ("rtol_x", "rtol_RH", "cond_range", "schedule"))( rtol_x=1e-6, rtol_RH=1e-4, cond_range=(1e-4 * si.second, 1 * si.second), diff --git a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py index 19355e57e5..a83f348b51 100644 --- a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py +++ b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py @@ -27,12 +27,7 @@ ) @pytest.mark.parametrize("rtol_x", (1e-7,)) @pytest.mark.parametrize("scheme", ("PySDM",)) -def test_single_supersaturation_peak( - scheme, - rtol_x, - rtol_RH, - plot=False -): +def test_single_supersaturation_peak(scheme, rtol_x, rtol_RH, plot=False): # arrange products = ( PySDM_products.WaterMixingRatio(unit="g/kg", name="ql"), @@ -58,13 +53,7 @@ def test_single_supersaturation_peak( builder = Builder(backend=CPU(), n_sd=n_sd) builder.set_environment(env) builder.add_dynamic(AmbientThermodynamics()) - builder.add_dynamic( - Condensation( - adaptive=True, - rtol_x=rtol_x, - rtol_RH=rtol_RH - ) - ) + builder.add_dynamic(Condensation(adaptive=True, rtol_x=rtol_x, rtol_RH=rtol_RH)) r_dry, concentration = ConstantMultiplicity(spectrum).sample(n_sd) v_dry = builder.formulae.trivia.volume(radius=r_dry) @@ -123,4 +112,3 @@ def test_single_supersaturation_peak( n_max = signal.argrelextrema(np.asarray(output["RH"]), np.greater)[0].shape[0] n_min = signal.argrelextrema(np.asarray(output["RH"]), np.less)[0].shape[0] assert n_max == 1 and n_min == 0 - From 6c90dff08448a3973e4377688bf85b5a0bb0e568 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sat, 9 Jul 2022 00:34:11 +0200 Subject: [PATCH 17/23] fix GPU backend condensation signature --- PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py | 1 + 1 file changed, 1 insertion(+) diff --git a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py index ff8c96525a..951421c624 100644 --- a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py +++ b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py @@ -332,6 +332,7 @@ def condensation( rhod, thd, qv, + RH, dv, prhod, pthd, From aa0b0b4176d47cffa9f0c9683ac368b63a1500e7 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sat, 9 Jul 2022 01:08:04 +0200 Subject: [PATCH 18/23] fix BDF solver API --- .../test_helpers/scipy_ode_condensation_solver.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py b/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py index 547fb9b9f2..13437082ef 100644 --- a/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py +++ b/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py @@ -42,6 +42,7 @@ def _condensation( rhod=particulator.environment["rhod"].data, thd=particulator.environment["thd"].data, qv=particulator.environment["qv"].data, + RH=particulator.environment["RH"].data, dv_mean=particulator.environment.dv, prhod=particulator.environment.get_predicted("rhod").data, pthd=particulator.environment.get_predicted("thd").data, @@ -178,14 +179,15 @@ def solve( # pylint: disable=too-many-arguments,too-many-locals thd, qv, rhod, + __, dthd_dt, dqv_dt, drhod_dt, m_d_mean, - __, ___, - dt, ____, + dt, + _____, ): n_sd_in_cell = len(cell_idx) y0 = np.empty(n_sd_in_cell + idx_x) From bcb9b61875529c05cf89920d12bc1a855acf9401 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sat, 9 Jul 2022 21:00:27 +0200 Subject: [PATCH 19/23] default rtol_RH=1e-6; bump examples req; reduce significant digit number in arabas_and_shima_2017/test_conservation --- PySDM/dynamics/condensation.py | 2 +- .../parcel/arabas_and_shima_2017/test_conservation.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/PySDM/dynamics/condensation.py b/PySDM/dynamics/condensation.py index e10973a968..f80af81bc6 100644 --- a/PySDM/dynamics/condensation.py +++ b/PySDM/dynamics/condensation.py @@ -10,7 +10,7 @@ DEFAULTS = namedtuple("_", ("rtol_x", "rtol_RH", "cond_range", "schedule"))( rtol_x=1e-6, - rtol_RH=1e-4, + rtol_RH=1e-6, cond_range=(1e-4 * si.second, 1 * si.second), schedule="dynamic", ) diff --git a/tests/smoke_tests/parcel/arabas_and_shima_2017/test_conservation.py b/tests/smoke_tests/parcel/arabas_and_shima_2017/test_conservation.py index 462b1f9a8c..e5d40d4254 100644 --- a/tests/smoke_tests/parcel/arabas_and_shima_2017/test_conservation.py +++ b/tests/smoke_tests/parcel/arabas_and_shima_2017/test_conservation.py @@ -47,7 +47,7 @@ def test_water_mass_conservation(settings_idx, mass_of_dry_air, scheme, coord): # Assert ql = liquid_water_mixing_ratio(simulation) qt = simulation.particulator.environment["qv"].to_ndarray() + ql - significant = 6 if scheme == "GPU" else 14 # TODO #540 + significant = 6 if scheme == "GPU" else 12 # TODO #540 np.testing.assert_approx_equal(qt, qt0, significant) if scheme != "SciPy": assert simulation.particulator.products["S_max"].get() >= output["S"][-1] From 6a300bc9047a1f4a469afc08fb4582f9577e91f8 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 7 Jun 2023 22:20:27 +0200 Subject: [PATCH 20/23] updating examples code --- .../Arabas_and_Shima_2017/settings.py | 2 +- .../Arabas_and_Shima_2017/simulation.py | 2 +- .../fig_5_SCIPY_VS_ADAPTIVE.py | 2 +- .../Bartman_et_al_2021/demo.ipynb | 2 +- .../Bartman_et_al_2021/demo_fig2.ipynb | 10 +++++----- .../Bartman_et_al_2021/demo_fig3.ipynb | 2 +- .../Morrison_and_Grabowski_2007/common.py | 2 +- examples/PySDM_examples/Pyrcel/simulation.py | 4 ++-- .../Shipway_and_Hill_2012/settings.py | 2 +- .../Shipway_and_Hill_2012/simulation.py | 2 +- .../Szumowski_et_al_1998/gui_settings.py | 18 +++++++++--------- .../Szumowski_et_al_1998/simulation.py | 2 +- .../PySDM_examples/Yang_et_al_2018/fig_2.ipynb | 10 +++++----- .../PySDM_examples/Yang_et_al_2018/settings.py | 2 +- .../Yang_et_al_2018/simulation.py | 2 +- 15 files changed, 32 insertions(+), 32 deletions(-) diff --git a/examples/PySDM_examples/Arabas_and_Shima_2017/settings.py b/examples/PySDM_examples/Arabas_and_Shima_2017/settings.py index e2ef220425..a4fa31a0e9 100644 --- a/examples/PySDM_examples/Arabas_and_Shima_2017/settings.py +++ b/examples/PySDM_examples/Arabas_and_Shima_2017/settings.py @@ -37,7 +37,7 @@ def __init__( self.n_output = 500 self.rtol_x = condensation.DEFAULTS.rtol_x - self.rtol_thd = condensation.DEFAULTS.rtol_thd + self.rtol_RH = condensation.DEFAULTS.rtol_RH self.coord = "volume logarithm" self.dt_cond_range = condensation.DEFAULTS.cond_range diff --git a/examples/PySDM_examples/Arabas_and_Shima_2017/simulation.py b/examples/PySDM_examples/Arabas_and_Shima_2017/simulation.py index 52f7db4b5e..3e69187862 100644 --- a/examples/PySDM_examples/Arabas_and_Shima_2017/simulation.py +++ b/examples/PySDM_examples/Arabas_and_Shima_2017/simulation.py @@ -34,7 +34,7 @@ def __init__(self, settings, backend=CPU): builder.add_dynamic( Condensation( rtol_x=settings.rtol_x, - rtol_thd=settings.rtol_thd, + rtol_RH=settings.rtol_RH, dt_cond_range=settings.dt_cond_range, ) ) diff --git a/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_5_SCIPY_VS_ADAPTIVE.py b/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_5_SCIPY_VS_ADAPTIVE.py index f87e1cf2df..1436eac3dc 100644 --- a/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_5_SCIPY_VS_ADAPTIVE.py +++ b/examples/PySDM_examples/Bartman_2020_MasterThesis/fig_5_SCIPY_VS_ADAPTIVE.py @@ -34,7 +34,7 @@ def data(n_output, rtols, schemes, setups_num): for settings_idx in range(setups_num): settings = setups[settings_idx] settings.rtol_x = rtol - settings.rtol_thd = rtol + settings.rtol_RH = rtol settings.n_output = n_output simulation = Simulation( settings, backend=CPU if scheme == "CPU" else GPU diff --git a/examples/PySDM_examples/Bartman_et_al_2021/demo.ipynb b/examples/PySDM_examples/Bartman_et_al_2021/demo.ipynb index b5bb39d47b..9f7cf9f1e0 100644 --- a/examples/PySDM_examples/Bartman_et_al_2021/demo.ipynb +++ b/examples/PySDM_examples/Bartman_et_al_2021/demo.ipynb @@ -60,7 +60,7 @@ "settings.simulation_time = .175 * settings.spin_up_time\n", "settings.output_interval = 1 * si.minute\n", "settings.condensation_rtol_x = 1e-6\n", - "settings.condensation_rtol_thd = 5e-7\n", + "settings.condensation_rtol_RH = 5e-7\n", "\n", "settings.condensation_dt_cond_range = (.25*si.s, settings.dt)\n", "settings.coalescence_dt_coal_range = settings.condensation_dt_cond_range\n", diff --git a/examples/PySDM_examples/Bartman_et_al_2021/demo_fig2.ipynb b/examples/PySDM_examples/Bartman_et_al_2021/demo_fig2.ipynb index feefdbd56a..f27f1bbc99 100644 --- a/examples/PySDM_examples/Bartman_et_al_2021/demo_fig2.ipynb +++ b/examples/PySDM_examples/Bartman_et_al_2021/demo_fig2.ipynb @@ -45,10 +45,10 @@ "tol = .5e-6\n", "m = 2\n", "runs = (\n", - " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_thd': tol}},\n", - " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_thd': tol*m}},\n", - " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_thd': tol*m*m}},\n", - " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_thd': tol*m*m*m}},\n", + " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_RH': tol}},\n", + " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_RH': tol*m}},\n", + " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_RH': tol*m*m}},\n", + " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_RH': tol*m*m*m}},\n", " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': False, 'condensation_substeps': 128}},\n", " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': False, 'condensation_substeps': 32}},\n", " {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': False, 'condensation_substeps': 8}},\n", @@ -87,7 +87,7 @@ " settings.simulation_time = settings.spin_up_time * (1 if 'CI' not in os.environ else .1)\n", " settings.output_interval = settings.dt\n", " settings.condensation_rtol_x = 1e-6\n", - " settings.condensation_rtol_thd = -1\n", + " settings.condensation_rtol_RH = -1\n", " settings.condensation_schedule = 'dynamic'\n", " settings.kappa = .8\n", " \n", diff --git a/examples/PySDM_examples/Bartman_et_al_2021/demo_fig3.ipynb b/examples/PySDM_examples/Bartman_et_al_2021/demo_fig3.ipynb index 75d155730c..169b9b3be3 100644 --- a/examples/PySDM_examples/Bartman_et_al_2021/demo_fig3.ipynb +++ b/examples/PySDM_examples/Bartman_et_al_2021/demo_fig3.ipynb @@ -88,7 +88,7 @@ " settings.output_interval = settings.dt\n", " settings.condensation_adaptive = True\n", " settings.condensation_rtol_x = 1e-6\n", - " settings.condensation_rtol_thd = 2e-5/7/7\n", + " settings.condensation_rtol_RH = 2e-5/7/7\n", " settings.condensation_schedule = 'dynamic'\n", " settings.kappa = .8\n", " \n", diff --git a/examples/PySDM_examples/Morrison_and_Grabowski_2007/common.py b/examples/PySDM_examples/Morrison_and_Grabowski_2007/common.py index 7cda18e28f..faa8d2040a 100644 --- a/examples/PySDM_examples/Morrison_and_Grabowski_2007/common.py +++ b/examples/PySDM_examples/Morrison_and_Grabowski_2007/common.py @@ -21,7 +21,7 @@ def __init__(self, formulae: Formulae): const = formulae.constants self.condensation_rtol_x = condensation.DEFAULTS.rtol_x - self.condensation_rtol_thd = condensation.DEFAULTS.rtol_thd + self.condensation_rtol_RH = condensation.DEFAULTS.rtol_RH self.condensation_adaptive = True self.condensation_substeps = -1 self.condensation_dt_cond_range = condensation.DEFAULTS.cond_range diff --git a/examples/PySDM_examples/Pyrcel/simulation.py b/examples/PySDM_examples/Pyrcel/simulation.py index 330936cb26..dd966db8e0 100644 --- a/examples/PySDM_examples/Pyrcel/simulation.py +++ b/examples/PySDM_examples/Pyrcel/simulation.py @@ -13,7 +13,7 @@ class Simulation(BasicSimulation): def __init__( - self, settings, products=None, scipy_solver=False, rtol_thd=1e-10, rtol_x=1e-10 + self, settings, products=None, scipy_solver=False, rtol_RH=1e-6, rtol_x=1e-10 ): env = Parcel( dt=settings.timestep, @@ -27,7 +27,7 @@ def __init__( builder = Builder(n_sd=n_sd, backend=CPU(formulae=settings.formulae)) builder.set_environment(env) builder.add_dynamic(AmbientThermodynamics()) - builder.add_dynamic(Condensation(rtol_thd=rtol_thd, rtol_x=rtol_x)) + builder.add_dynamic(Condensation(rtol_RH=rtol_RH, rtol_x=rtol_x)) volume = env.mass_of_dry_air / settings.initial_air_density attributes = { diff --git a/examples/PySDM_examples/Shipway_and_Hill_2012/settings.py b/examples/PySDM_examples/Shipway_and_Hill_2012/settings.py index 5773dcf62e..e6e881290b 100644 --- a/examples/PySDM_examples/Shipway_and_Hill_2012/settings.py +++ b/examples/PySDM_examples/Shipway_and_Hill_2012/settings.py @@ -127,7 +127,7 @@ def drhod_dz(z_above_reservoir, rhod): self.mpdata_settings = {"n_iters": 3, "iga": True, "fct": True, "tot": True} self.condensation_rtol_x = condensation.DEFAULTS.rtol_x - self.condensation_rtol_thd = condensation.DEFAULTS.rtol_thd + self.condensation_rtol_RH = condensation.DEFAULTS.rtol_RH self.condensation_adaptive = True self.condensation_update_thd = False self.coalescence_adaptive = True diff --git a/examples/PySDM_examples/Shipway_and_Hill_2012/simulation.py b/examples/PySDM_examples/Shipway_and_Hill_2012/simulation.py index 1712e202ca..9d0d70c5a7 100644 --- a/examples/PySDM_examples/Shipway_and_Hill_2012/simulation.py +++ b/examples/PySDM_examples/Shipway_and_Hill_2012/simulation.py @@ -69,7 +69,7 @@ def zZ_to_z_above_reservoir(zZ): self.builder.add_dynamic( Condensation( adaptive=settings.condensation_adaptive, - rtol_thd=settings.condensation_rtol_thd, + rtol_RH=settings.condensation_rtol_RH, rtol_x=settings.condensation_rtol_x, update_thd=settings.condensation_update_thd, ) diff --git a/examples/PySDM_examples/Szumowski_et_al_1998/gui_settings.py b/examples/PySDM_examples/Szumowski_et_al_1998/gui_settings.py index 9487692ac4..1140224812 100644 --- a/examples/PySDM_examples/Szumowski_et_al_1998/gui_settings.py +++ b/examples/PySDM_examples/Szumowski_et_al_1998/gui_settings.py @@ -80,16 +80,16 @@ def __init__(self, settings): description="simulation time $[s]$", ) self.ui_condensation_rtol_x = IntSlider( - value=np.log10(settings.condensation_rtol_thd), - min=-9, - max=-3, + value=np.log10(settings.condensation_rtol_RH), + min=-7, + max=-1, description="log$_{10}$(rtol$_x$)", ) - self.ui_condensation_rtol_thd = IntSlider( - value=np.log10(settings.condensation_rtol_thd), + self.ui_condensation_rtol_RH = IntSlider( + value=np.log10(settings.condensation_rtol_RH), min=-9, max=-3, - description="log$_{10}$(rtol$_\\theta$)", + description="log$_{10}$(rtol$_{RH}$)", ) self.ui_condensation_adaptive = Checkbox( value=settings.condensation_adaptive, @@ -274,8 +274,8 @@ def condensation_rtol_x(self): return 10**self.ui_condensation_rtol_x.value @property - def condensation_rtol_thd(self): - return 10**self.ui_condensation_rtol_thd.value + def condensation_rtol_RH(self): + return 10**self.ui_condensation_rtol_RH.value @property def condensation_adaptive(self): @@ -371,7 +371,7 @@ def box(self): self.ui_dt, self.ui_simulation_time, self.ui_condensation_rtol_x, - self.ui_condensation_rtol_thd, + self.ui_condensation_rtol_RH, self.ui_condensation_adaptive, self.ui_coalescence_adaptive, self.ui_displacement_rtol, diff --git a/examples/PySDM_examples/Szumowski_et_al_1998/simulation.py b/examples/PySDM_examples/Szumowski_et_al_1998/simulation.py index 9596abef9c..68a0b1c6b7 100644 --- a/examples/PySDM_examples/Szumowski_et_al_1998/simulation.py +++ b/examples/PySDM_examples/Szumowski_et_al_1998/simulation.py @@ -58,7 +58,7 @@ def reinit(self, products=None): kwargs["substeps"] = (self.settings.condensation_substeps,) condensation = Condensation( rtol_x=self.settings.condensation_rtol_x, - rtol_thd=self.settings.condensation_rtol_thd, + rtol_RH=self.settings.condensation_rtol_RH, adaptive=self.settings.condensation_adaptive, dt_cond_range=self.settings.condensation_dt_cond_range, schedule=self.settings.condensation_schedule, diff --git a/examples/PySDM_examples/Yang_et_al_2018/fig_2.ipynb b/examples/PySDM_examples/Yang_et_al_2018/fig_2.ipynb index c44be0ac50..6d34a1a516 100644 --- a/examples/PySDM_examples/Yang_et_al_2018/fig_2.ipynb +++ b/examples/PySDM_examples/Yang_et_al_2018/fig_2.ipynb @@ -97,7 +97,7 @@ " result['T'] = np.array(output[\"T\"])\n", " result['n'] = settings.n / (settings.mass_of_dry_air * si.kilogram)\n", " result['dt_max'] = settings.dt_max\n", - " result['rtol_thd'] = settings.rtol_thd\n", + " result['rtol_RH'] = settings.rtol_RH\n", " result['rtol_x'] = settings.rtol_x\n", " result['dt_cond_max'] = output['dt_cond_max']\n", " result['dt_cond_min'] = output['dt_cond_min']\n", @@ -134,8 +134,8 @@ "for i, output in enumerate(outputs):\n", " dt_max = output['dt_max']\n", " rtol_x = output['rtol_x']\n", - " rtol_thd = output['rtol_thd']\n", - " tols = f'rtol_x = {rtol_x}, rtol_thd = {rtol_thd}'\n", + " rtol_RH = output['rtol_RH']\n", + " tols = f'rtol_x = {rtol_x}, rtol_RH = {rtol_RH}'\n", " \n", " ax[i,1].set_title(f'dt_max = {dt_max}, '+tols)\n", " \n", @@ -195,8 +195,8 @@ " tols = f'tolerance = {rtol}'\n", " else:\n", " rtol_x = output['rtol_x']\n", - " rtol_thd = output['rtol_thd']\n", - " tols = f'rtol_x = {rtol_x}, rtol_thd = {rtol_thd}'\n", + " rtol_RH = output['rtol_RH']\n", + " tols = f'rtol_x = {rtol_x}, rtol_RH = {rtol_RH}'\n", " \n", " hist = output['r_bins_values']\n", " hist = (hist[:,0:-1] + hist[:,1:])/2\n", diff --git a/examples/PySDM_examples/Yang_et_al_2018/settings.py b/examples/PySDM_examples/Yang_et_al_2018/settings.py index 98d54216a4..72ea0df5d7 100644 --- a/examples/PySDM_examples/Yang_et_al_2018/settings.py +++ b/examples/PySDM_examples/Yang_et_al_2018/settings.py @@ -44,7 +44,7 @@ def __init__( self.coord = "VolumeLogarithm" self.adaptive = True self.rtol_x = condensation.DEFAULTS.rtol_x - self.rtol_thd = condensation.DEFAULTS.rtol_thd + self.rtol_RH = condensation.DEFAULTS.rtol_RH self.dt_cond_range = condensation.DEFAULTS.cond_range self.T0 = 284.3 * si.kelvin diff --git a/examples/PySDM_examples/Yang_et_al_2018/simulation.py b/examples/PySDM_examples/Yang_et_al_2018/simulation.py index 8c9d77bcb2..be7ec1b7b2 100644 --- a/examples/PySDM_examples/Yang_et_al_2018/simulation.py +++ b/examples/PySDM_examples/Yang_et_al_2018/simulation.py @@ -36,7 +36,7 @@ def __init__(self, settings, backend=CPU): condensation = Condensation( adaptive=settings.adaptive, rtol_x=settings.rtol_x, - rtol_thd=settings.rtol_thd, + rtol_RH=settings.rtol_RH, dt_cond_range=settings.dt_cond_range, ) builder.add_dynamic(condensation) From 4634f3e2ffe1ef4ed88da0664219e89769d36c68 Mon Sep 17 00:00:00 2001 From: Oleksii Bulenok Date: Fri, 8 Jul 2022 11:59:40 +0200 Subject: [PATCH 21/23] Change thd_tol to RH_tol --- .../test_single_supersaturation_peak.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py index a83f348b51..38c92d0802 100644 --- a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py +++ b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py @@ -27,7 +27,12 @@ ) @pytest.mark.parametrize("rtol_x", (1e-7,)) @pytest.mark.parametrize("scheme", ("PySDM",)) -def test_single_supersaturation_peak(scheme, rtol_x, rtol_RH, plot=False): +def test_single_supersaturation_peak( + scheme, + rtol_x, + rtol_RH, + plot=False +): # arrange products = ( PySDM_products.WaterMixingRatio(unit="g/kg", name="ql"), @@ -53,7 +58,13 @@ def test_single_supersaturation_peak(scheme, rtol_x, rtol_RH, plot=False): builder = Builder(backend=CPU(), n_sd=n_sd) builder.set_environment(env) builder.add_dynamic(AmbientThermodynamics()) - builder.add_dynamic(Condensation(adaptive=True, rtol_x=rtol_x, rtol_RH=rtol_RH)) + builder.add_dynamic( + Condensation( + adaptive=True, + rtol_x=rtol_x, + rtol_RH=rtol_RH + ) + ) r_dry, concentration = ConstantMultiplicity(spectrum).sample(n_sd) v_dry = builder.formulae.trivia.volume(radius=r_dry) From b20dbc4da1c711debc9216870fe08199e1529e8e Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sat, 9 Jul 2022 00:15:55 +0200 Subject: [PATCH 22/23] making pylint & black happy; bumping test-time-req for examples --- .../test_single_supersaturation_peak.py | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py index 38c92d0802..a83f348b51 100644 --- a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py +++ b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py @@ -27,12 +27,7 @@ ) @pytest.mark.parametrize("rtol_x", (1e-7,)) @pytest.mark.parametrize("scheme", ("PySDM",)) -def test_single_supersaturation_peak( - scheme, - rtol_x, - rtol_RH, - plot=False -): +def test_single_supersaturation_peak(scheme, rtol_x, rtol_RH, plot=False): # arrange products = ( PySDM_products.WaterMixingRatio(unit="g/kg", name="ql"), @@ -58,13 +53,7 @@ def test_single_supersaturation_peak( builder = Builder(backend=CPU(), n_sd=n_sd) builder.set_environment(env) builder.add_dynamic(AmbientThermodynamics()) - builder.add_dynamic( - Condensation( - adaptive=True, - rtol_x=rtol_x, - rtol_RH=rtol_RH - ) - ) + builder.add_dynamic(Condensation(adaptive=True, rtol_x=rtol_x, rtol_RH=rtol_RH)) r_dry, concentration = ConstantMultiplicity(spectrum).sample(n_sd) v_dry = builder.formulae.trivia.volume(radius=r_dry) From e442ee8143eb05cc14eb22a076c663f7010e962e Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 7 Jun 2023 22:42:52 +0200 Subject: [PATCH 23/23] todo labels, pylint --- PySDM/backends/impl_numba/methods/condensation_methods.py | 2 +- PySDM/dynamics/condensation.py | 2 +- .../abdul_razzak_ghan_2000/test_single_supersaturation_peak.py | 1 + 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/PySDM/backends/impl_numba/methods/condensation_methods.py b/PySDM/backends/impl_numba/methods/condensation_methods.py index c997ff167c..941a774d14 100644 --- a/PySDM/backends/impl_numba/methods/condensation_methods.py +++ b/PySDM/backends/impl_numba/methods/condensation_methods.py @@ -424,7 +424,7 @@ def calculate_ml_new( # pylint: disable=too-many-arguments,too-many-statements, lambdaK = phys_lambdaK(T, p) lambdaD = phys_lambdaD(DTp, T) for drop in cell_idx: - if v[drop] < 0: # TODO: use _unfrozen from freezing_methods.py + if v[drop] < 0: # TODO #1086: use _unfrozen from freezing_methods.py continue x_old = x(v[drop]) r_old = radius(v[drop]) diff --git a/PySDM/dynamics/condensation.py b/PySDM/dynamics/condensation.py index f80af81bc6..406cdcbc1d 100644 --- a/PySDM/dynamics/condensation.py +++ b/PySDM/dynamics/condensation.py @@ -62,7 +62,7 @@ def register(self, builder): adaptive=self.adaptive, fuse=32, multiplier=2, - RH_rtol=1e-7, # TODO: use rtol_RH? + RH_rtol=1e-7, # TODO #868: use rtol_RH? max_iters=self.max_iters, ) builder.request_attribute("critical volume") diff --git a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py index a83f348b51..8c5b5c0999 100644 --- a/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py +++ b/tests/smoke_tests/parcel/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py @@ -28,6 +28,7 @@ @pytest.mark.parametrize("rtol_x", (1e-7,)) @pytest.mark.parametrize("scheme", ("PySDM",)) def test_single_supersaturation_peak(scheme, rtol_x, rtol_RH, plot=False): + # pylint: disable=too-many-locals # arrange products = ( PySDM_products.WaterMixingRatio(unit="g/kg", name="ql"),