Why can I not achieve consistency?

I am setting up an experiment to validate DualSPHysics's wave propagation, which I later want to use for a big WEC simulation.

I do however, have some problems in achieving any consistency.

Setup is the following (2D):

  • T = 1.42s
  • λ = 3.14 m
  • d = 1.682m
  • H = 0.11m
  • d/λ = 0.535 (deep water)
  • H/λ = 0.035 (linear wave)

I now conducted 24 runs with different values of H/dp (5, 10, 15, 20, 25, 30) and different formulations of the viscosity term (SPS, α=0.01, α=0.001, α=0.0001).

An example of a phase averaged period (by upcrossing analysis) is given in the following for SPS and 3 resolutions.

From the plot it is evident that the numerical solution is in fact diverging from the analytical one with increasing H/dp. Similar behaviour can be observed for the other viscosity formulations. To get an overview of all runs I calculated the "normalized (by H) root mean square deviation" (NRMSD) at WG1 and plotted it.

Again I am very confused by the fact that my numerical solution is not converging towards the analytical one with increasing spatial resolution. In fact it seems to perform worse for e.g. α=0.001 and SPS.

So I have two questions:

  • Why can I not achieve consistency? Meaning numerical solution converges towards analytical one for better resolution.
  • Why does my numerical wave seem to be shifted "up left" in comparison to the analytical one (see averaged period plot). I have been trying to figure this out for quite a while now but can't make any progress.

Example .xml file and Run output for SPS H/dp=5:


  • Dear,

    I see you are using 1st order generation and comparing with 1t order theoretical time series. Am I wrong?

    I see that 2nd order prescribes a slightly sharper wave crest and differences also in the trough. If my understanding is correct, please use 2nd order generation.

    Besides, so not use Laminar+SPS in 2D, since it's development was meant for 3D. I would rather recommend to use artificial viscosity with alpha=0.01. Rota Roselli et al. (2018) proved that alpha=0.01 was optimal value for wave propagation, together with H/dp=10.

  • You are right. I was generating linear first order waves. Checking Le Mehaute Diagram again, I found that I am actually in the second order domain:

    That was a very silly mistake of mine, but sadly, after redoing the simulations with second order generation I am observing very similar results:

    Again I am confused about the fact that my results are not getting better with increasing H/dp, and seem to be the best for H/dp=15. From my understanding this is what was defined as consistent behavior (error going down with better discretization). But I realized I had a mistake with terminology here:

    • Consistency: The truncation errors of a numerical approximation of a differential equation are going to zero with dx->0 and dt->0
    • Stability: A numerical solution method is said to be stable if it does not magnify the errors that appear in the course of numerical solution process.
    • Convergence: A numerical method is said to be convergent if the solution of the discretized equations tends to the exact solution of the differential equation as the grid spacing tends to zero. Usually that would be the case for a consistent and stable scheme. (convergence = consistency + stability)

    Thus, my original question has to be adapted to "Why can I not achieve convergence?".

    Adding up on this I read your suggestion about Roselli (2018) and found the following sentence: "However, convergence has not yet been completely proven in SPH and the choice of this parameter influences the contribution of artificial viscosity terms.". Also I read that there are consistency problems in SPH, due to the anti-symetry in the discretized kernel sampling (heterogeneously placed particles within the kernel support). Therefore, I assume I can not achieve convergence, because consistency is not fulfilled.

    I again plotted the free surface elevation for some example H/dp. But I found that by the averaging procedure, which I used in the original question, I am in deed plotting without the phase error. So this time the plots are with phase error. The main error of the NRMSD plot from above therefore, seems to be coming from this phase shift.

    alpha = 0.01

    https://forums.dual.sphysics.org/uploads/434/1YOJVFEABVAW.png There was an error displaying this embed.

    SPS (which seems to be performing quite well also in 2D):

    https://forums.dual.sphysics.org/uploads/021/X61V5W2E1LJH.png There was an error displaying this embed.

    So my questions:

    • Is it even expected that I am closer to the analytical solution for higher H/dp?
    • Why do I observe the shift in the wave?
  • edited June 28

    Sorry I forgot to quote. I am not sure if any notification is sent if I am not quoting.

Sign In or Register to comment.