[Simbelmynë]

Hey,

I am running a gravity driven flow between two plates in order to test the order of the numerical scheme. Does it matter what initial condition I start with? I guess my question is about how many iterations/time-steps that are needed in order to judge the convergence behavior of a method. If the initial condition is the true solution then it converges to machine precision within one time-step and never changes (the time derivatives are zero for all times, which seems correct). If I use an initial condition far from the true solution then it takes a relatively long time. Right now I set a condition close to the true solution and let it run for approximately 20000 time-steps on both the coarse and fine grid. Is this OK?

Second question is about local and global order. Which type is more common to report?

Regards

[FMDenaro]

Originally Posted by Simbelmynë

Hey,

I am running a gravity driven flow between two plates in order to test the order of the numerical scheme. Does it matter what initial condition I start with? I guess my question is about how many iterations/time-steps that are needed in order to judge the convergence behavior of a method. If the initial condition is the true solution then it converges to machine precision within one time-step and never changes (the time derivatives are zero for all times, which seems correct). If I use an initial condition far from the true solution then it takes a relatively long time. Right now I set a condition close to the true solution and let it run for approximately 20000 time-steps on both the coarse and fine grid. Is this OK?

Second question is about local and global order. Which type is more common to report?

Regards

Despite of the initial condition, if your exact solution is steady then you must ensure that all the time derivatives are smaller than a threshold low value.
The time required for reaching the steady state depends both on the initial condition and on the physics of the problem (The Re number is an estimation of the order of magnitute of the steady time).
I suggest using the measure of the error in a Linf norm.

[Simbelmynë]

Originally Posted by FMDenaro

Despite of the initial condition, if your exact solution is steady then you must ensure that all the time derivatives are smaller than a threshold low value.
The time required for reaching the steady state depends both on the initial condition and on the physics of the problem (The Re number is an estimation of the order of magnitute of the steady time).
I suggest using the measure of the error in a Linf norm.

I have used the average L1-norm of the time-derivative(s) as a measure of convergence. Not sure if it is correct to average it, but it seems odd to define a convergence criteria that is dependent on the grid size, so I used averaging. I set the criteria to 1e-8 in this case.

If I set the exact value as initial condition then the average L1-norm becomes exactly zero, which seems a bit odd considering machine round-off errors. Anyways, in this case the solution does not change at all with time, so I guess this would constitute a true numerical steady-state.

[FMDenaro]

Originally Posted by Simbelmynë

I have used the average L1-norm of the time-derivative(s) as a measure of convergence. Not sure if it is correct to average it, but it seems odd to define a convergence criteria that is dependent on the grid size, so I used averaging. I set the criteria to 1e-8 in this case.

If I set the exact value as initial condition then the average L1-norm becomes exactly zero, which seems a bit odd considering machine round-off errors. Anyways, in this case the solution does not change at all with time, so I guess this would constitute a true numerical steady-state.

Starting with the exact solution as initial condition is a sort of "consistence" test and is a good fact you get an unchanged solution with zero time-derivative.
However, I think that the Linf norm is better suited because it is also local, can give you indication of where possible errors are greater.

[Simbelmynë]

Originally Posted by FMDenaro

Starting with the exact solution as initial condition is a sort of "consistence" test and is a good fact you get an unchanged solution with zero time-derivative.
However, I think that the Linf norm is better suited because it is also local, can give you indication of where possible errors are greater.

Thank you I will try the Linf-norm. I also just tried to lower the criteria from 1e-8 to 1e-20 and now even the very coarse grid (33x33) converges to the true solution, even if I start from a zero initial field. How can I possibly determine the order of the method in this case? Perhaps gravity driven flow between plates is a poor method to use for determining the order of a method?

I use Fortran (with latest gfortran compiler) and have all reals at double precision.

[FMDenaro]

Originally Posted by Simbelmynë

Thank you I will try the Linf-norm. I also just tried to lower the criteria from 1e-8 to 1e-20 and now even the very coarse grid (33x33) converges to the true solution, even if I start from a zero initial field. How can I possibly determine the order of the method in this case? Perhaps gravity driven flow between plates is a poor method to use for determining the order of a method?

I use Fortran (with latest gfortran compiler) and have all reals at double precision.

usually, analytical solutions are best suited to check the convergence order...
5-6 grids should be sufficient to reache a monotonic slope in the error curve

[Simbelmynë]

Originally Posted by FMDenaro

usually, analytical solutions are best suited to check the convergence order...
5-6 grids should be sufficient to reache a monotonic slope in the error curve

I use the analytical solution. I am not sure how I can check the order of the method. It seems the solution is actually deteriorating with an increasing grid size, perhaps this is connected to an accumulating round-off error?

Example (analytical value is 0.03):

17x17 - 0.029999999999997
33x33 - 0.029999999999989
65x65 - 0.029999999999955

If I use the following expression to calculate the order,

p=log_2(abs((u_c-u_a)/(u_f-u_a)))

where u_c and u_f is the coarse and fine solutions respectively and u_a is the analytical solution, then I get negative orders at around -2...

So the solution is getting worse at a second order rate as I refine the grid?

[FMDenaro]

I strongly suggest:

1) use the time step of the finest grid for all the other grids
2) be careful in using the correct convergence criteria for the iterative solvers
3) for each grid compute the Linf error and then plot the errors versus the grid size in a double logarithmic scale

[Simbelmynë]

Originally Posted by FMDenaro

I strongly suggest:

1) use the time step of the finest grid for all the other grids
2) be careful in using the correct convergence criteria for the iterative solvers
3) for each grid compute the Linf error and then plot the errors versus the grid size in a double logarithmic scale

Thank you for some good suggestions!

1) This was a good suggestion. Now I have the following results

17x17 - 0.029999999999942
33x33 - 0.029999999999944
65x65 - 0.029999999999945

2) Seems OK. Tried lowering the criteria for the iterative solver by 3 orders of magnitude and the result does not change at all.

[FMDenaro]

Originally Posted by Simbelmynë

Thank you for some good suggestions!

1) This was a good suggestion. Now I have the following results

17x17 - 0.029999999999942
33x33 - 0.029999999999944
65x65 - 0.029999999999945

2) Seems OK. Tried lowering the criteria for the iterative solver by 3 orders of magnitude and the result does not change at all.

you have crearly a problem in the code ...check where the max error is located

[Simbelmynë]

Originally Posted by FMDenaro

you have crearly a problem in the code ...check where the max error is located

Here is a plot of the error, abs(u-u_a). It is only one line of 17 but I think it gives a good indication where errors are largest in the domain. Could you elaborate why there is clearly a problem in the code?

[FMDenaro]

what you should plot is the max(err) vs. the size h obtained on each grid

[Simbelmynë]

Originally Posted by FMDenaro

what you should plot is the max(err) vs. the size h obtained on each grid

This is basically what I have done (the maximum error is in the center and all values I have reported is from the center), but I put it in a log log plot here. h2 is a second order slope.

Is it common to compare errors that are so close to machine precision?

[FMDenaro]

Originally Posted by Simbelmynë

This is basically what I have done (the maximum error is in the center and all values I have reported is from the center), but I put it in a log log plot here. h2 is a second order slope.

Is it common to compare errors that are so close to machine precision?

No, you can see that your error does not scale at all...
this is why I wrote you have something wrong in the test you are performing.

[Simbelmynë]

Originally Posted by FMDenaro

No, you can see that your error does not scale at all...
this is why I wrote you have something wrong in the test you are performing.

OK. The test case is gravity driven flow between two plates. Boundary conditions are walls and periodic. the 17x17 grid converges to a maximum error of 1e-13 and the code behaves very well in the other test cases I have run (cavity flow, backward facing step). If I check the order of the code in a more complex flow it is approximately 2, but then I use the GCI method and 3 different grid sizes.

1) Do you think the test case here is viable to find out the order of the method?
2) If not, what method do you suggest?

Thanks.

[mprinkey]

The solution to the continuous problem is parabolic in space, right? If your spatial discretization algorithm can perfectly capture quadratic behavior, the "error scaling" is independent of the number of samples. That may be what you are seeing. Consider the "error" associated with linear interpolation across a LINEAR data set. It doesn't matter what two data points you chose, they are all collinear and give the same interpolation results (within round-off) independent of deltaX. Or another way to look at it, what is the Taylor series expansion of F(x) = x^2? Above 2nd order, they are all zero. There are no higher-order terms.

If you want to do error scaling tests, you need to have underlying solution behaviors that will produce H.O.T.s in discretization scheme so that you can measure them. You may just not have that in such a simple case. It would show that your scheme is consistent UP TO quadratic interpolation, but nothing about higher orders.

[FMDenaro]

I agree, if the exact solution is just a linear stable temperature profile, then the test is not suitable to be adopted to check the order of accuracy of the scheme. What about the velocity field, is it vanishing or has a constant streamwise value??

[Simbelmynë]

Thank you for your input!

Originally Posted by mprinkey

The solution to the continuous problem is parabolic in space, right? If your spatial discretization algorithm can perfectly capture quadratic behavior, the "error scaling" is independent of the number of samples. That may be what you are seeing. Consider the "error" associated with linear interpolation across a LINEAR data set. It doesn't matter what two data points you chose, they are all collinear and give the same interpolation results (within round-off) independent of deltaX. Or another way to look at it, what is the Taylor series expansion of F(x) = x^2? Above 2nd order, they are all zero. There are no higher-order terms.

If you want to do error scaling tests, you need to have underlying solution behaviors that will produce H.O.T.s in discretization scheme so that you can measure them. You may just not have that in such a simple case. It would show that your scheme is consistent UP TO quadratic interpolation, but nothing about higher orders.

This sounds reasonable. Is this common for FV/FD methods? If so then I guess the test problem is not so useful. Would a backward facing step against experimental data serve as a better test in your oppinion?

Originally Posted by FMDenaro

I agree, if the exact solution is just a linear stable temperature profile, then the test is not suitable to be adopted to check the order of accuracy of the scheme. What about the velocity field, is it vanishing or has a constant streamwise value??

I do not solve for temperature at all. I'm not sure I understand the question. I solve for the parabolic velocity profile between two plates (2d incompressible NS). It is constant in the streamwise direction, i.e. the parabolic profile is the same over the length of the domain.

[FMDenaro]

Originally Posted by Simbelmynë

Hey,

I am running a gravity driven flow between two plates in order to test the order of the numerical scheme. Does it matter what initial condition I start with? I guess my question is about how many iterations/time-steps that are needed in order to judge the convergence behavior of a method. If the initial condition is the true solution then it converges to machine precision within one time-step and never changes (the time derivatives are zero for all times, which seems correct). If I use an initial condition far from the true solution then it takes a relatively long time. Right now I set a condition close to the true solution and let it run for approximately 20000 time-steps on both the coarse and fine grid. Is this OK?

Second question is about local and global order. Which type is more common to report?

Regards

then what do you mean for gravity driven? I was thinking about a buoyancy flow....

hovewer, if you are running the simple Poiseuille test case, then the second derivative of u along the height is not zero and you should be able to check for the convergence of the error curve...
are you using Re=1?

[Simbelmynë]

Originally Posted by FMDenaro

then what do you mean for gravity driven? I was thinking about a buoyancy flow....

Well, I compared against an analytical solution... Anyways, you are correct the test case is Poiseuille flow (although not in cylindrical coordinates).

Originally Posted by FMDenaro

hovewer, if you are running the simple Poiseuille test case, then the second derivative of u along the height is not zero and you should be able to check for the convergence of the error curve...
are you using Re=1?

Re=20 in the above simulations.