# Time stepping formulas: velocity prediction with sympleptic scheme

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?

## Comments

"During the predictor stage the values of

accelerationand density" should be "During the predictor stage the values ofpositionand density". Is that correct?Yes, thank you

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

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

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

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.

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 Communications187 (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.

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.

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.

@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

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

@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.

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

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

notchanged 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.

@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

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: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;Cited LiteratureComputer Physics Communications187 (February 2015): 204–16. https://doi.org/10.1016/j.cpc.2014.10.004.Computers & Fluids179 (January 2019): 632–44. https://doi.org/10.1016/j.compfluid.2018.11.020.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.Reports on Progress in Physics68, no. 8 (1 August 2005): 1703–59. https://doi.org/10.1088/0034-4885/68/8/R01.Fluid Mechanics and the SPH Method. Oxford University Press, 2015. https://global.oup.com/academic/product/fluid-mechanics-and-the-sph-method-9780198744238Tertiary sourcesHello 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

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

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)

andmodified the last one (32c) into (33d). The algorithm has just changed. More simply, I would recognize it as a very earnestsingle-step two-stage scheme.Equations (33) advance the solution (position,velocity) over a single time step with a first-order

forward Euler schemefor the first half time step (33a-b), and with a second-ordermidpoint rulefor 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 anEuler-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.

Oops, correction. This sentence

Equations (33) advance the solution (position,velocity) over a single time step with a first-order

forward Euler schemefor the first half time step (33a-b), and with a second-ordermidpoint rulefor 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 schemefor the first half time step (33a-b), and with a second-ordermidpoint rulefor the~~second half~~whole time step (33c-d)As a note, the difference between (33c) and (33d) is that

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);midpoint velocityv(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.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.

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