[petsc-users] An issue of order of accuracy.

[Barry Smith]

On Dec 10, 2013, at 10:33 AM, Alan Z. Wei <zhenglun.wei at gmail.com> wrote:

> Dear Dr. Smith,
> Indeed, non-uniform grid with the Poisson has only 1st order; my docx.
> file shows that also. However, it should recover to 'exact solution'
> with only small round-off error if the 'exact solution' is a polynomial
> function of degree 2 or less. In my tests, I constructed a quadratic
> function and the round-off error is quite large. I wonder if my usage of
> ksp solver or pc has some problems.

   Are you sure that it reproduces a quadratic exactly? With pencil and paper apply the differencing to a quadratic function what does it generate for the right hand side?

   Barry

> 
> thanks,
> Alan
>> 
>>  Alan,
>> 
>>    For non-uniform grid with the Poisson the order is no longer 2nd order. It is only for uniform grid that you get “magic” cancelation that makes it second order.
>> 
>>    Barry
>> 
>> On Dec 10, 2013, at 12:45 AM, Alan <zhenglun.wei at gmail.com> wrote:
>> 
>>> Thank you, Dr. Smith.
>>> I 'cooked up' a 4th-order polynomial and find out that the Poisson
>>> solver with uniform Cartesian mesh is 2nd order of accuracy. Also, if
>>> the solution is 3rd or less order polynomial, the L2Norm shows very
>>> small round-off error, which is usually 1e-7 if the -ksp_rtol is 1-e7.
>>> Following this idea, I re-derived the coefficients for the Poisson
>>> equation for solver with non-uniform grid, which is attached with a
>>> .docx file. From the equation, it indicates that at least 3rd or higher
>>> order polynomial should be constructed for the solution to detect the
>>> real order of accuracy. In other words, for solution of 2nd or lower
>>> order polynomials, the L2Norm should show very small round-off error.
>>> However, the reality is not and here is a brief description of my tests.
>>> The test focuses on the three-dimensional Poisson solver with
>>> non-uniform Cartesian mesh, which is modified from
>>> /src/ksp/ksp/example/tutorial/ex45.c. For simplicity, in this test, only
>>> x-direction has non-uniform grid, y- and z- are uniform grid. I have 3
>>> sets of tests, which are recorded in TestRun attached. The difference
>>> among these 3 sets of tests are the -ksp_rtol. Each set has 3 different
>>> commands to run the program with different ksp solver (GMRES or
>>> richardson) and pc type (GAMG and Hypre boomeramg).
>>> As I recorded in TestRun, while the ksp_rtol is 1e-7, the L2Norm is
>>> fairly large. As Dr. Smith explained, tight algebraic tolerances are
>>> needed to eliminate the algebraic error. Therefore, I changed ksp_rtol
>>> to 10-e15. However, the L2Norm can only reach 1e-7. Compared with the
>>> Poisson solver with uniform Cartesian mesh (which the L2Norm exhibits
>>> 1e-7 round-off error as -ksp_rtol = 1e-7), the L2Norm from the
>>> non-uniform Poisson solver is relatively high. Is this normal?
>>> Moreover, the Residual norm shown for cases with ksp_rtol = 1e-15 is
>>> around 7, 28 or even 81, which are far away from the real ksp_rtol
>>> imported. I monitored the ksp iterations with ksp_monitor and found out
>>> that these solvers usually iterate around 20 times. Why wouldn't them
>>> continue to iterate until the 1e-15 is achieved?
>>> Based on my observation, neither solver/pc works fine for this Poisson
>>> solver with non-uniform mesh. Is there any other option to made some
>>> improvements?
>>> At last, for the Poisson solver with non-uniform Cartesian grid, is
>>> there any better way to prove its validity?
>>> 
>>> sincerely appreciate,
>>> Alan
>>> 
>>>> Alan,
>>>> 
>>>>  I changed your initial grid size to 10 to get faster solve times and get what is below.  Note that your “exact solution” is quadratic, since the method is second order this means that not only is the “exact solution” an exact solution to the PDE, it is also an exact solution to the algebraic equations (for any grid size) hence the L2Norm of the error is only due to the round off of the algebraic solution, not due to any discretization error. In general, when testing for discretization error you always need to “cook up” an exact solution that is not completely represented in the approximation space. You also need to use a really tight algebraic tolerance to eliminate the algebraic error from the computation
>>>> 
>>>> 
>>>> 
>>>> ~/Src/petsc/test-dir  master $ ./ex45 -pc_type mg  -ksp_rtol 1e-12 
>>>> mx = 10, my = 10, mz =10, mm = 1, nn = 1, pp = 1
>>>> Residual norm 1.81616e-12
>>>> L2Norm = 1.107359e-12
>>>> Total Time Elapsed: 0.048599
>>>> ~/Src/petsc/test-dir  master $ ./ex45 -pc_type mg  -ksp_rtol 1e-12 -da_refine 1
>>>> mx = 19, my = 19, mz =19, mm = 1, nn = 1, pp = 1
>>>> Residual norm 3.36741e-12
>>>> L2Norm = 1.037148e-12
>>>> Total Time Elapsed: 0.183398
>>>> ~/Src/petsc/test-dir  master $ ./ex45 -pc_type mg  -ksp_rtol 1e-12 -da_refine 2
>>>> mx = 37, my = 37, mz =37, mm = 1, nn = 1, pp = 1
>>>> Residual norm 1.09476e-11
>>>> L2Norm = 2.330658e-12
>>>> Total Time Elapsed: 1.180839
>>>> ~/Src/petsc/test-dir  master $ ./ex45 -pc_type mg  -ksp_rtol 1e-12 -da_refine 3
>>>> mx = 73, my = 73, mz =73, mm = 1, nn = 1, pp = 1
>>>> Residual norm 3.19809e-11
>>>> L2Norm = 2.278763e-12
>>>> Total Time Elapsed: 10.819450
>>>> ~/Src/petsc/test-dir  master $ ./ex45 -pc_type mg   -da_refine 3
>>>> mx = 73, my = 73, mz =73, mm = 1, nn = 1, pp = 1
>>>> Residual norm 0.000103197
>>>> L2Norm = 1.011806e-05
>>>> Total Time Elapsed: 7.250106
>>>> 
>>>> 
>>>> On Dec 4, 2013, at 5:25 PM, Alan Z. Wei <zhenglun.wei at gmail.com> wrote:
>>>> 
>>>>> Dear all, 
>>>>>   I hope you had a great Thanksgiving. 
>>>>>   Currently, I tested the order of accuracy for /src/ksp/ksp/tutorial/example/ex45.c. Since the 2nd-order discretization is used in this program and ksp solver is converged to 10^-7, I expected that the solution should provides a 2nd-order in L2 norm. However, as I tested (even with a Laplace equation), the L2 norm slope is much less than 2. Sometime, if the grid size is reduced, the L2 norm increases. Could anyone help me about this issue, please?
>>>>> 
>>>>> Here is the L2 norm outputted:
>>>>> 
>>>>> Grid	L2 norm (10^-8)
>>>>> 0.05	4.36242
>>>>> 0.025	2.20794
>>>>> 0.0125	7.02749
>>>>> 0.00625	12.64   
>>>>>   Once the grid size is reduced to half, the number of the grid will be multiplied by 8 in order to keep the same size of the computational domain.
>>>>>   The code is also attached. It is from ex45.c with very little modifications. 
>>>>> 
>>>>> thanks in advance,
>>>>> Alan
>>>>> 
>>>>> <ex45.c><TestRun.txt>
>>> <ex45.c><makefile.txt><TestRun.txt><Non-Uniform Poisson.docx>
>