On Dec 10, 2013, at 10:33 AM, Alan Z. Wei <[email protected]> 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 <[email protected]> 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 <[email protected]> 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> >
