Added example based on Loftus and Wordsworth 2021 figure 2 #1648

lursz · 2025-06-10T14:44:56Z

Added the notebook and utils for recreating figure 2 from Loftus and Wordsworth 2021.
The changes:

Implemented a comprehensive simulation setup with support for multi-core processing and custom settings,
Introduced the AlienParcel class representing the different experimental environments.
Added tests for accuracy and consistency for the example.

slayoo · 2025-06-10T16:08:07Z

@kaitlyn-loftus, FYI - new example being developed for PySDM - kudos @lursz, @emmacware, @piotrkubala, @Hevagog & @mzsuetam!

kaitlyn-loftus · 2025-06-10T20:10:55Z

great news! happy to share any code or try to help out!!

there are also a bunch of functions for properties of a CO2 atmosphere here: https://github.com/kaitlyn-loftus/serious-cirrus/blob/main/src/%CE%9C%CE%A6.jl

kubijaku · 2025-06-16T09:41:24Z

examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py

+class Earth(Planet):
+    g_std: float = 9.82 * si.metre / si.second**2
+    T_STP: float = 290 * si.kelvin
+    p_STP: float = 1.01325 * 1e6 * si.pascal


do we need to hardcode all constants? cannot we used one set in library with better approximation f.e. scipy.constants?

scipy doesn't have mars, saturn, etc. constants, so for single convention and for readability sake I'd argue it is not needed. I see no value added when we have to specifically import scipy only in the Earth planet and not for others

hi all, guess I am subscribed to email notifications on this thread.

I think p_STP should be

p_STP: float = 1.01325 * 1e5 * si.pascal

not

p_STP: float = 1.01325 * 1e6 * si.pascal

hi all, guess I am subscribed to email notifications on this thread.

Github notifies anyone who has been @-mentioned in a PR

I think p_STP should be

p_STP: float = 1.01325 * 1e5 * si.pascal

not

p_STP: float = 1.01325 * 1e6 * si.pascal

thanks!
we can also use si.hPa

Good catch, we'll fix it right away

As @lursz's right away turns out to take some time, i fixed it here

kubijaku · 2025-06-16T09:49:17Z

examples/PySDM_examples/Loftus_and_Wordsworth_2021/planet.py

+    T_STP: float = 300 * si.kelvin
+    p_STP: float = 1.01325 * 1e6 * si.pascal
+    RH_zref: float = 0.75
+    dry_molar_conc_H2: float = 0


add some comments to clarify what are those variables and why they have such a values (they are just constants or you toke them from another paper)?

Thanks for the proposition, fixed here.

makowskitomasz

Good job

makowskitomasz · 2025-06-17T08:13:21Z

tests/examples_tests/Loftus_and_Wordsworth_2021/unit_test.py

Current tests cover only valid input scenarios. To improve robustness, consider adding edge case tests, such as:

r_wet < 0

mass_of_dry_air < 0

Extremely low/high values for Tcloud or pcloud

with pytest.raises(ValueError): Settings( ... )

makowskitomasz · 2025-06-17T08:24:13Z

examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py

+    def advance_parcel_vars(self):
+        """
+        Compute new values of displacement, dry-air density and volume,
+        and write them to self._tmp and self.mesh.dv
+        """
+        dt = self.particulator.dt
+        formulae = self.particulator.formulae
+        T = self["T"][0]
+        p = self["p"][0]
+
+        dz_dt = -self.particulator.attributes["terminal velocity"].to_ndarray()[0]
+        water_vapour_mixing_ratio = (
+            self["water_vapour_mixing_ratio"][0]
+            - self.delta_liquid_water_mixing_ratio / 2
+        )
+
+        drho_dz = formulae.hydrostatics.drho_dz(
+            p=p,
+            T=T,
+            water_vapour_mixing_ratio=water_vapour_mixing_ratio,
+            lv=formulae.latent_heat_vapourisation.lv(T),
+            d_liquid_water_mixing_ratio__dz=(
+                self.delta_liquid_water_mixing_ratio / dz_dt / dt
+            ),
+        )
+        drhod_dz = drho_dz
+
+        self.particulator.backend.explicit_euler(self._tmp["z"], dt, dz_dt)
+        self.particulator.backend.explicit_euler(
+            self._tmp["rhod"], dt, dz_dt * drhod_dz
+        )
+
+        self.mesh.dv = formulae.trivia.volume_of_density_mass(
+            (self._tmp["rhod"][0] + self["rhod"][0]) / 2, self.mass_of_dry_air
+        )


This method is nearly identical to Parcel.advance_parcel_vars(), differing only in the computation of dz_dt. Consider extracting dz_dt logic into a protected method like _compute_dz_dt() in the base class and overriding it here. This would eliminate duplication and improve maintainability.

Good idea 👍.
Fixed here

MarcinZ20 · 2025-06-17T12:11:24Z

tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py

Parametrize the test

Firstly, in my opinion, the test is thorough and checks the critical part of your model. Great work on that! 👍🏼

With that being said, I've noticed that to achieve a better readability, you could make use of @pytest.mark.parametrize decorator and parametrize the test_figure_2_replication_accuracy function.

@pytest.mark.parametrize(...) def test_figure_2_replication_accuracy(self): ## code

To make it work, I believe that you should also somehow extract the sampling logic (that happens directly in the test) to the @pytest.fixture:

@pytest.fixture def sample_groundtruth_cases(): ## code

You can read some more on the use of fixtures here: pytest-fixtures

From my point of view, these changes will result in:

Cleaner code (makes the code more 'pythonic')

Separation of concerns (decoupling sampling logic from the assertion logic)

Easier test readability and debugging (each parametrized case will show as a separate test)

Done in this commit

MarcinZ20 · 2025-06-17T12:43:58Z

examples/PySDM_examples/Loftus_and_Wordsworth_2021/figure_2.ipynb

Code documentation

As the notebook is nicely constructed and easy to read, I think that some additional context as to what this example represents would be much appreciated 🙏🏼

As an example, you could perhaps provide some information about the data being used here:

radius_array = np.logspace(-4.5, -2.5, 50) * si.m RH_array = np.linspace(0.25, 0.99, 50) output_matrix = np.full((len(RH_array), len(radius_array)), np.nan)

I think that some short, markdown instructions/explanations would be more than enough for potential users, to dive into the code more quickly.

Here is an example of a notebook from PySDM about Kinetic Limitations in which the comments explain the use-case.

olkko · 2025-06-17T14:53:59Z

examples/PySDM_examples/Loftus_and_Wordsworth_2021/simulation.py

+        }
+
+        self.save(output)
+        while self.particulator.environment["z"][0] > 0 and output["r"][-1] > 1e-6:


I think it would be helpful to replace the magic number 1e-6 in stopping condition with a named constant so that its meaning is clear and it can be easily adjusted if needed. Also you can add comments explaining why this threshold was chosen.

That's true. Now the name
MIN_DROPLET_RADIUS
is used here.

haarmeggido

Overall very good job guys

haarmeggido · 2025-06-16T19:13:35Z

examples/PySDM_examples/Loftus_and_Wordsworth_2021/parcel.py

+        p = self["p"][0]
+
+        dz_dt = -self.particulator.attributes["terminal velocity"].to_ndarray()[0]
+        water_vapour_mixing_ratio = (


Using hardcoded index [0] suggests that the parcel assumes a single particle or layer and can behave incorrectly if any changes occur. Could the length of arrays be asserted here?

This is actually the code we copied from PySDM's parcel.py file.
As I am implementing the solution proposed by @makowskitomasz (using _compute_dz_dt function; this commit), this PR will not be a good place for fixing that potential issue.

haarmeggido · 2025-06-16T19:21:07Z

tests/examples_tests/Loftus_and_Wordsworth_2021/accuracy_test.py

+                        Zcloud,
+                        pcloud,
+                    )
+                except Exception as e:


The exception logs a warning but does not assert fail, should it be handled with pytest.fail() instead/in addition?

Yeah, you're absolutely right

Done in this commit

slayoo · 2025-06-18T07:01:33Z

CI fails de to pre-commit being neither installed or run:

black....................................................................Failed
- hook id: black
- files were modified by this hook

reformatted tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py

All done! ✨ 🍰 ✨
1 file reformatted, 837 files left unchanged.

trim trailing whitespace.................................................Passed
fix end of files.........................................................Passed
debug statements (python)................................................Passed
check json...............................................................Passed

For help, see: https://github.com/open-atmos/python-dev-hints/wiki/Pre-commit-FAQ

* simulation setup * uv lock * introduced all equations * oblate spheroid impl * AlienParcel class extending Parcel * PySDM package internal imports fix * Formulae bug fixes * Planet is now dataclass instead of dict * drop falling without pancake shape * cloudbase fix * figure settings * multicore support for simulation (#2) * parallel v1 * thread count fix * parallel * small parallel fixes * working multithread solution * Multi-thread + vector plots * black, black used everywhere --------- Co-authored-by: Piotr Kubala <pkubala@student.agh.edu.pl> * pip install joblib * Tests (#1) * tests I guess * uv fix * planet impl * ground truth * gen plot function * unit tests, monoticity tests, sanity tests * tests are now in tests * added ground truth files * removed leftovers * npy files * tests tidyup * tests adapted for current structure but still useless * renamed test to comply with pytest convention * fix ground_truth generation script * refactor unittests * added accuracy tests * Added saturation test for cloud base in unit tests * Refactored imports and clean up unused code in simulation and parcel modules * Enhanced accuracy tests by adding assertions for output completeness and consistency in unit tests * Refactored calculation of simulated mass fraction evaporation point to initialize variable and improved error handling --------- Co-authored-by: Hevagog <qkawiak@gmail.com> Co-authored-by: Mateusz Mazur <mazurm@student.agh.edu.pl> * Notebook not running fix - added back necessary imports m8 u removed too much * black formatter * black for ipynb * Update PySDM/physics/particle_shape_and_density/oblate_spheroid.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> * Updated tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py Using dir shadowed the built-in dir() function. Renamed this variable (e.g., root_dir) to avoid confusion. Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> * Fixed pressure values in Planet subclasses to use correct scientific notation * Updated method signature for proper binding. Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> * Refactored Settings class to remove unused 'coord' parameter and updated diffusion_coordinate assignment * Removed unused Lofus et al. 2021 implementations from particle shape and density, terminal velocity, and ventilation modules * plots warnings fixes * removed files rogers and yau * notebook rerun * notebook fix * minor linter fixes --------- Co-authored-by: emmacware <ecware@ucdavis.edu> Co-authored-by: Piotr Kubala <pkubala@student.agh.edu.pl> Co-authored-by: Hevagog <qkawiak@gmail.com> Co-authored-by: Mateusz Mazur <mazurm@student.agh.edu.pl>

…lanetary pressure units

Introduced const for MIN DROPLET RADIUS

piotrkubala · 2025-06-23T21:25:18Z

CI fails de to pre-commit being neither installed or run:

black....................................................................Failed
- hook id: black
- files were modified by this hook

reformatted tests/examples_tests/Loftus_and_Wordsworth_2021/ground_truth/gen_figure.py

All done! ✨ 🍰 ✨
1 file reformatted, 837 files left unchanged.

trim trailing whitespace.................................................Passed
fix end of files.........................................................Passed
debug statements (python)................................................Passed
check json...............................................................Passed

For help, see: https://github.com/open-atmos/python-dev-hints/wiki/Pre-commit-FAQ

Fixed.

codecov · 2025-06-24T20:07:37Z

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 85.66%. Comparing base (07512fe) to head (dba682d).

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #1648   +/-   ##
=======================================
  Coverage   85.66%   85.66%           
=======================================
  Files         403      403           
  Lines        9771     9773    +2     
=======================================
+ Hits         8370     8372    +2     
  Misses       1401     1401

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

slayoo mentioned this pull request Jun 10, 2025

Physics of Falling Raindrops in Diverse Planetary Atm. slayoo/AGH-MPS-2025#4

Open

kubijaku reviewed Jun 16, 2025

View reviewed changes

makowskitomasz reviewed Jun 17, 2025

View reviewed changes

MarcinZ20 reviewed Jun 17, 2025

View reviewed changes

olkko reviewed Jun 17, 2025

View reviewed changes

haarmeggido approved these changes Jun 17, 2025

View reviewed changes

lursz and others added 3 commits June 19, 2025 22:35

fix exception messages and more fixtures in tests

29168e0

refactor: streamlined dz_dt computation in Parcel class and updated p…

f55d91a

…lanetary pressure units

mzsuetam force-pushed the rain-on-me branch from 88fcb6e to f55d91a Compare June 19, 2025 20:35

mzsuetam and others added 7 commits June 19, 2025 22:57

style: improved code formatting and readability

bb15488

🔧 PR fixes

3f4f9a8

Introduced const for MIN DROPLET RADIUS

some linter fixes

7354f93

accuracy test linter

1d29d0c

accuracy linter fixed

70464d6

linter fixed

736038c

pre-commit run

217919e

piotrkubala added 5 commits June 24, 2025 01:48

linter other fixes

b8f819d

bibliography fix

e0874d5

parcel small fix

4d0299a

linter fixes

7f114a7

last linter

dba682d

Added example based on Loftus and Wordsworth 2021 figure 2 #1648

Are you sure you want to change the base?

Added example based on Loftus and Wordsworth 2021 figure 2 #1648

Uh oh!

Conversation

lursz commented Jun 10, 2025

Uh oh!

slayoo commented Jun 10, 2025

Uh oh!

kaitlyn-loftus commented Jun 10, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

makowskitomasz left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mzsuetam Jun 19, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Parametrize the test

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Code documentation

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

haarmeggido left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mzsuetam Jun 19, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

slayoo commented Jun 18, 2025

Uh oh!

piotrkubala commented Jun 23, 2025

Uh oh!

codecov bot commented Jun 24, 2025

Codecov Report

Uh oh!

Uh oh!

kaitlyn-loftus commented Jun 10, 2025 •

edited

Loading

mzsuetam Jun 19, 2025 •

edited

Loading

mzsuetam Jun 19, 2025 •

edited

Loading