Time stepping formulas: velocity prediction with sympleptic scheme

sph_tudelft_nl

Referring to the Users' Guide and the subsection 3.7.2 on the sympleptic scheme (https://github.com/DualSPHysics/DualSPHysics/wiki/3.-SPH-formulation#37-time-stepping), what is the precise formula used to calculate the velocity at the time instant n+1/2?

This could be an addition to formula (33) if I am not mistaken. Corrections welcome.

Also, the text "During the predictor stage the values of acceleration and density" should be "During the predictor stage the values of position and density". Is that correct?

Alex

"During the predictor stage the values of acceleration and density" should be "During the predictor stage the values of position and density". Is that correct?

Yes, thank you

sph_tudelft_nl

Glad to help. Please one more question: what is the expression for v(n+1/2) eventually implemented in the code?

Asalih3d

I am also interested in that. Tried to search the code but couldn't find the correct place where the calculation was performed.

Kind regards

jmdalonso

You can find the step calculation using Symplectic in JSphCpuSingle::ComputeStep_Sym(). This method calls JSphCpu::ComputeSymplecticPre() for predictor step and JSphCpu::ComputeSymplecticCorr() for corrector step of Symplectic. You can find similar methods in the GPU code.

Best regards

sph_tudelft_nl

Thanks for all contributions, but I would like to bring the topic back to the original question.

What is the precise formula used to calculate the velocity at the time instant n+1/2?

Or:

What is the expression for v(n+1/2) implemented in the code?

A description of it lacks in the Users' Guide, namely around formula 33, if I am not mistaken.

An answer of the developers' is the most authoritative. This piece of information provides an important explanation of the behaviour of simulations. If deemed applicable, an edit of the Users' Guide would help users understand what DualSPhysics does under the hood.

Thanks again for addressing this post.

sph_tudelft_nl

Apologies for urging a clarification once again. The same issue is present in Crespo, A.J.C., J.M. Domínguez, B.D. Rogers, M. Gómez-Gesteira, S. Longshaw, R. Canelas, R. Vacondio, A. Barreiro, and O. García-Feal. 2015. ‘DualSPHysics: Open-Source Parallel CFD Solver Based on Smoothed Particle Hydrodynamics (SPH)’. Computer Physics Communications 187 (February): 204–16. https://doi.org/10.1016/j.cpc.2014.10.004.

Directions to an accurate report, and its inclusion in the Users' Guide, will be most gratefully appreciated. Thanks in advance for this.

sph_tudelft_nl

Is it possible that the "predictor/corrector" time integration for fluid particles implemented in the code follows the algorithm below? Would someone from the Developers Team please confirm/correct?

Thanks in advance for looking into this and letting know.

sph_tudelft_nl

This flow possibly improves the analysis one post up.

Would someone with a fresh pair of eyes and knowledge of the code please double-check whether this analysis is correct?

Thanks in advance for looking into this and letting know.

Asalih3d

@sph_tudelft_nl

I just started looking into this too, and I agree with you - it is quite confusingly setup and I am not sure either what is actually done.

EDIT: From your source, I cannot figoure out how "a_n" is computed? I assumte it is acceleration, but what acceleration exactly?

EDIT2: Think I got it. It is the current acceleration of the particle.

Kind regards

Asalih3d

To my understanding this is how it should be done in theory;

I have no clue why the DualSPHysics code uses the (2-x)/(2+x) in calculating density, but what I have tried to outline is in theory what should happen. If there is a mistake, please do correct me.

Kind regards

sph_tudelft_nl

@Asalih3d Thanks for taking the bother of making an independent analysis. In the meantime I have done a little work on this (really a little) and I hope that I can share part of it some time soon.

Asalih3d

No problem! Looking forward to seeing your little work, whenever it is ready :-)

In the mean time if anyone has more knowledge about the theory / why it is a bit different in DualSPHysics code I would love to hear as well.

Kind regards

sph_tudelft_nl

Regarding question by @Asalih3d on the time integration of density, I would like to add that combining the two expressions of formula (1.6) in my scheme above gives:

Please double check for safety.

If this is correct, the factor before D*delta t is unity when rho(n)=rho(n+1/2) that is when is the density has not changed in the first half time step, which seems a special situation to me. So this factor, which is compounded in that epsilon quantity, does play a role in the density update in most of the cases, when a change of density is needed according to the first stage (formerly, the predictor step). So, E&OE, this formula is a modified midpoint rule updating the density from time level n to time level n+1.

In a shallow search I could not find any reference to this epsilon expression. It would be helpful if developers could give directions to the rational of this expression. Thanks for that.

Asalih3d

@sph_tudelft_nl

Thanks for looking more into it, but I don't feel in a position to confirm/deny, since I don't even know what the correct method is yet.. :-)

Tried reading in Monaghan's 2005 paper, but I could not find the relevant part I assume.

Kind regards, hope we will understand it better some time

sph_tudelft_nl

This is the result of the analysis stimulated by the first questions up in this discussion and partly reported above. The list below consists of statements to critique (also in the sense that they might be wrong) as well as of misunderstandings or understatements to clear. Unfortunately, we at Delft cannot prioritize the lines of investigation that opened up in this occasion. Therefore, we are glad to share these findings with the general community; others may want to correct or integrate to advance the knowledge record.

With reference to Section 2.7 on Time stepping in the article by Crespo et al (2015) and, specifically, to updating the state of fluid particles:

1. The property of symplecticness, or symplecticity, applies equally well to the Verlet scheme in Section 2.7.1; see Section 2 of Leimkuhler et al. (1996);
2. The description of the ‘symplectic predictor-corrector’ in Section 2.7.2 contains a few oversights concerning formulas and names;
3. The computer program associated to the article is DualSPHysics version 3.00 (17-12-2013), stored at https://data.mendeley.com/datasets/mcgc3636xr/1). The predictor-corrector formulas implemented in the subroutines JSphCpu::ComputeSymplecticPre and JSphCpu::ComputeSymplecticCorr differ from those described in Section 2.7.2;
4. The ‘symplectic predictor/corrector’ scheme coded in versions 4.4.003 (10-04-2019) and 5.0.100 (17-04-2020) differs from that of version 3 and, in turn, from the description in the article;
5. Leimkuhler et al. (1996) contains neither the symplectic predictor/corrector formulas published in the article nor those implemented in the codes;
6. Ditto for Monaghan (2005);
7. The naming ‘symplectic predictor/corrector’ does not apply to the algorithm coded in DualSPHysics 4.4 and 5.0. That algorithm is, rather, a single-step two-stage method: stage one is the integration over the first half step, stage two is the integration over the second half step;
8. Among the literature listed in https://dual.sphysics.org/index.php/references/ and published since 2015, discounting the works dealing with incompressible SPH and/or using a Verlet scheme only, only Green et al. (2019) report the time integration method effectively implemented in the code;
9. Possibly, it has still to be shown that the single-step two-stage scheme implemented in the code is symplectic;
10. Canonical symplecticness concerns the time integration of position and velocity (Violeau, 2015, Section 5.4); time integration of density should be explained separately or in the light of a non-canonical treatment of symplecticness;
11. The current version of the Users’ Guide probably features the same strengths and weaknesses as the article.

Cited Literature

• Crespo, A.J.C., J.M. Domínguez, B.D. Rogers, M. Gómez-Gesteira, S. Longshaw, R. Canelas, R. Vacondio, A. Barreiro, and O. García-Feal. ‘DualSPHysics: Open-Source Parallel CFD Solver Based on Smoothed Particle Hydrodynamics (SPH)’. Computer Physics Communications 187 (February 2015): 204–16. https://doi.org/10.1016/j.cpc.2014.10.004.
• Green, M.D., R. Vacondio, and J. Peiró. ‘A Smoothed Particle Hydrodynamics Numerical Scheme with a Consistent Diffusion Term for the Continuity Equation’. Computers & Fluids 179 (January 2019): 632–44. https://doi.org/10.1016/j.compfluid.2018.11.020.
• Leimkuhler, B.J., S. Reich, and R.D. Skeel. ‘Integration Methods for Molecular Dynamics’. In Mathematical Approaches to Biomolecular Structure and Dynamics, edited by Jill P. Mesirov, Klaus Schulten, and De Witt Sumners, 82:161–85. The IMA Volumes in Mathematics and Its Applications. New York, NY: Springer New York, 1996. https://doi.org/10.1007/978-1-4612-4066-2_10.
• Monaghan, J.J. ‘Smoothed Particle Hydrodynamics’. Reports on Progress in Physics 68, no. 8 (1 August 2005): 1703–59. https://doi.org/10.1088/0034-4885/68/8/R01.
• Violeau, D. Fluid Mechanics and the SPH Method. Oxford University Press, 2015. https://global.oup.com/academic/product/fluid-mechanics-and-the-sph-method-9780198744238

Tertiary sources

• Predictor-corrector method: https://en.wikipedia.org/wiki/Predictor%E2%80%93corrector_method
• Symplectic integrator: https://en.wikipedia.org/wiki/Symplectic_integrator

ben

Hello Giordano,

The scheme implemented in DualSPHysics is a symplectic Position Verlet scheme. I can recommend it as an instructive exercise to prove it. Alas, the code and its terminology is written for convenience, not as a self-contained tutorial in symplectic methods. Some future clarifying statements in the guide-wiki might help remove queries from users.

Best regards,

Ben

ben

Dear Asalih3d,

The density is a scalar and obeys the same evolution as position. The expression for density timestepping used in the code comes from Parshikov et al. 2000 which is identical to a mid-point symplectic after some rearranging and when first implemented in in 2008, it was found to alleviate some precision issues that are no longer relevant.

I hope this helps. Best regards,

Ben

sph_tudelft_nl

https://forums.dual.sphysics.org/discussion/comment/3758#Comment_3758

Thanks @ben for pursuing this topic and for the work made at revising the description of the time stepping in the User’s Guide aka wiki. This makes it possible to understand and recognize what DualSPHysics does.

In section 3.2.7 of the wiki Formula (32) is recognizably a position Verlet scheme: its drift-kick-drift procedure is almost iconic. (Jargon-buster: drift is the displacement update after a ‘drifting velocity’, and kick is the velocity update after a ‘kicking acceleration’.)

However, I doubt that one should keep on calling the algorithm (33) a position Verlet scheme after having added a second equation (33b) and modified the last one (32c) into (33d). The algorithm has just changed. More simply, I would recognize it as a very earnest single-step two-stage scheme.

Equations (33) advance the solution (position,velocity) over a single time step with a first-order forward Euler scheme for the first half time step (33a-b), and with a second-order midpoint rule for the second half time step (33c-d). If I am not wrong, its overall accuracy is still second-order. I refrain from anticipating the symplectic properties of (33) but I am not understating either they do not apply. Nonetheless, I would plainly call the second method used by DualSPHysics an Euler-and-midpoint rule/scheme/method, or the like, rather than a symplectic position Verlet scheme.

As to the code naming, in science and technology any naming based on conventions rather than on convenience greatly helps community members understand one another and save time for action. Thanks for any little effort put into this.

sph_tudelft_nl

Oops, correction. This sentence

Equations (33) advance the solution (position,velocity) over a single time step with a first-order forward Euler scheme for the first half time step (33a-b), and with a second-order midpoint rule for the second half time step (33c-d)

should be

Equations (33) advance the solution (position,velocity) over a single time step with a first-order forward Euler scheme for the first half time step (33a-b), and with a second-order midpoint rule for the second half whole time step (33c-d)

As a note, the difference between (33c) and (33d) is that

• the midpoint acceleration (force per unit mass) F(n+1/2) is computed from the SPH interactions -- so indirectly via the midpoint position and velocity determined with (33a-b);
• the midpoint velocity v(n+1/2) is computed from an average of (33b) and (33c).

Both (33c-d) follow the pseudo equation: updated value = old value + rate of change at the mid time level x whole time step size.

ngkhaiching2000

Hi DualSPHysics developer:

Can I get your advice on the following sympletic scheme implemented in DualSPHysics

(1) Is the force term F_n+0.5 calculated from rho_n+0.5, OR rho_n+1 (Eq. 34b)?

(2) For the term R_n at the predictor stage (Eq 34a), is it obtained from the value of R_n+0.5 obtained from the previous correction stage? 

Please refer to the green arrow below.

Thanks.

Alex

@gfourtakas when you have time, can you check this, please?