Comparison against the 'py-pde' package based on the example: Diffusion equation with spatial dependence

[arekpaterak]

In this pull request, we extend the examples with a reproduction of the "Diffusion equation with spatial dependence" example from the py-pde library (the original example).
Besides showcasing PyMPDATA capabilities on a new example, this enables us to compare both libraries.

To reproduce the example, we have had to extend PyMPDATA itself to handle heterogeneous diffusivity.

---

We have started the work on this project in a separate repository with a different project structure so the change history doesn't fully reflect real contributions. Then we moved out changes to a fork of PyMPDATA.

• Adrian (@MyNameIsArko) primarily contributed to reproducing the example with PyMPDATA and making necessary extensions to the library.
• Wiktor (https://github.com/WiktorProsowicz) and Jan (https://github.com/yancostrishevsky) prepared the core logic of the comparison.
• Arek (@arekpaterak) integrated it with PyMPDATA library to make this direct pull request.

[michalszc]

Cell ID: 51c60e1d

It looks like the time measurement is not isolated to the solution's execution, as it includes compilation. To avoid this, consider compiling the code ahead of time - see this guide for details.

[michalszc]

As demonstrated during the laboratory classes, increasing the number of corrective iterations beyond 3 yields minimal, if any, improvement in solution accuracy and can often reduce precision, while increasing computational cost.

[arekpaterak]

Thanks for the suggestion, we changed this with the commit 373636f

[michalszc]

Missing newline at end of file

[szpajak]

Cell ID: f2ac8a37

The atol=0.5 argument in np.allclose is much larger than the actual observed max difference (about 0.12). For more meaningful validation, it think it should be smaller.

[arekpaterak]

Thanks for the suggestion, we fixed it with the commit 9daeb12

[szpajak]

Check if correct, because this gives cell-centered dx. For standard node-centered grid, use (xmax - xmin)/(args.grid_points - 1)

[WiktorProsowicz]

My friend, this is an enormously useful remark. It has been fixed in 6321be4

[szpajak]

Import is not used and can be deleted.

[arekpaterak]

Good catch, it is fixed with the commit 5dab1eb

[berkaydoruk]

you might wanna add a small assert to this check numerically. For example, assert abs(kymo[-1].sum() - kymo[0].sum()) < 1e-5 also havent seen in test_similarity file. Controling this can be important when simulating physical processes like diffusion.

[berkaydoruk]

This function do many things at once like from simulation setup and execution to kymograph processing and figure creation. To improve modularity, you may split it into smaller helper functions (e.g. one for running the solver, one for postprocessing the kymograph, one for plotting).

[PrzemyslawMaresz]

Consider adding a Parameter Sensitivity Study

The current notebook provides a compelling comparison for a single set of simulation parameters (grid_points = 64, dt = 1e-3). To make the analysis more robust and showcase the characteristics of the new implementation, it would be beneficial to add a brief parameter sensitivity study.

For example, you could run the comparison in a loop for a few different grid resolutions and plot the results. This would allow you to show:

• Convergence of the method: How does the RMSE between the two solutions decrease as the grid becomes finer?
• Performance scaling: How does the performance ratio (pympdata_time / py_pde_time) change with the problem size?

Demonstrating that the error decreases with higher resolution would significantly strengthen the claim that the PyMPDATA implementation is a valid alternative.

[PrzemyslawMaresz]

Question on the Discretization of Interfacial Diffusivity

In the __make_heterogeneous_laplacian function, the diffusivity at the cell face is calculated using an arithmetic mean.

While this is a straightforward approach, for problems involving heterogeneous media, using a harmonic mean for the interfacial diffusivity is often recommended to ensure the continuity of the flux and improve the physical accuracy of the model.

Have you considered this alternative? For diffusivity profiles with sharp gradients, the harmonic mean can provide more robust and accurate results. It might be worth comparing both approaches.

[szpajak]

I think it would be better to replace this hard-coded mu_coeff either with named constant value or to pass it as argument (with default value).

[szpajak]

You should change extra: Dict[str, any] to extra: Dict[str, Any] and import Any from typing.

[yancostrishevsky]

After an exhaustive analysis of the code and many hours of dedicated effort by the development team, we have successfully resolved this issue. d840f92

[MoudarSuliman]

The current implementation for the heterogeneous Laplacian returns early with return when options.non_zero_mu_coeff is False, but it might be clearer to raise a NotImplementedError or state in the docstring that heterogeneous diffusion is only supported when this option is enabled.

[MoudarSuliman]

Handling of edge cases (such as what happens if the diffusivity profile contains zeros or negative values) isn’t explicitly addressed—consider adding input validation or at least a warning to ensure consistent behavior for all user inputs.

[slayoo]

CI reports this failure (https://github.com/open-atmos/PyMPDATA/actions/runs/15725758755/job/44315180539?pr=578):

  File "/home/runner/work/PyMPDATA/PyMPDATA/./PyMPDATA/options.py", line 37, in Options
    dtype: np.float32 | np.float64 = np.float64,
TypeError: unsupported operand type(s) for |: 'type' and 'type'