- Email: [email protected]
Debugging Hints
Appendix B Debugging Hints The hints in this appendix are grouped per program, but assume that the reader has started by building (and debugging) tl~e sin~pler programs as M G l d . e and M G 2 d . e before attempting the more complicated ones. The authors strongly suggest that even if the readers objective is to build only tile E H L 2 d program, the reader should start to program and debug the M g l d and M G 2 d codes, until they perform as the examples given in chapter 3. Then the reader shouht use the M G 2 d code and modify it to obtain the H L 2 d code, debug it completely before moving on to the D R Y 2 d code and finally to the E H L 2 d code. Even if the reader looses a little time doing the extra programming, he or she will eventually save a lot of time debugging. Needless to say that the gradual approach ensures a much better understanding, which will be helpful when extending the programs to solve more complicated problems. The other basic rule which should be followed throughout the hints is that once dealing with a certain problem, the simplest possible configuration should be chosen, before trying the more complicated ones. Only if those sinlple configurations are fully debugged, the reader should move on and replace them one by one, by more complicated routines.
B.1
MGld.c & MG2d.c
9relax 0 Tile first task is to ensure that tile relaxation routine is correct, see Section 2.6, 3.9 and 3.10. Use man5," relaxation sweeps to solve the problem on level 1. Does tile residual norm diminish correctly (Figure 3.14 and 3,16)? If not check that tile residllal in tile point which has just been relaxed is zero. If not check if the residual is correctly computed, if necessary use the fact that the analytical solution is known. Does the solution obtained resemble the analytical solution? If the relaxation on level 1 works correctly check the convergence of relaxation on level 2, then 3. Convergence becomes very slow on finer grids, and many relaxation sweeps are necessary (Figure 3.14).
If relaxation works correctly oll all grids, it is time to start testing two level cycles. Start on level 2, relax twice, coarsen, relax ten times oil level 1, refine, relax on level 2. Has this cycle reduced ttle residual norm by a factor of 10? 255
A P P E N D I X B. DEBUGGING HINTS
256
If not, study the residual norm during the cycle, it should be continuous when going from level 2 to 1. Is it reduced by a factor of at least 100 on grid 17 A useful test is to modify coarsen_f0 slightly. Modify the coarse grid right hand side f by setting the fine grid residuals r~ to zero (3.74). The coarse grid solution u H should give zero coarse grid residuals, and the fine grid correction should also be zero. This test effectively decouples the fine grid relaxation from the coarse grid correction cycle, and deviations from the described behaviour point directly at the faulty routine. This test also permits a (graphical) comparison of the approximate solutions u on fine and coarse grid. It should be pointed out that a transfer error in a single point is sometimes sufficient to hinder convergence. In order to keep the transfer routines as simple as possible, start by using injection for u and r (Section 3.3). Once this works fine replace them by full weighting. Another useful tool is to start the 2 level cycle from a solution which has already converged on level 2, either through the use of many relaxation sweeps oil this grid, or through the use of the analytical solution. If the single cycle, 2 level program works correctly, it is time to study the many cycle behaviour. Each V-cycle should reduce the residual norm by a factor of 10, until the machine precision is reached (10-15). If tile residual norm does not decrease as predicted, this can be caused by some rather subtle errors, sometimes in a single point. Reread the points outlined above. If the 2 level many cycle behaviour is correct, it is time to study the 3 level one cycle behaviour, using the same techniques as described above for 2 levels. Then follows the level 3, many cycle test etc. When a program passes the level 2 many cycle test, the majority of the problems have been resolved. The last step is to introduce the FMG cycle, FMG interpolation and the approximate error norm calculation. The final test are performed using FMG with 1, 2 and 3 V-cycles, checking the discretization error and the numerical error (Section 3.10.3).
B.2
HL2d.c
Debugging the H L 2 d program uses the same steps as debugging M G 2 d . First it is necessary to check the proper functioning of the relaxation as was described previously. In addition, the correct implementation on different levels of the function H has to be checked. When studying the V-cycle behaviour it is advised to use the problem without cavitation first, and to establish tile correct convergence of this linear problem. Tile results should be compared with those obtained in Section 4.5. Because of the long and narrow domain, attention has to be paid on the coarsest grid with respect to the FMG high order interpolation. Tile fourth order interpolation requires at least four points in every direction. When adding the cavitation condition, it is very important to maintain a boundary which does not move between the different grids. This excludes the full weighting operators for u and r (Section 3.3), at least in the vicinity of the cavitation boundary. Furthermore, the prerelax 0 routine can be suppressed initially and replaced by relax() itself. Once good convergence is established, a significant amount of work can be gained
B.3. DRY2D.C
257
by using prerelaxO, but for debugging the simplicity of using a single routine is important. Finally the cavitating results should be comparable to the ones of Section 4.5.2.
B.3
DRY2d.c
Debugging the D R Y 2 d program requires a dual approach. The fast integration and the fast solver should be debugged separately. The order is not very important. To debug the fast integration start by checking the classical integration. When using a Hertzian pressure distribution, does one obtain the film thickness profile with nearly zero values for X 2 + )-2 < 17 If so is this also the case on levels 2 and 3? Does the diseretization error reduce as predicted? The second step is the implementation of the fast integration itself, starting with the second order and correcting over the entire domain rnl = m2 = 1000. The transfer routines in x and y direction should be programmed only once, using 'cut and paste' to obtain the second routine. Calculating the integral on level 2 using a single level deep, should give exactly the same result as the classical integration (up to the level of the machine precision). Large errors in all points normally stem from programming errors in tillK10 or domain errors in the summations. Correct values in even points but large errors in odd points hint at errors in the interpolation and correction routines. Finally, fine grid integrals with many levels deep should be computed, and should again give the same result as the classical integrals up to the level of the machine precision. If all these tests are correct, the small values of m l and m2 should be used, but this results rapidly in relatively large errors. A next step involves programming the fourth order routines. This involves the use of the 'ghost' points and req~lires a careful analysis of the summation boundaries. The steps described above should again be followed, paying attention to errors that tend to creep in from the boundaries when using more levels deep. This type of error hints at incorrect treatment of the 'ghost' points or errors in the summation boundaries. Finally, the small values of m l and m2 should be implemented, resulting in larger errors, but already an appreciable gain in computing speed is obtained. The sixth order routines, and finally the eight order routines follow the same rules, where the 'cut and paste' technique really pays of, in limiting the programming and debugging effort. When the fast integration work correctly, it is time to program and debug the fast solver, starting once again with the relaxation and the classical integration. Since this routine uses the distributed relaxation, it is important to test if tile five point relaxation does indeed result in zero residuals. On finer levels the relaxation becomes unstable when starting from a zero initial solution. Thus the debugging should be carried out starting from the Hertzian pressure distribution. In a first time, the H0 constant should be frozen at - 1 and not be changed. When this results in a correct convergence, the Ho value can be varied to account for the force balance equation, where overall convergence should only be slightly influenced by these changes in H0. When the fast solver works correctly, it can be combined with the fast integration, possibly using tile large values of ml - m2 = 1000 to obtain exactly the same results as
258
APPENDIX B. DEBUGGING HINTS
with the classical integration. Small differences can thus be easily detected and traced back to the defective routine. Furthern~ore, this allows a direct comparison between the second, fourth, sixth and eight order fast integration. When all errors are removed, the small values of m l and m2 should be introduced again, and the gain in computing time can be evaluated.
B.4
EHL2d.c
The main change from the DRY2d code to the EHL2d code is in the equations and the relaxation. In that case debugging the EHL2d.c code starts with debugging the relaxation. First construct the case where only Gauss-Seitlel line relaxation is used. This leads to the simplest system of equations. For a low load case, the relaxation should converge to machine error residuals. If the relaxation diverges, check if the system of equations is solved exactly. If the residual stalls it may be due to an error in the system of equations near the cavitation boundary. When the Gauss-Seidel line relaxation on a single grid works correctly one can proceed to the mixed Gauss-Seidel/Jacobi relaxation. Again this mixed relaxation is debugged by checking single grid convergence. Apply the relaxation on asingle grid to the model cases 3 1 = 20, L - 10 and ] l I = 200, L = 10. Compare the performance and the result with the figures given in Chapter 6. Also values of the minimum and central fihn thickness can be compared with values given in Figure 6.15. Is the value roughly correct ? If the single grid performance is satisfactory one can proceed to the multigrid case. First check if a cycle converges for 2 levels, many cycles, then for 3 levels, many cycles, etc. Compare the results with the output given in Chapter 6 for the model cases. If the cycle does not converge check the coarse grid correction by not transferring residuals of the Reynolds equation in the routine coarscn_f. In that case the residuals on all coarser grids should be zero. and the solver should act as if it performs relaxations on a single grid. Naturally for this test H0 should be fixed. If one has obtained a code that can reproduce to a good extend the results presented in Chapter 6 one can start to apply it to other cases.
B.1
MGld.c & MG2d.c
9relax 0 Tile first task is to ensure that tile relaxation routine is correct, see Section 2.6, 3.9 and 3.10. Use man5," relaxation sweeps to solve the problem on level 1. Does tile residual norm diminish correctly (Figure 3.14 and 3,16)? If not check that tile residllal in tile point which has just been relaxed is zero. If not check if the residual is correctly computed, if necessary use the fact that the analytical solution is known. Does the solution obtained resemble the analytical solution? If the relaxation on level 1 works correctly check the convergence of relaxation on level 2, then 3. Convergence becomes very slow on finer grids, and many relaxation sweeps are necessary (Figure 3.14).
If relaxation works correctly oll all grids, it is time to start testing two level cycles. Start on level 2, relax twice, coarsen, relax ten times oil level 1, refine, relax on level 2. Has this cycle reduced ttle residual norm by a factor of 10? 255
A P P E N D I X B. DEBUGGING HINTS
256
If not, study the residual norm during the cycle, it should be continuous when going from level 2 to 1. Is it reduced by a factor of at least 100 on grid 17 A useful test is to modify coarsen_f0 slightly. Modify the coarse grid right hand side f by setting the fine grid residuals r~ to zero (3.74). The coarse grid solution u H should give zero coarse grid residuals, and the fine grid correction should also be zero. This test effectively decouples the fine grid relaxation from the coarse grid correction cycle, and deviations from the described behaviour point directly at the faulty routine. This test also permits a (graphical) comparison of the approximate solutions u on fine and coarse grid. It should be pointed out that a transfer error in a single point is sometimes sufficient to hinder convergence. In order to keep the transfer routines as simple as possible, start by using injection for u and r (Section 3.3). Once this works fine replace them by full weighting. Another useful tool is to start the 2 level cycle from a solution which has already converged on level 2, either through the use of many relaxation sweeps oil this grid, or through the use of the analytical solution. If the single cycle, 2 level program works correctly, it is time to study the many cycle behaviour. Each V-cycle should reduce the residual norm by a factor of 10, until the machine precision is reached (10-15). If tile residual norm does not decrease as predicted, this can be caused by some rather subtle errors, sometimes in a single point. Reread the points outlined above. If the 2 level many cycle behaviour is correct, it is time to study the 3 level one cycle behaviour, using the same techniques as described above for 2 levels. Then follows the level 3, many cycle test etc. When a program passes the level 2 many cycle test, the majority of the problems have been resolved. The last step is to introduce the FMG cycle, FMG interpolation and the approximate error norm calculation. The final test are performed using FMG with 1, 2 and 3 V-cycles, checking the discretization error and the numerical error (Section 3.10.3).
B.2
HL2d.c
Debugging the H L 2 d program uses the same steps as debugging M G 2 d . First it is necessary to check the proper functioning of the relaxation as was described previously. In addition, the correct implementation on different levels of the function H has to be checked. When studying the V-cycle behaviour it is advised to use the problem without cavitation first, and to establish tile correct convergence of this linear problem. Tile results should be compared with those obtained in Section 4.5. Because of the long and narrow domain, attention has to be paid on the coarsest grid with respect to the FMG high order interpolation. Tile fourth order interpolation requires at least four points in every direction. When adding the cavitation condition, it is very important to maintain a boundary which does not move between the different grids. This excludes the full weighting operators for u and r (Section 3.3), at least in the vicinity of the cavitation boundary. Furthermore, the prerelax 0 routine can be suppressed initially and replaced by relax() itself. Once good convergence is established, a significant amount of work can be gained
B.3. DRY2D.C
257
by using prerelaxO, but for debugging the simplicity of using a single routine is important. Finally the cavitating results should be comparable to the ones of Section 4.5.2.
B.3
DRY2d.c
Debugging the D R Y 2 d program requires a dual approach. The fast integration and the fast solver should be debugged separately. The order is not very important. To debug the fast integration start by checking the classical integration. When using a Hertzian pressure distribution, does one obtain the film thickness profile with nearly zero values for X 2 + )-2 < 17 If so is this also the case on levels 2 and 3? Does the diseretization error reduce as predicted? The second step is the implementation of the fast integration itself, starting with the second order and correcting over the entire domain rnl = m2 = 1000. The transfer routines in x and y direction should be programmed only once, using 'cut and paste' to obtain the second routine. Calculating the integral on level 2 using a single level deep, should give exactly the same result as the classical integration (up to the level of the machine precision). Large errors in all points normally stem from programming errors in tillK10 or domain errors in the summations. Correct values in even points but large errors in odd points hint at errors in the interpolation and correction routines. Finally, fine grid integrals with many levels deep should be computed, and should again give the same result as the classical integrals up to the level of the machine precision. If all these tests are correct, the small values of m l and m2 should be used, but this results rapidly in relatively large errors. A next step involves programming the fourth order routines. This involves the use of the 'ghost' points and req~lires a careful analysis of the summation boundaries. The steps described above should again be followed, paying attention to errors that tend to creep in from the boundaries when using more levels deep. This type of error hints at incorrect treatment of the 'ghost' points or errors in the summation boundaries. Finally, the small values of m l and m2 should be implemented, resulting in larger errors, but already an appreciable gain in computing speed is obtained. The sixth order routines, and finally the eight order routines follow the same rules, where the 'cut and paste' technique really pays of, in limiting the programming and debugging effort. When the fast integration work correctly, it is time to program and debug the fast solver, starting once again with the relaxation and the classical integration. Since this routine uses the distributed relaxation, it is important to test if tile five point relaxation does indeed result in zero residuals. On finer levels the relaxation becomes unstable when starting from a zero initial solution. Thus the debugging should be carried out starting from the Hertzian pressure distribution. In a first time, the H0 constant should be frozen at - 1 and not be changed. When this results in a correct convergence, the Ho value can be varied to account for the force balance equation, where overall convergence should only be slightly influenced by these changes in H0. When the fast solver works correctly, it can be combined with the fast integration, possibly using tile large values of ml - m2 = 1000 to obtain exactly the same results as
258
APPENDIX B. DEBUGGING HINTS
with the classical integration. Small differences can thus be easily detected and traced back to the defective routine. Furthern~ore, this allows a direct comparison between the second, fourth, sixth and eight order fast integration. When all errors are removed, the small values of m l and m2 should be introduced again, and the gain in computing time can be evaluated.
B.4
EHL2d.c
The main change from the DRY2d code to the EHL2d code is in the equations and the relaxation. In that case debugging the EHL2d.c code starts with debugging the relaxation. First construct the case where only Gauss-Seitlel line relaxation is used. This leads to the simplest system of equations. For a low load case, the relaxation should converge to machine error residuals. If the relaxation diverges, check if the system of equations is solved exactly. If the residual stalls it may be due to an error in the system of equations near the cavitation boundary. When the Gauss-Seidel line relaxation on a single grid works correctly one can proceed to the mixed Gauss-Seidel/Jacobi relaxation. Again this mixed relaxation is debugged by checking single grid convergence. Apply the relaxation on asingle grid to the model cases 3 1 = 20, L - 10 and ] l I = 200, L = 10. Compare the performance and the result with the figures given in Chapter 6. Also values of the minimum and central fihn thickness can be compared with values given in Figure 6.15. Is the value roughly correct ? If the single grid performance is satisfactory one can proceed to the multigrid case. First check if a cycle converges for 2 levels, many cycles, then for 3 levels, many cycles, etc. Compare the results with the output given in Chapter 6 for the model cases. If the cycle does not converge check the coarse grid correction by not transferring residuals of the Reynolds equation in the routine coarscn_f. In that case the residuals on all coarser grids should be zero. and the solver should act as if it performs relaxations on a single grid. Naturally for this test H0 should be fixed. If one has obtained a code that can reproduce to a good extend the results presented in Chapter 6 one can start to apply it to other cases.