7.5 A new mixed–variable iterative procedure
7.5.2 Examples
In this section the performance of the new iterative procedure is compared to the con-ventional mixed–ψ procedure. The example involves the infiltration of water into a one dimensional 60 cm column of soil with van Genuchten parameters ks = 0.00922 cm/s, θr = 0.102, θs = 0.368, α = 0.0355 cm−1, n = 2.0 and m = 1−1/n= 0.5. These pa-rameters, which correspond to what is shown in Figure 7.3, have been used by numerous authors in evaluating the performance of numerical procedures, see e.g. [4, 38, 31]. In all cases the full Newton method is applied, i.e. the exact Jacobian is used, and in all tests the time step is taken to vary as
∆tn+1= 2p
tn∆t0+ ∆t0 (7.33)
such that the entire time step selection procedure is controlled by the initial time step
∆t0. This time step selection results in the front propagating at approximately the same
Department of Civil Engineering - Technical University of Denmark 65
Numerical solution procedures 7.5 A new mixed–variable iterative procedure
Initial point: θθθ0n+1=θθθn,ψψψ0n+1=ψψψ0n Iterations: j= 0, . . . , jmax
Residual and Jacobian: rjn+1,Jjn+1 Increment: ∆ψψψjn+1=−(Jjn+1)−1rjn+1 Prediction: θeθθj+1n+1=θθθjn+1+
·∂θθθ
∂ψψψ
¸
∆ψψψjn+1
Variable check: k= 1, . . . ,number of nodes if(θeθθj+1n+1)k < θTOL
(θj+1n+1)k = (θej+1n+1)k
(ψj+1n+1)k=ψ((θej+1n+1)k) else
(ψj+1n+1)k= (ψn+1j )k+ (∆ψjn+1)k
(θj+1n+1)k =θ((ψen+1j+1)k) end variable check
end iteration
Table 7.2Iterative procedure for the Richards equation
discrete speed throughout the time stepping, i.e. for arbitrary times tn and tn+1, and corresponding front positionszn+1 andzn, the discrete velocity of the front is given by
Vf=∆zn+1
∆tn+1
=zn+1−zn tn+1−tn
'const., n= 0, . . . , nmax (7.34) whereVfdepends only on ∆t0.
In all the examples the switch between the two methods of update is made when the predicted effective saturation reaches a value close to unity, i.e. when
See= eθ−θr θs−θr
> STOLe (7.35)
where a value of STOLe = 0.9995 has been used in all examples. Furthermore, in all examples we have assumed vanishing compressibility, i.e. c0= 0.
Example 1
In the first example constant pressure head boundary conditions are imposed as ψtop =
−75 cm andψbottom=−1000 cm and the initial pressure head isψinit=−1000 cm
through-7.5 A new mixed–variable iterative procedure Numerical solution procedures
Figure 7.4Iterative procedure for the Richards equation
out the domain as in [4, 38, 31]. The efficiency of the method under these conditions is shown in Table 7.3, and the moisture content and pressure head profiles illustrated in Figure 7.5. As can be seen the new method clearly outperforms the conventional mixed–
ψ method. Not only are the iteration counts lower, but whereas the mixed–ψ method fails at large values of ∆t0the new method seems to handle this without any significant problems.
Mixed–ψ New method
∆t0(s) No. Step No. Iter Iter/Step No. Iter Iter/Step
1000.1 1000 2979 332.99 2010 2.01
1001.0 1317 1264 333.99 2947 2.99
1010.0 1100 1698 336.98 2404 4.04
1050.0 1045 1463 32.5 2199 4.42
1100.0 1032 1941 29.4 2165 5.12
1500.0 1015 † † 2097 6.47
1000.0 1010 † † 2070 7.00
Table 7.3Comparison of conventional mixed–ψscheme and new iterative procedure. Ex-ample 1,ψinit=−1000cm,ψtop=−75cm andψbottom=−1000cm. †Convergence not achieved after 200 iterations.
Department of Civil Engineering - Technical University of Denmark 67
Numerical solution procedures 7.5 A new mixed–variable iterative procedure
0 0.1 0.2 0.3
0 10 20 30 40 50 60
Effective saturation
Distance from bottom (cm)
t = 1,445 s
t = 8,820 s
t = 22,445 s
t = 42,320 s
t = 68,445 s
t = 100,820 s
(a)
-1000 -800 -600 -400 -200 0
10 20 30 40 50 60
Pressure head (cm)
Distance from bottom (cm)
t = 1,445 s t = 8,820 s t = 22,445 s t = 42,320 s
t = 68,445 s t = 100,820 s
(b)
Figure 7.5Moisture content (a) and pressure head (b) profiles. Example 1,ψtop=−75cm andψbottom=−1000cm,∆t0= 5.0s.
Example 2
Next a somewhat more challenging example is considered. The material and geometry is identical to that of the previous example, but the bottom is now a no–flow boundary, qbottom= 0 and at the top a pressure ofψtop= +100 cm, i.e. the infiltration is carried out with a constant free water level of 100 cm at the top of the column.
The solution statistics are given in Table 7.4. Compared to the previous example, some-what higher iteration counts are seen. In Figure 7.6 the reason for this is shown: because of the high pressure at the top, the front are much sharper and furthermore, the varia-tions in saturation within the column is much greater. Again the mixed–ψ method has problems with convergence for large time steps, although quite surprisingly, it does con-verge for the two largest values of ∆t0. On the other hand, the new method converges regardless of the time step, and although the number of iterations per time step increases as the time step increases, the total number of iterations decreases.
Example 3
As a final example we now consider the case where the soil is initially very dry corre-sponding toψinit=−10,000 cm, and otherwise all other parameters are as in the previous example. Such conditions are known to cause severe problems with the convergence of the mixed–ψmethod [19]. The results shown in Table 7.5 confirm this. The new method, however, is largely unaffected as compared to the previous example.
Together, these three examples clearly demonstrated the superiority of the new method over the conventional mixed–ψ method. The comparisons made are perhaps not entirely fair in that an automatic adjustment of the time step could surely bring the mixed–ψ
7.5 A new mixed–variable iterative procedure Numerical solution procedures
Mixed–ψ New method
∆t0(s) No. Step No. Iter Iter/Step No. Iter Iter/Step
10.0005 775 3207 304.14 2877 233.66
10.0010 548 2630 304.80 2192 234.00
10.0100 174 † † 2929 235.34
10.1000 755 † † 2447 238.12
11.0000 718 3321 17.8 2219 12.2
10.0000 776 3180 30.0 2140 23.3
Table 7.4Comparison of conventional mixed–ψscheme and new iterative procedure. Ex-ample 2, ψinit = −1000cm, ψtop = 100cm and qbottom = 0. †Convergence not achieved after 200 iterations.
0 0.2 0.4 0.6 0.8 1 0
10 20 30 40 50 60
Effective saturation
Distance from bottom (cm)
t = 5.76 s t = 29.16 s t = 70.56 s t = 129.96 s t = 207.36 s t = 302.76 s
(a)
-1000 -750 -500 -250 0 0
10 20 30 40 50 60
Pressure head (cm)
Distance from bottom (cm)
t = 5.76 s t = 29.16 s t = 70.56 s t = 129.96 s t = 207.36 s t = 302.76 s
(b)
Figure 7.6Moisture content (a) and pressure head (b) profiles. Example 2,ψtop= 100cm andqbottom= 0,∆t0= 0.01s.
method to converge. Moreover, for very small time steps the difference between the two methods is not too great, and thus, if the goal is to produce very accurate solutions the mixed–ψ method should be adequate. Thus, the main virtue of the new method lies in its ability to perform very well for large time steps. Moreover, the method should also be ideally suited to produce steady state solutions, where the initial guess is usually quite critical for the chances of convergence.
Finally, it should be emphasized that the new method can be coded by some quite insignif-icant modifications to existing mixed–ψcodes. Furthermore, the additional cost involved with the new method as compared to the conventional one is also insignificant since the
Department of Civil Engineering - Technical University of Denmark 69
Numerical solution procedures 7.5 A new mixed–variable iterative procedure
Mixed–ψ New method
∆t0(s) No. Step No. Iter Iter/Step No. Iter Iter/Step
10.0001 1733 6547 3.78 5495 293.17
10.0002 1225 4964 4.05 4076 293.33
10.0005 1775 † † 2837 293.66
10.0015 1548 † † 2195 294.01
10.0105 1174 † † 2959 295.52
10.1005 1755 † † 2483 298.78
11.0005 1718 † † 2254 14.1
10.0005 1786 † † 2174 29.0
Table 7.5Comparison of conventional mixed–ψ scheme and new iterative procedure. Ex-ample 3,ψinit= −10,000cm, ψtop = 100cm andqbottom = 0. †Convergence not achieved after 200 iterations.
bulk of the computations will be concerned with the computation of the pressure head increments, for which the two methods require the solution of the same set of equations.