CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Time integration order of accuracy when solving acoustic wave equation using FEM

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 1, 2017, 15:35
Post Time integration order of accuracy when solving acoustic wave equation using FEM
  #1
New Member
 
Join Date: Jul 2016
Posts: 28
Rep Power: 10
lzhaok6 is on a distinguished road
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?
Attached Images
File Type: jpg Screen Shot 2017-05-01 at 2.26.53 PM.jpg (43.0 KB, 13 views)
lzhaok6 is offline   Reply With Quote

Old   May 1, 2017, 16:18
Default
  #2
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,893
Rep Power: 73
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
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.
juliom likes this.
FMDenaro is online now   Reply With Quote

Old   May 1, 2017, 17:15
Default
  #3
New Member
 
Join Date: Jul 2016
Posts: 28
Rep Power: 10
lzhaok6 is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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.
lzhaok6 is offline   Reply With Quote

Old   May 1, 2017, 17:32
Default
  #4
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,893
Rep Power: 73
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
You wrote:

"Specifically, I want to propagate an exponential decay shock wave in a 1D shock tube "
FMDenaro is online now   Reply With Quote

Old   May 1, 2017, 17:49
Default
  #5
New Member
 
Join Date: Jul 2016
Posts: 28
Rep Power: 10
lzhaok6 is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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?
lzhaok6 is offline   Reply With Quote

Old   May 1, 2017, 18:15
Default
  #6
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,893
Rep Power: 73
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
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 is online now   Reply With Quote

Old   May 1, 2017, 18:35
Default
  #7
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,893
Rep Power: 73
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by lzhaok6 View Post
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.
FMDenaro is online now   Reply With Quote

Old   May 1, 2017, 18:56
Default
  #8
New Member
 
Join Date: Jul 2016
Posts: 28
Rep Power: 10
lzhaok6 is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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.
lzhaok6 is offline   Reply With Quote

Old   May 1, 2017, 19:03
Default
  #9
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,893
Rep Power: 73
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by lzhaok6 View Post
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.
FMDenaro is online now   Reply With Quote

Old   May 1, 2017, 21:06
Default
  #10
New Member
 
Join Date: Jul 2016
Posts: 28
Rep Power: 10
lzhaok6 is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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? )

Last edited by lzhaok6; May 1, 2017 at 23:50.
lzhaok6 is offline   Reply With Quote

Old   May 2, 2017, 03:59
Default
  #11
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,893
Rep Power: 73
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
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
FMDenaro is online now   Reply With Quote

Old   May 2, 2017, 15:19
Default
  #12
New Member
 
Join Date: Jul 2016
Posts: 28
Rep Power: 10
lzhaok6 is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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.
lzhaok6 is offline   Reply With Quote

Old   May 2, 2017, 15:28
Default
  #13
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,893
Rep Power: 73
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
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?
FMDenaro is online now   Reply With Quote

Old   May 2, 2017, 17:28
Default
  #14
New Member
 
Join Date: Jul 2016
Posts: 28
Rep Power: 10
lzhaok6 is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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.
lzhaok6 is offline   Reply With Quote

Old   May 2, 2017, 17:46
Default
  #15
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,893
Rep Power: 73
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
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.
FMDenaro is online now   Reply With Quote

Old   May 2, 2017, 18:03
Default
  #16
New Member
 
Join Date: Jul 2016
Posts: 28
Rep Power: 10
lzhaok6 is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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.
lzhaok6 is offline   Reply With Quote

Old   May 2, 2017, 18:08
Default
  #17
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,893
Rep Power: 73
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by lzhaok6 View Post
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...
FMDenaro is online now   Reply With Quote

Old   May 3, 2017, 02:21
Default
  #18
New Member
 
Join Date: Jul 2016
Posts: 28
Rep Power: 10
lzhaok6 is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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.
lzhaok6 is offline   Reply With Quote

Old   May 3, 2017, 03:52
Default
  #19
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,893
Rep Power: 73
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by lzhaok6 View Post
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
FMDenaro is online now   Reply With Quote

Old   May 3, 2017, 10:31
Default
  #20
New Member
 
Join Date: Jul 2016
Posts: 28
Rep Power: 10
lzhaok6 is on a distinguished road
Thank you for your suggestions, Prof. Denaro. I'll try the one-step analysis then.
lzhaok6 is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
conjugate heat transfer in OpenFOAM skuznet OpenFOAM Running, Solving & CFD 99 March 16, 2017 06:07
High Courant Number @ icoFoam Artex85 OpenFOAM Running, Solving & CFD 11 February 16, 2017 14:40
Floating point exception error lpz_michele OpenFOAM Running, Solving & CFD 53 October 19, 2015 03:50
Cannot run the code properly: very large time step continuity error crst15 OpenFOAM Running, Solving & CFD 9 December 14, 2014 19:17
How to write k and epsilon before the abnormal end xiuying OpenFOAM Running, Solving & CFD 8 August 27, 2013 16:33


All times are GMT -4. The time now is 16:58.