[TurbJet]

Hi,

I am writing a simple 1D compressible Euler solver. The time integration I use is RK4. Now I am trying to do time-convergence check. But I am not sure about the exact procedure of doing this.

Does anyone can help me out?

[FMDenaro]

Originally Posted by TurbJet

Hi,

I am writing a simple 1D compressible Euler solver. The time integration I use is RK4. Now I am trying to do time-convergence check. But I am not sure about the exact procedure of doing this.

Does anyone can help me out?

Not sure about your question. For "time convergence" are you talking about a steady state?

[TurbJet]

Originally Posted by FMDenaro

Not sure about your question. For "time convergence" are you talking about a steady state?

No. Now I am solving 1D wave equation with periodic BC, so the wave will keep propagating. Thus, it will be unsteady.

What I am doing now is transport the wave with 1 period of time, and then compare with the initial solution. But I am guessing this is not a valid way to do this. But I can't find another for doing this.

Or maybe I can directly verify the convergence order of my RK4 module?

[FMDenaro]

So, you have a linear wave equation? You can compare the exact solution at any time as f(x,t)=f(x-ut,0). However, if you want to see the slope of the error curve ONLY for the time accuracy, you have to use a very refined spatial grid to eliminate the effect of the spatial term in the local truncation error.

[TurbJet]

Originally Posted by FMDenaro

So, you have a linear wave equation? You can compare the exact solution at any time as f(x,t)=f(x-ut,0). However, if you want to see the slope of the error curve ONLY for the time accuracy, you have to use a very refined spatial grid to eliminate the effect of the spatial term in the local truncation error.

Yes, this is what I am looking for and it's I have done. But I transport the wave further, through the periodic boundary and back to its initial position.

The tricky part is, since I only evaluate the interior grid points for spatial derivatives, when I am doing this "transform" x-ut, it's a little bit tricky to find this t, I have to exclude some points. In this case, my time convergence looks terrible.

Maybe I'll just transport the wave within the domain?

[FMDenaro]

Originally Posted by TurbJet

Yes, this is what I am looking for and it's I have done. But I transport the wave further, through the periodic boundary and back to its initial position.

The tricky part is, since I only evaluate the interior grid points for spatial derivatives, when I am doing this "transform" x-ut, it's a little bit tricky to find this t, I have to exclude some points. In this case, my time convergence looks terrible.

Maybe I'll just transport the wave within the domain?

Well, i think you have some bugs in the code, there is no problem using periodic BC.s. Just consider the periodicity lenght in (x-ut).

[TurbJet]

Originally Posted by FMDenaro

Well, i think you have some bugs in the code, there is no problem using periodic BC.s. Just consider the periodicity lenght in (x-ut).

It just came to me, is it possible that the error comes from: after traveling a long way, the wave profile got diffused out due to numerical dissipation and diffusion introduced by small time step size? So maybe I should make the wave travel a much shorter distance?


Also, I can't have 4th order convergence for my spatial derivative (I use 4th order scheme). I already set time step size small (~ 1E-5), temporal error should not be a problem. Even so, I get resulting convergence order ~O(1) (man this is a disaster !).

So I am completely lost. Is there anything else could affect spatial convergence?

[praveen]

Here is a suggestion. You will fix some grid size h and there is a time step coming from cfl condition, lets say it is dt.

Compute solution using dt, dt/2, dt/4, dt/8
Lets call these solutions u1, u2, u3, u4

Now we expect

|u1 - u2| = C(dt)^p
|u2 - u3| = C(dt/2)^p
etc.

So

p = log(|u2-u3|/|u1-u2|)/log(2)

All solutions u1, u2, u3, u4, ... will have similar spatial errors since they are on same mesh h. So |u1-u2| will cancel that spatial error and you can then measure the time integration error.

Please tell us what you find if you try this.

[TurbJet]

Originally Posted by praveen

Here is a suggestion. You will fix some grid size h and there is a time step coming from cfl condition, lets say it is dt.

Compute solution using dt, dt/2, dt/4, dt/8
Lets call these solutions u1, u2, u3, u4

Now we expect

|u1 - u2| = C(dt)^p
|u2 - u3| = C(dt/2)^p
etc.

So

p = log(|u2-u3|/|u1-u2|)/log(2)

All solutions u1, u2, u3, u4, ... will have similar spatial errors since they are on same mesh h. So |u1-u2| will cancel that spatial error and you can then measure the time integration error.

Please tell us what you find if you try this.

Thanks for the advice. But I tried this one, the order is around O(1/2). Maybe I messed up my time scheme. Anyway, I'll check my time integration again.

BTW, is there similar way like this to verify my spatial convergence?

[praveen]

What was the order of your error norm |u1-u2| when you performed time step convergence test ? If it is already close to machine zero, you will not be able to measure the time step convergence rate, since errors are dominated by round off errors.

When you perform spatial convergence test, you have to reduce time step also due to cfl condition.

If you get correct rates in this spatial convergence study, then you can take it that both your spatial and time schemes are fine.

[TurbJet]

Originally Posted by praveen

What was the order of your error norm |u1-u2| when you performed time step convergence test ? If it is already close to machine zero, you will not be able to measure the time step convergence rate, since errors are dominated by round off errors.

When you perform spatial convergence test, you have to reduce time step also due to cfl condition.

If you get correct rates in this spatial convergence study, then you can take it that both your spatial and time schemes are fine.

The norm is still away from machine precision, which ~1E-5. Maybe my time step size is too large. I'll try it again.

UPDATE: it's still not good, p ~ 1.

[FMDenaro]

I suggest some hints:


1) start using a smooth solution, that is like sin(x)
2) start checking the one-step error, that is compute the slope of the errors only after 1 time step. It must provide the correct accuracy if you have no bug
3) Start the convergence analysis using constant cfl for all the grids




post here the error curve you get

[TurbJet]

Originally Posted by FMDenaro

I suggest some hints:


1) start using a smooth solution, that is like sin(x)
2) start checking the one-step error, that is compute the slope of the errors only after 1 time step. It must provide the correct accuracy if you have no bug
3) Start the convergence analysis using constant cfl for all the grids




post here the error curve you get

Sure, I'll definitely try this and I'll get back to you.

However, I am not sure about your 2nd point: here you are talking about testing spatial convergence (with constant time step size and varying grid spacing), am I right?

[FMDenaro]

Originally Posted by TurbJet

Sure, I'll definitely try this and I'll get back to you.

However, I am not sure about your 2nd point: here you are talking about testing spatial convergence (with constant time step size and varying grid spacing), am I right?

No. Let me explain better this issue. The first thing to do is to assess that your discretization is consistent, without bugs. Now that means to assess that the local truncation error (LTE) vanishes as dt,h ->0. Using an exact solution, you actually are measuring a norm of the discretization error (DE), not of the LTE.
In a one-time analysis, you have O(DE)=dt*O(LTE), hence your LTE is of the order of DE divided by dt.
However, immagine the LTE as a function of dt and h. The consistence of the scheme requires you simultaneously make dt and h tending to zero. This can be obtained if you use the same cfl value on all the grids. That is dt and h must be diminished accordingly.
If this test will give you the correct error slope then you can check separately the spatial and temporal accuracy. That is, for the temporal accuracy, you can fix a very small value h (you need to "cancel" any effect of the spatial error) and, taking h fixed, you can vary the dt. This is a test at variable cfl.

[praveen]

There is another issue that some schemes will have larger errors if you use a very small cfl.

For linear advection equation, most schemes give exact solution with cfl=1. If you keep h fixed and reduce dt, then your cfl -> 0 and your errors can actually go up !!!

[TurbJet]

Originally Posted by FMDenaro

I suggest some hints:


1) start using a smooth solution, that is like sin(x)
2) start checking the one-step error, that is compute the slope of the errors only after 1 time step. It must provide the correct accuracy if you have no bug
3) Start the convergence analysis using constant cfl for all the grids




post here the error curve you get

Ok, I tried with one fixed extreme small time step size dt = 1E-6, and start with grid size of dx = 1E-6 and 1E-7; my convective velocity is exactly 1 and the function I am transporting is on [0 , 1].

I follow your advice by transporting it with only one time step, so after one time step, the the left point of the function should hit x = 1E-6. Thus I evaluate the value at x = 1E-6 which should be close to sin(0)=0.

But what I get from these two grid sizes are 1E-7 and 1E-8, only first order apparently. (I don't want to go to finer grid, the computation cost is way to high for my crapy laptop)

I also tried lowering dx and dt simultaneously by setting CFL=1 (hope this will not cause other problem), and I start with the same grid size as above. The results are pretty much the same.


Is it possible the error comes from the dispersion error, cause the spatial scheme I use for hyperbolic part is 4th-order central difference.

[FMDenaro]

Originally Posted by TurbJet

Ok, I tried with one fixed extreme small time step size 1E-6, and start with grid size of 1E-6 and 1E-7; my convective velocity is exactly 1 and the function I am transporting is sin(2pi x) on [0 , 1].

I follow your advice by transporting it with only one time step, so after one time step, the the left point of the function should hit x = 1E-6. Thus I evaluate the value at x = 1E-6 which should be close to sin(0)=0.

But what I get from these two grid sizes are 1E-7 and 1E-8, only first order apparently. (I don't want to go to finer grid, the computation cost is way to high for my crapy laptop)

I also tried lowering dx and dt simultaneously by setting CFL=1 (hope this works), and I start with the same grid size as above. The results are pretty much the same.


Is it possible the error comes from the dispersion error, cause the spatial scheme I use for hyperbolic part is 4th-order central difference.

I suggested the one-step analysis exactly to avoid the onset of any numerical instability issue.After 1 time step you must measure the correct convergence slope, irrespectly of oscillations.
Just to assess that you are computing exactly the error norm, try using the simple FTUS scheme (first order in time and space) and then the Lax-Wendroff scheme (second order time and space). If these schemes provide the correct convergence slope you have a bugs in your 4th order discretization.


PS: using the constant cfl you can start the analysis also using larger grid size, you should be able to get at least 5-6 values for the error curve. I suggest using the max norm so you can check the location of the max error

[TurbJet]

Originally Posted by FMDenaro

I suggested the one-step analysis exactly to avoid the onset of any numerical instability issue.After 1 time step you must measure the correct convergence slope, irrespectly of oscillations.
Just to assess that you are computing exactly the error norm, try using the simple FTUS scheme (first order in time and space) and then the Lax-Wendroff scheme (second order time and space). If these schemes provide the correct convergence slope you have a bugs in your 4th order discretization.


PS: using the constant cfl you can start the analysis also using larger grid size, you should be able to get at least 5-6 values for the error curve. I suggest using the max norm so you can check the location of the max error

Well I tried FTUS for one time step with dx=1E-3 and dt=1E-3, with CFL=1, the dispersion still happens at those regions where were originally zeros. So I can get a good convergence order as well.
(curiously, why CFL=1 has no significant diffusion? I thought upwind scheme is well-known for its diffusivity.)

Should I focus only on those region where the wave lies, and excluding those "zeros" regions, when computing the order?

I don't have time for trying Lax-Wendroff now. I'll give it a shot tomorrow and I'll come back and update this reply.

[FMDenaro]

Originally Posted by TurbJet

Well I tried FTUS for one time step with dx=1E-3 and dt=1E-3, with CFL=1, the transport is exact in one step
(curiously, why CFL=1 has no significant diffusion? I thought upwind scheme is well-known for its diffusivity.)

I don't have time for trying Lax-Wendroff now. I'll give it a shot tomorrow and I'll come back and update this reply.

At cfl=1 the FTUS gives the exact solution! You should check the accuracy at cfl<1, only this way you see the effect of the numerical diffusion.

Again, to assess that you have no bugs, for any scheme you test, it is required you obtain several values for the error slope. You cannot say nothing only by a couple of points.

[TurbJet]

Originally Posted by FMDenaro

At cfl=1 the FTUS gives the exact solution! You should check the accuracy at cfl<1, only this way you see the effect of the numerical diffusion.

Again, to assess that you have no bugs, for any scheme you test, it is required you obtain several values for the error slope. You cannot say nothing only by a couple of points.

I didn't get to test Lax-Wendroff. But for the upwind scheme, the convergence order for L_inf is almost zero (with fixed dt); The reason I think is due to the dispersion at these non-smooth regions (highlighted by red circle).

However, lowering dx/dt at the same time (where CFL = 0.9) can give the 1st order convergence in L_inf.

upwind.jpg