[lzhaok6]

I'm solving a 1D linear acoustic wave equation using Finite Element Method for spatial discretization (linear element) and second order Finite Difference Method (explicit Newmark scheme) for temporal discretization. Specifically, I want to propagate an exponential decay shock wave in a 1D shock tube (use boundary condition to imposing the wave profile). The governing equation and boundary condition are attached.

As we know, the dispersion is an inherent property of FEM when solving wave propagation problem and using CFL number equals to one seems to give the best solution since the numerical wave speed then equals the physical wave speed.

My question is:
How to examine the "observed order of accuracy" of temporal discretization in this problem? When I used CFL=1, there seems to be no dispersion in the solution and the wave profile after a certain time of propagation agrees well with the theoretical profile. Then, I tried to fix the mesh size and decrease the time step, which means CFL number is reduced and smaller than 1. However, right after the CFL number drops below 1, I start to see spurious oscillation right after the wave front which increases the overall error in the solution. Thus, the accuracy is even degraded when the time step size is decreased. As a result, I cannot retrieve the theoretical order of accuracy of time discretization.

I saw people recommend to keep CFL number fixed while investigating the time integration accuracy, which means mesh size is decreased as time step is decreased. Is this a reasonable way to examine the temporal order of accuracy? Did I misunderstand something?

[FMDenaro]

Theoretically, you cannot perform an accuracy test for a case where the solution is singular... the accuracy test has the aim to evaluate the order of the magnituted in the leading term appearing in the local truncation error. But if the derivative does not exist that is not correct.
The accuracy order analysis are performed with smooth solutions.

You can find some details in the book of Leveque.

[lzhaok6]

Originally Posted by FMDenaro

Theoretically, you cannot perform an accuracy test for a case where the solution is singular... the accuracy test has the aim to evaluate the order of the magnituted in the leading term appearing in the local truncation error. But if the derivative does not exist that is not correct.
The accuracy order analysis are performed with smooth solutions.

You can find some details in the book of Leveque.

Hi FMDenaro,

Thanks for your response. I don't actually understand why the solution is singular... By using 2nd order Newmark temporal discretization, I do have the truncation error term. Theoretically, if I fix the mesh size and decrease the time step, I should get the correct order of accuracy. However, the thing is the dispersion of spatial discretization seems to deteriorate the overall accuracy. As a result, I'm not able to observe the theoretically observed order of accuracy.

[FMDenaro]

You wrote:

"Specifically, I want to propagate an exponential decay shock wave in a 1D shock tube "

[lzhaok6]

Originally Posted by FMDenaro

You wrote:

"Specifically, I want to propagate an exponential decay shock wave in a 1D shock tube "

The shock wave pressure is added as the boundary condition on the left end of the shock tube.
I believe your point is that the problem is caused by the discontinuity of the boundary condition at the first time step. Is that correct?

[FMDenaro]

I don't understand ...you have a linear PDE describing a couple of travelling waves at +/-c. The initial/boundary conditions define the presence or not of the travelling discontinuity along x. Therefore, if you set a discontinuity it will travel and decay in time.
It's the same whe you work with the first order linear equation df/dt+c*df/dx=0 and set a discontinuity in condition for the function f. It travels along x and you cannot define there the derivative.

[FMDenaro]

Originally Posted by lzhaok6

The shock wave pressure is added as the boundary condition on the left end of the shock tube.
I believe your point is that the problem is caused by the discontinuity of the boundary condition at the first time step. Is that correct?

Why do not use an initial smooth function like a gaussian centered at x=0? you can just let to evolve towards +/- x direction. This way you can check the orde of accuracy in time and space.

[lzhaok6]

Originally Posted by FMDenaro

Why do not use an initial smooth function like a gaussian centered at x=0? you can just let to evolve towards +/- x direction. This way you can check the orde of accuracy in time and space.

Thank you very much for your suggestion. I've tried to impose a sin wave on the boundary. However, the observed order of accuracy is still not as expected. The L2 norm of error is like: 0.0032 when CFL=0.1; 0.0031 when CFL=0.2; 0.0055 when CFL=0.4; 0.0051 when CFL=0.8.

I think the overall accuracy is influenced by the dispersion of spatial discretization of finite element method. The dispersion would increase as the mismatch between the numerical speed and physical speed is enlarged (i.e. when CFL deviates from 1). However, I'm not sure if this explanation is theoretically sound.

[FMDenaro]

Originally Posted by lzhaok6

Thank you very much for your suggestion. I've tried to impose a sin wave on the boundary. However, the observed order of accuracy is still not as expected. The L2 norm of error is like: 0.0032 when CFL=0.1; 0.0031 when CFL=0.2; 0.0055 when CFL=0.4; 0.0051 when CFL=0.8.

I think the overall accuracy is influenced by the dispersion of spatial discretization of finite element method. The dispersion would increase as the mismatch between the numerical speed and physical speed is enlarged (i.e. when CFL deviates from 1). However, I'm not sure if this explanation is theoretically sound.

First, set the initial condition at x=0 and fix a domain quite long from -L to +L such that you can evolve the solution without entering BC.s.
Then, to check the accuracy order you have to take CFL=constant and use several refined grids. On each grid compute the L_inf norm at a fixed time. The slope of the error curve versus the grid size will give you the accuracy order.

Numerical dispersion and/or dissipation depend on the combined discretization in time and space. Dispersion is not caused by the fact that CFL is not 1.

[lzhaok6]

Originally Posted by FMDenaro

First, set the initial condition at x=0 and fix a domain quite long from -L to +L such that you can evolve the solution without entering BC.s.
Then, to check the accuracy order you have to take CFL=constant and use several refined grids. On each grid compute the L_inf norm at a fixed time. The slope of the error curve versus the grid size will give you the accuracy order.

Numerical dispersion and/or dissipation depend on the combined discretization in time and space. Dispersion is not caused by the fact that CFL is not 1.

Since I want to examine the order of accuracy of time integration, I suppose I should use different refinement level of time step and change the dx at the same time to keep the CFL constant. This would result in the same result as what you've suggested. I'll do the experiment and post what I observe in this thread.

Is there any particular kind of reason of fixing CFL number in order to examine the order of accurate of temporal/spatial discretization? (To keep the physics the same? )

[FMDenaro]

The order of the accuracy can be evaluated if all the terms of the local truncation error vanish simultaneously, therefore one performs a grid refinement by taking CFL=constant.

If you want check only the dependence on the time step, then you need to make the spatial dependence not relevant. This is obtained by using only one spatial grid with the mesh size very very small. Then, you can perform the analysis by refining the time step (that is you work by reducing the CFL). It is fundamental you set a very very small mesh size

[lzhaok6]

Originally Posted by FMDenaro

The order of the accuracy can be evaluated if all the terms of the local truncation error vanish simultaneously, therefore one performs a grid refinement by taking CFL=constant.

If you want check only the dependence on the time step, then you need to make the spatial dependence not relevant. This is obtained by using only one spatial grid with the mesh size very very small. Then, you can perform the analysis by refining the time step (that is you work by reducing the CFL). It is fundamental you set a very very small mesh size

Prof. Denaro,

I tried your suggestion. Since I want to examine the temporal order of accuracy, I have to keep my mesh size fixed. Seems like the dispersion problem is greatly relieved by using a Gaussian function as the waveform. The attached thumbnail is the observed order of accuracy I got from the L_inf norm. It seems like the order of accuracy is not perfectly 2nd order as I expected. But at least the solution is converging.

[FMDenaro]

First, I strongly suggest to represent the numerical solution only with symbols, not with continue lines. That will help the visualization.
How fine is the spatial grid?

Then, I do not understand the plot of the accuracy order. What are you representing? The accuracy order is estimated by the asymptotical slope of the curve obtained by plotting the error versus the dt in a double logarithic scale. What are the coordinates in your plot?

[lzhaok6]

Originally Posted by FMDenaro

First, I strongly suggest to represent the numerical solution only with symbols, not with continue lines. That will help the visualization.
How fine is the spatial grid?

Then, I do not understand the plot of the accuracy order. What are you representing? The accuracy order is estimated by the asymptotical slope of the curve obtained by plotting the error versus the dt in a double logarithic scale. What are the coordinates in your plot?

The spatial mesh size is 0.005m. The order of accuracy is calculated by ln(L(k+1)/L(k))/ln(r), in which L is the error norm, k is the refinement level, and r is the refinement rate (r=2 in this case). The X axis is the relative grid spacing in logarithmic scale (h=1 means the finest grid); Y axis is the order of accuracy calculated by the equation above. The result is actually the discrete order of accuracy without curve fitting.

I plotted the log-log plot of error norm vs grid spacing. It seems like the order of accuracy is between 1st and 2nd order.

[FMDenaro]

But for grid size you are denoting the time step, right? why the scale is 10 up to 10^2? If you have a second order spatial accuracy you get in the LTE a term O(10^-5). I suggest to use a more refined spatial grid.
Furthermore, the order of accuracy is reached only asymptotically, oscillations are quite usual before, so I suggest to produce yet solutions for smaller time steps.

[lzhaok6]

Originally Posted by FMDenaro

But for grid size you are denoting the time step, right? why the scale is 10 up to 10^2? If you have a second order spatial accuracy you get in the LTE a term O(10^-5). I suggest to use a more refined spatial grid.
Furthermore, the order of accuracy is reached only asymptotically, oscillations are quite usual before, so I suggest to produce yet solutions for smaller time steps.

Yes, the grid spacing stands for the time step. However, it is a relative time step. For instance, h=1 stands for dt=1e-5; h=2 stands for dt=2e-5; h=4 stands for dt=4e-5... Right now, I'm fixing spatial mesh size. So I suppose the ideal order of accuracy is still 2nd order.

The spatial discretization is only linear since I'm using a basic linear element.

Do you mean that the oscillation is a normal phenomenon? I'll try to use smaller time step.

[FMDenaro]

Originally Posted by lzhaok6

Yes, the grid spacing stands for the time step. However, it is a relative time step. For instance, h=1 stands for dt=1e-5; h=2 stands for dt=2e-5; h=4 stands for dt=4e-5... Right now, I'm fixing spatial mesh size. So I suppose the ideal order of accuracy is still 2nd order.

The spatial discretization is only linear since I'm using a basic linear element.

Do you mean that the oscillation is a normal phenomenon? I'll try to use smaller time step.

Use a finer spatial grid step (for example 10^-3) and check the slope at smaller time steps.
The convergence to the correct slope is never reached monotonically...

[lzhaok6]

Originally Posted by FMDenaro

Use a finer spatial grid step (for example 10^-3) and check the slope at smaller time steps.
The convergence to the correct slope is never reached monotonically...

I tried to use 4000 elements (0.0025m per element) and further decreased the time step size. It seems like the order of accuracy gradually converges to 1st order. Theoretically, the order of accuracy should be 2nd order. I think I might have misimplemented something in the scheme that finally degraded the precision. However, I believe there is another possibility that the numerical property of the problem is changed when CFL number is different, which makes the observed order of accuracy differ from the theoretical one.

[FMDenaro]

Originally Posted by lzhaok6

I tried to use 4000 elements (0.0025m per element) and further decreased the time step size. It seems like the order of accuracy gradually converges to 1st order. Theoretically, the order of accuracy should be 2nd order. I think I might have misimplemented something in the scheme that finally degraded the precision. However, I believe there is another possibility that the numerical property of the problem is changed when CFL number is different, which makes the observed order of accuracy differ from the theoretical one.

Yes, the slope seems to indicate also less than first order of accuracy.
The CFL (the ratio dt/dx) says the way in which you evaluate the local truncation error that is a function of the discretization steps. In the 2D space, you are fixing dx and you are evaluating the LTE for dt diminuishing.

You could also perform the one-step analysis to check the accuracy

[lzhaok6]

Thank you for your suggestions, Prof. Denaro. I'll try the one-step analysis then.