the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 05 Jun 2026
Pysammos 1.0.0: a discrete-to-continuum transformation Python tool to analyse the rheology of granular materials
Abstract. Granular flow processes involving different interstitial fluids and coupling regimes are widespread in both natural systems (e.g., landslides, rock avalanches, river sediment transport) and industrial applications (e.g., aggregates in concrete manufacturing, powder technology, animal feed). Despite this, the behaviour of complex granular flows is not fully understood. Modelling of granular media through software packages that couple the Discrete Element Method with Computational Fluid Dynamics (DEM-CFD) enables a comprehensive description of granular small-scale mechanics through simulations of particle-particle and particle-fluid interactions at high temporal and spatial resolution. These approaches are pivotal in the development of constitutive models that represent the bulk rheology of granular media. While DEM-CFD simulations provide particle-scale information (e.g., particle velocities and forces), extracting continuum fields requires a discrete-to-continuum (D2C) transformation that applies a mathematical approach termed coarse-graining. Although some DEM software packages include built-in D2C capabilities, others, such as MFiX-DEM, do not. Consequently, users are often required to develop custom D2C workflows or adapt simulation outputs to the requirements of other D2C tools. Hence, we introduce Pysammos, a Python package that performs discrete-to-continuum transformations and is designed to be user-friendly, open-source, and computationally efficient. Pysammos is able to process polydisperse granular mixtures of any particle shape, while also offering the option to analyse different particle phases separately. It post-processes output files from MFiX-DEM software and produces vtkhdf outputs ready for visualisation in ParaView, as well as a more generic h5 format for further data analysis. Pysammos is able to operate on standard desktop computers as well as on HPC systems. Finally, we showcase a variety of exemplar applications such as sediment erosion, crystals and magma in a conduit, bedload transport and impact cratering.
-
Notice on discussion status
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
-
Preprint
(3958 KB)
-
Supplement
(415 KB)
-
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(3958 KB) - Metadata XML
-
Supplement
(415 KB) - BibTeX
- EndNote
- Final revised paper
Journal article(s) based on this preprint
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2026-2591', Alexandre Sac-Morane, 15 Jun 2026
-
AC1: 'Reply on RC1', Claudia Elijas-Parra, 11 Sep 2026
We thank the referee 1 for their positive comments and thoughtful suggestions below. Each point has been addressed individually, with changes to the manuscript text shown in orange text. Unchanged manuscript associated with these changes is shown in blue text. The responses to each of the comments are provided in the supplement file "RC1_comments_response.pdf".
-
AC1: 'Reply on RC1', Claudia Elijas-Parra, 11 Sep 2026
-
RC2: 'Comment on egusphere-2026-2591', Francois Guillard, 13 Aug 2026
This paper describe a coarse-graining software, Pysammos, dedicated to the spatial averaging of discrete element modelling data, in particular the one produces by the software MFiX-DEM. The principle of the method and software are explained, as well as its specificities, and its performance are compared with other coarse-graining codes and with different geometries and simulation conditions.
The paper is overall well written and detailed, and can be useful as an entry point to coarse-graining, while proposing additional features not often present in already existing solutions. Some of the specific aspects, such as the phase identification and the validation with non-spherical particles, can be useful in a broad range of contexts. I only have some minor comments and questions below.1/ In the simulations of particles suspended in a magma, it is not clear how the flow is created, and if it is a one-way or two-ways interaction between the fluid and particles. It would be interesting as well to provide the volume fraction field, which may have some spatial variations. Also, the simulations are indicated as 2D, but the particles seems to be three-dimensional?
2/ The results on the stress under the pile (4.4.4) are quite surprising. Experimental observations show essentially hydrostatic pressure, sometimes with a dip of pressure below the center, but both appear quite different from the results reported. Even if frictional redistribution of the stress occurs (similar to the Janssen effect), the weight of the grains still needs to be supported somehow. It would be good to also look at the shear stresses to assess whether this frictional redistribution does occur. Also, Fig 13a may be better plotted with a linear scale, since again the expectation would be a linear stress increase with depth on average.
3/ The calculation of the contact stress involves integrating the windowing function over the branch vector, as shown in equation 5. Is this approximated through averaging the values at the center, as is commonly done, or is the integral numerically solved? How? This may explain some of the differences pointed out in σc later on in the paper.
4/ There are many different clustering algorithm that exist, how was the particular one used chosen?
5/ line 144: the mathematical detail of the Dirac delta function is inaccurate. Technically, the Dirac function is infinity at the particle centre, and 0 otherwise, such as its integral over the whole space is 1.
6/ Equations 5 and C4: the tensor product, or Kronecker product that is used for the vector in those equations, should be specified. It is usually denoted with ⊗.
7/ line 193: the use of only equilibrium states to justify the replacement of ensemble averaging with time averaging appears quite restrictive. Most granular systems are out of equilibrium, and time averaging can still be meaningfully performed is the systems are at steady state.
8/ Fig. 11b: the scale shows velocities going up to 10^30 m/s, which is much higher than the speed of light. There is probably an issue with the scale.
9/ Fig 15 caption: the explanation of the "greyish" colour is incorrect. Ty is not zero at these locations, but undefined, since the volume fraction is zero (no particles). As indicated by equation C13, this would correspond to dividing a zero kinetic stress tensor trace by the zero density at that location, which numerically should result in a Not a Number.
10/ line 885: it would be good to define D right after eq. C13, rather than on line 895.
11/ eq C19: the equation uses d̄, which is not defined. Is it the same as D43 (eq. C6)?
Citation: https://doi.org/10.5194/egusphere-2026-2591-RC2 -
AC2: 'Reply on RC2', Claudia Elijas-Parra, 11 Sep 2026
We thank the referee 2 for their positive comments and thoughtful suggestions below. Each point has been addressed individually, with changes to the manuscript text shown in orange text. Unchanged manuscript associated with these changes is shown in blue text. The responses to the referee comments are provided in the supplement file "RC2_comments_response.pdf".
-
AC2: 'Reply on RC2', Claudia Elijas-Parra, 11 Sep 2026
Peer review completion
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2026-2591', Alexandre Sac-Morane, 15 Jun 2026
This contribution presents Pysammos, a software that can be used to post-process the results of any Discrete Element Modelizations (DEM). Indeed, the discrete particle data obtained with the DEM can be transformed into continuous fields, more accessible for visualization purposes and deeper analysis. The architecture, the abilities, the validity, and the performances of Pysammos are depicted in this paper. This method appears user-friendly, computationally efficient, and flexible to the DEM software used.
The article is well written and well structured.
Even if I am positive about this article, I have a few queries to consider during this minor revision.
L. 55: I would change “force displacement law” to “contact law”.
L. 116: There is no verb in “Similarly, to enhance DEM analysis by extracting continuum fields without the need to handle inner-level source code”. Could you reformulate this sentence to ease the reading?
L. 149: I would not refer to Eq. 9 solely, as you use multiple functions for the coarse-graining function \Psi. Is a reference to Section 3.4.3 more accurate?
Eq. 5: Could you check this Equation? In the current version, an individual contact is computed twice. I guess that the second sum should start at j=i+1, in agreement with the formulation depicted in Weinhart et al. (2012).
L. 170: It appears in L. 385 that you made a different choice than other D2C software (MercuryDPM, for instance) for the particle search (particle vs contact point). A justification of this assumption can be relevant here.
L. 264: I would start with the Section “Data reading and handling” as the format of the input data is described. Please consider reorganizing your Sections.
Figure 6: Could you check this plot? Indeed, the Heavyside function does not verify the formulation detailed in Eq. 11. The function should be =0 for w<r, while the plot shows a non-null value for w<r<2w. This verification should be conducted for the Lucy and Gaussian functions.
Figure 6: I guess that W(r) should be Psi(r). If not, could you specify the definition of W(r)?
Figure B2: Same remark as Figure 6: I guess that W(r) should be Psi(r). If not, could you specify the definition of W(r)?
L. 529: This depth-independence of the vertical pressure is controversial. I agree with the Authors that Fig. 13 does not show a clear increase in the pressure. However, the relation P=rho.g.h is well established in the literature. I guess this relation is retrieved for the configuration depicted if a larger height h or a larger density rho is considered. Could you add limitations to the conclusion raised concerning the homogeneous vertical pressure with depth?
Figure 12: Could you consider modifying the scale used for the right map? The color assigned to the value -6.6 does not seem to be employed in the map. The range used for the scale can be shorter.
L. 574: As discussed with Section 4.4.3, the default value of w (=0.75*d_43) may not be the accurate choice to make. Maybe a default value of w that depends on the mean distance of the particles is a more resilient selection. In particular, if this code is used by a standard user (see the definition of the profile in Section 3.1.3). This remark is more a hint to the Authors, the text of the article does not have to be changed. I have not tried; maybe it is a bad idea.
L. 683: The consideration of the particle-wall interactions seems in Equation 5 to be already conducted (with the particles between N+1 and N+Nb). If it is finally not the case, could you delete this part of the Equation?
Citation: https://doi.org/10.5194/egusphere-2026-2591-RC1 -
AC1: 'Reply on RC1', Claudia Elijas-Parra, 11 Sep 2026
We thank the referee 1 for their positive comments and thoughtful suggestions below. Each point has been addressed individually, with changes to the manuscript text shown in orange text. Unchanged manuscript associated with these changes is shown in blue text. The responses to each of the comments are provided in the supplement file "RC1_comments_response.pdf".
-
AC1: 'Reply on RC1', Claudia Elijas-Parra, 11 Sep 2026
-
RC2: 'Comment on egusphere-2026-2591', Francois Guillard, 13 Aug 2026
This paper describe a coarse-graining software, Pysammos, dedicated to the spatial averaging of discrete element modelling data, in particular the one produces by the software MFiX-DEM. The principle of the method and software are explained, as well as its specificities, and its performance are compared with other coarse-graining codes and with different geometries and simulation conditions.
The paper is overall well written and detailed, and can be useful as an entry point to coarse-graining, while proposing additional features not often present in already existing solutions. Some of the specific aspects, such as the phase identification and the validation with non-spherical particles, can be useful in a broad range of contexts. I only have some minor comments and questions below.1/ In the simulations of particles suspended in a magma, it is not clear how the flow is created, and if it is a one-way or two-ways interaction between the fluid and particles. It would be interesting as well to provide the volume fraction field, which may have some spatial variations. Also, the simulations are indicated as 2D, but the particles seems to be three-dimensional?
2/ The results on the stress under the pile (4.4.4) are quite surprising. Experimental observations show essentially hydrostatic pressure, sometimes with a dip of pressure below the center, but both appear quite different from the results reported. Even if frictional redistribution of the stress occurs (similar to the Janssen effect), the weight of the grains still needs to be supported somehow. It would be good to also look at the shear stresses to assess whether this frictional redistribution does occur. Also, Fig 13a may be better plotted with a linear scale, since again the expectation would be a linear stress increase with depth on average.
3/ The calculation of the contact stress involves integrating the windowing function over the branch vector, as shown in equation 5. Is this approximated through averaging the values at the center, as is commonly done, or is the integral numerically solved? How? This may explain some of the differences pointed out in σc later on in the paper.
4/ There are many different clustering algorithm that exist, how was the particular one used chosen?
5/ line 144: the mathematical detail of the Dirac delta function is inaccurate. Technically, the Dirac function is infinity at the particle centre, and 0 otherwise, such as its integral over the whole space is 1.
6/ Equations 5 and C4: the tensor product, or Kronecker product that is used for the vector in those equations, should be specified. It is usually denoted with ⊗.
7/ line 193: the use of only equilibrium states to justify the replacement of ensemble averaging with time averaging appears quite restrictive. Most granular systems are out of equilibrium, and time averaging can still be meaningfully performed is the systems are at steady state.
8/ Fig. 11b: the scale shows velocities going up to 10^30 m/s, which is much higher than the speed of light. There is probably an issue with the scale.
9/ Fig 15 caption: the explanation of the "greyish" colour is incorrect. Ty is not zero at these locations, but undefined, since the volume fraction is zero (no particles). As indicated by equation C13, this would correspond to dividing a zero kinetic stress tensor trace by the zero density at that location, which numerically should result in a Not a Number.
10/ line 885: it would be good to define D right after eq. C13, rather than on line 895.
11/ eq C19: the equation uses d̄, which is not defined. Is it the same as D43 (eq. C6)?
Citation: https://doi.org/10.5194/egusphere-2026-2591-RC2 -
AC2: 'Reply on RC2', Claudia Elijas-Parra, 11 Sep 2026
We thank the referee 2 for their positive comments and thoughtful suggestions below. Each point has been addressed individually, with changes to the manuscript text shown in orange text. Unchanged manuscript associated with these changes is shown in blue text. The responses to the referee comments are provided in the supplement file "RC2_comments_response.pdf".
-
AC2: 'Reply on RC2', Claudia Elijas-Parra, 11 Sep 2026
Peer review completion
Journal article(s) based on this preprint
Data sets
Example granular material simulation data for the Pysammos Python package C. Elijas-Parra et al. https://doi.org/10.5281/zenodo.19351802
Model code and software
Pysammos v1.0.0 C. Elijas-Parra et al. https://doi.org/10.5281/zenodo.19355667
Pysammos - GitHub repository C. Elijas https://github.com/Claudia-Elijas/pysammos
Pysammos - Documentation C. Elijas https://claudia-elijas.github.io/pysammos/
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 206 | 157 | 39 | 402 | 57 | 18 | 28 |
- HTML: 206
- PDF: 157
- XML: 39
- Total: 402
- Supplement: 57
- BibTeX: 18
- EndNote: 28
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Claudia Elijas-Parra
Eric C. P. Breard
John P. Morrissey
P. J. Zrelak
Mark Naylor
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(3958 KB) - Metadata XML
-
Supplement
(415 KB) - BibTeX
- EndNote
- Final revised paper
This contribution presents Pysammos, a software that can be used to post-process the results of any Discrete Element Modelizations (DEM). Indeed, the discrete particle data obtained with the DEM can be transformed into continuous fields, more accessible for visualization purposes and deeper analysis. The architecture, the abilities, the validity, and the performances of Pysammos are depicted in this paper. This method appears user-friendly, computationally efficient, and flexible to the DEM software used.
The article is well written and well structured.
Even if I am positive about this article, I have a few queries to consider during this minor revision.
L. 55: I would change “force displacement law” to “contact law”.
L. 116: There is no verb in “Similarly, to enhance DEM analysis by extracting continuum fields without the need to handle inner-level source code”. Could you reformulate this sentence to ease the reading?
L. 149: I would not refer to Eq. 9 solely, as you use multiple functions for the coarse-graining function \Psi. Is a reference to Section 3.4.3 more accurate?
Eq. 5: Could you check this Equation? In the current version, an individual contact is computed twice. I guess that the second sum should start at j=i+1, in agreement with the formulation depicted in Weinhart et al. (2012).
L. 170: It appears in L. 385 that you made a different choice than other D2C software (MercuryDPM, for instance) for the particle search (particle vs contact point). A justification of this assumption can be relevant here.
L. 264: I would start with the Section “Data reading and handling” as the format of the input data is described. Please consider reorganizing your Sections.
Figure 6: Could you check this plot? Indeed, the Heavyside function does not verify the formulation detailed in Eq. 11. The function should be =0 for w<r, while the plot shows a non-null value for w<r<2w. This verification should be conducted for the Lucy and Gaussian functions.
Figure 6: I guess that W(r) should be Psi(r). If not, could you specify the definition of W(r)?
Figure B2: Same remark as Figure 6: I guess that W(r) should be Psi(r). If not, could you specify the definition of W(r)?
L. 529: This depth-independence of the vertical pressure is controversial. I agree with the Authors that Fig. 13 does not show a clear increase in the pressure. However, the relation P=rho.g.h is well established in the literature. I guess this relation is retrieved for the configuration depicted if a larger height h or a larger density rho is considered. Could you add limitations to the conclusion raised concerning the homogeneous vertical pressure with depth?
Figure 12: Could you consider modifying the scale used for the right map? The color assigned to the value -6.6 does not seem to be employed in the map. The range used for the scale can be shorter.
L. 574: As discussed with Section 4.4.3, the default value of w (=0.75*d_43) may not be the accurate choice to make. Maybe a default value of w that depends on the mean distance of the particles is a more resilient selection. In particular, if this code is used by a standard user (see the definition of the profile in Section 3.1.3). This remark is more a hint to the Authors, the text of the article does not have to be changed. I have not tried; maybe it is a bad idea.
L. 683: The consideration of the particle-wall interactions seems in Equation 5 to be already conducted (with the particles between N+1 and N+Nb). If it is finally not the case, could you delete this part of the Equation?