Acta Polytechnica Acta Polytechnica 53(4):347–358, 2013 © Czech Technical University in Prague, 2013 available online at http://ctn.cvut.cz/ap/ ALGORITHMS FOR SOLVING DARCIAN FLOW IN STRUCTURED POROUS MEDIA Michal Kuráža,∗, Petr Mayerb a Czech University of Life Sciences Prague, Faculty of Environmental Sciences, Department of Water Resources and Environmental Modeling b Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mathematics ∗ corresponding author: kuraz@fzp.czu.cz Abstract. This paper presents several algorithms that were implemented in DRUtES [1], a new open source project. DRUtES is a finite element solver for coupled nonlinear parabolic problems, namely the Richards equation with the dual porosity approach (modeling the flow of liquids in a porous medium). Mass balance consistency is crucial in any hydrological balance and contaminant transportation evaluations. An incorrect approximation of the mass term greatly depreciates the results that are obtained. An algorithm for automatic time step selection is presented, as the proper time step length is crucial for achieving accuracy of the Euler time integration method. Various problems arise with poor conditioning of the Richards equation: the computational domain is clustered into subregions separated by a wetting front, and the nonlinear constitutive functions cover a high range of values, while a very simple diagonal preconditioning method greatly improves the matrix properties. The results are presented here, together with an analysis. Keywords: Richards equation, time step adaptivity, distinct unsaturated hydraulic properties, problem conditioning, diagonal preconditioner. 1. Introduction It is important to be able to predict fluid movement in an unsaturated/saturated zone in many fields, rang- ing from agriculture, via hydrology, to technical ap- plications of dangerous waste disposal in deep rock formations. The mathematical model of unsaturated flow was originally published by L. A. Richards [2]. The Richards equation problem has undergone various investigations and numerical treatments. Its finite element solution was originally published by Neuman in 1970 for several engineering applications, e.g. dam seepage modeling, see [3]. The existence and the uniqueness of its solution was discovered 10 years later, by Alt & Luckhaus [4]. A fundamental work analyzing a mass conservation numerical method for the Richards equation was published in 1990 by Celia et al. [5].An analysis of the Richards equation prob- lem has been the subject of several works by Kačur, e.g. [6]. Various methods for numerical treatment of the problem have recently been published, see e.g. Kees et al. [7], Šolín et al. [8] or Kuráž et al. [9, 10]. An efficient algorithm for numerical solution of the Richards equation should handle temporal adaptivity. This topic has already been discussed by a variety of authors. The first work maintaining the truncation error of the mixed form of the Richards equation was published by Kavetski et al. [14]. Tocci et al. [15] proved that the h-based Richards’ equation with proper temporal adaptivity can provide a mass conser- vative numerical solution, and later Miller et al. [16] published an algorithm for temporal adaptation of the step and degree of approximation. This paper is organized as follows. The first part discusses problems related to inaccurate time integra- tion of the Richards equation. A new and efficient automated time step selection method is presented here, and its efficiency is evaluated on a simple bench- mark problem. The second part of the paper discusses poor conditioning of the Richards equation problem. A simple and extremely efficient matrix preconditioner based on diagonal scaling greatly improves the con- ditioning, and is presented here. Diagonal precon- ditioning, or more precisely Jacobi preconditioning, is one of the first preconditioning methods that ever appeared in computational history, and its theoretical background was originally published more a hundred years ago. We have not succeeded in our search for the original reference. This method is a very classical approach, and it is even a subject of early courses in applied mathematics. However, it has been found that this standard method is highly efficient for solv- ing the systems of linear equations that arise from numerical approximation of the Richards equation. It has been observed that this preconditioning method even enables very fast convergence of the conjugate 347 http://ctn.cvut.cz/ap/ Michal Kuráž, Petr Mayer Acta Polytechnica 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 -200 -150 -100 -50 0 50 θ [ -] capillary pressure head (h) [cm] claystone gravel concrete 1e-10 1e-08 1e-06 0.0001 0.01 1 100 10000 -200 -150 -100 -50 0 50 K (h ) [c m .d -1 ] capillary pressure head (h) [cm] claystone gravel concrete Figure 1. Plot of the constitutive functions for the Richards equation (1). Left – the van Genuchten function – retention curve – relation between water content θ and the capillary pressure head h (mass function). Right – the Mualem function – unsaturated hydraulic conductivity) gradient method for normal equations, despite pow- ering the condition number. Several more or less recent remarks on using the diagonal preconditioner for CFD and porous media transport problems can be found in [20] and [21],but we have not succeeded in a search for a publication demonstrating the high efficiency of this algorithm as observed. The efficiency of this method will be presented on the basis of one theoretical application and one practical application. However, during the last two years we have not dis- covered any Richards equation problem where this preconditioning method performed in a significantly different way than as presented here. 2. Mathematical Model The problem of Darcian flow in a structured porous medium is usually expressed by the dual porosity con- ceptual approach. The governing equations for vari- ably saturated Darcian flow and contaminant trans- port in the dual flow regime were originally published by Gerke, et al. [11]. The following section introduces the strong and weak formulation of dual porosity variably saturated flow. A very brief description is given of the constitu- tive relations and coefficients involved in this problem, and an interested reader can find more information in the references provided below. 2.1. Strong and Weak Formulations The mathematical problem to be solved is the Richards equation with the dual porosity conceptual approach. Let Ω be a bounded domain in 0 and measn−1(ΓNm) > 0. Let the same hold for the pair ΓDf , ΓNf . Let n denotes the outer unit normal to ∂Ω. For positive T let QT = Ω× [0, T ). The Richards equation in the dual regime was pre- sumed by Gerke & Genuchten [11] as Cm(hm)∂hm ∂t = ∇. (Km(hm)∇hm) + ∂Km(hm) ∂z + Γw 1− ωf , (x, t) ∈ QT , Cf (hf )∂hf ∂t = ∇. (Kf (hf )∇hf ) + ∂Kf (hf ) ∂z − Γw ωf , (x, t) ∈ QT , Γw = αw(hf − hm), (1) where the subscripts f and m denote the subsys- tem of fractures (macropores) and matrix blocks (mi- cropores), h is the capillary pressure head function [L], K(h) is the unsaturated hydraulic conductivity [LT−1], C(h) is the water retention capacity [L−1], usually defined as C(h) = θ(h)′ + θ(h) θS Ss, where θ(h) is the water content function [–], Ss is the specific aquifer storage [L−1], θS is the saturated water con- tent [–], ωf denotes the volume ratio of the fracture pore system related to the total pore system, and αw is the first order mass transfer coefficient [L−1 T−1] presumed as αw = β α2 DP Kaγw, (2) where β is the dimensionless geometry coefficient, αDP is the characteristic half width [L] of the matrix block, Ka is the effective hydraulic conductivity [LT−1] of the matrix at or near the fracture/matrix interface, and γw is the dimensionless scaling factor. The consti- tutive relation for function θ(h) was supplied by the van Genuchten law [12], and for the K(h) function it was supplied by the Mualem law [13], see Figure 1 for a plot of these functions. The initial conditions are stated as hm(x, 0) = h0 m(x) ∀x ∈ Ω, hf (x, 0) = h0 f (x) ∀x ∈ Ω, (3) 348 vol. 53 no. 4/2013 Algorithms for Solving Darcian Flow in Structured Porous Media and the boundary conditions hm(x, t) = hDm(x, t), (x, t) ∈ ΓDm × [0, T ), ∂hm ∂n = hNm(x, t), (x, t) ∈ ΓNm × [0, T ), hf (x, t) = hDf (x, t), (x, t) ∈ ΓDf × [0, T ), ∂hf ∂n = hNf (x, t), (x, t) ∈ ΓNf × [0, T ). (4) We assume that all coefficients are sufficiently smooth. Let E1(Ω) := { v ∈ C∞(Ω)n; supp v ∩ ΓDm = ∅ } , E2(Ω) := { v ∈ C∞(Ω)n; supp v ∩ ΓDf = ∅ } . (5) Let the linear space V1 be a closure of E1(Ω), and let V2 be a closure of E2(Ω), in the norm of W 1,2(Ω). Functions hm, hf ∈ C([0, T ); W 1,2(Ω)) ∩ C([0, T ); L∞(Ω))∩ C1([0, T ); L2(Ω)), such that hm(x, 0) = h0 m, hf (x, 0) = h0 f , for x ∈ Ω and hm(x, t) = hDm for x ∈ ΓDm × [0, T ) and similarly for hf , constitute a variational solution of system (1) together with initial condition and boundary conditions (3)–(4) iff∫ Ω Cm(hm)∂hm ∂t ϕ1dΩ + ∫ Ω (Km(hm)∇hm)∇ϕ1dΩ − ∫ Ω ∂Km(hm) ∂z ϕ1dΩ − ∫ Ω Γw 1− ωf ϕ1dΩ = ∫ ΓN m hNmϕ1dΓNm (6) and∫ Ω Cf (hf )∂hf ∂t ϕ2dΩ + ∫ Ω (Kf (hf )∇hf )∇ϕ2dΩ − ∫ Ω ∂Kf (hf ) ∂z ϕ2dΩ + ∫ Ω Γw ωf ϕ2dΩ = ∫ ΓN f hNf ϕ2dΓNf (7) hold for every ϕ1 ∈ V1, ϕ2 ∈ V2, and t ∈ [0, T ). The time derivative in equations (6)–(7) was han- dled with the backward time difference. The contin- uous unknown hn+1 m,f will be replaced by a piecewise- polynomial function (in our implementation the degree of the polynomial is p = 1) hn+1 m = Dm∑ j=1 am,jϕ1,j , hn+1 f = Df∑ j=1 af,jϕ2,j (8) and the test functions ϕ1 and ϕ2 with finite element test functions ϕi1 and ϕi2 to obtain∫ Ω Cm(hn+1 m )h n+1 m − hnm ∆t ϕi1dΩ + ∫ Ω ( Km(hn+1 m )∇hn+1 m ) ∇ϕi1dΩ − ∫ Ω K ′m(hn+1 m )∂h n+1 m ∂z ϕi1dΩ − ∫ Ω Γw 1− ωf ϕi1dΩ = ∫ ΓN m hn+1,N m ϕi1dΓNm (9) for all i = 1; 2; . . . ;Dm, Dm is the dimension of the finite element subspace of V1, and∫ Ω Cf (hn+1 f ) hn+1 f − hnf ∆t ϕi2dΩ + ∫ Ω ( Kf (hn+1 f )∇hn+1 f ) ∇ϕi2dΩ − ∫ Ω K ′f (hn+1 f ) ∂hn+1 f ∂z ϕi2dΩ − ∫ Ω Γw ωf ϕi2dΩ = ∫ ΓN f hn+1,N f ϕi2dΓNf (10) for all i = 1; 2; . . . ;Df , Df is the dimension of the finite element subspace of V2. The superscript n enumerates the time steps, and ∆t is the length of the time step. The nonlinear operator was treated by the standard Picard method. This method performs a simple lin- earization of nonlinear product terms by moving all of them but one to the previous iteration level. For equations (9)–(10) the following set is obtained∫ Ω Cm(hn+1 m,k ) hn+1 m,k+1 − hnm ∆t ϕi1dΩ + ∫ Ω ( Km(hn+1 m,k )∇hn+1 m,k+1 ) ∇ϕi1dΩ − ∫ Ω K ′m(hn+1 m,k ) ∂hn+1 m,k+1 ∂z ϕi1dΩ − ∫ Ω Γw 1− ωf ϕi1dΩ = ∫ ΓN m hn+1,N m ϕi1dΓNm (11) for all i = 1; 2; . . . ;Dm, and∫ Ω Cf (hn+1 f,k ) hn+1 f,k+1 − hnf ∆t ϕi2dΩ + ∫ Ω ( Kf (hn+1 f,k )∇hn+1 f,k+1 ) ∇ϕi2dΩ − ∫ Ω K ′f (hn+1 f,k ) ∂hn+1 f,k+1 ∂z ϕi2dΩ − ∫ Ω Γw ωf ϕi2dΩ = ∫ ΓN f hn+1,N f ϕi2dΓNf (12) for all i = 1; 2; . . . ;Df . Here the subscript k denotes the iteration level. The stopping criterion was ob- tained by looking at the distance of two consecutive 349 Michal Kuráž, Petr Mayer Acta Polytechnica Figure 2. The critical point of solution h was as- sumed from nodal values of the discretization mesh. Then the only acceptable solution in the next time level is a solution that does not exceed the ranges defined by some critical error of the first order Taylor series approximation of solution h at the previous time level. Picard approximations of the solutions in the both domains – ‖hn+1 m,k+1 − h n+1 m,k ‖2 + ‖hn+1 f,k+1 − h n+1 f,k ‖2 < ε, (13) where ε is the desired accuracy of the method. Equa- tions (11)–(12) constitute the system of linear alge- braic equations to be solved in each iteration. The solution of this problem is scattered into hm, and hf . However the scattered outputs hm, hf are not as interesting as the behavior of the capillary pressure head h, which is obtained by averaging the matrix domain m and the fast domain f : h = (1− ωf )hm + ωfhf . (14) One could also hardly find any real problem, where such scattering could be required in the boundary and initial condition setup. Thus the initial condi- tion and the Dirichlet boundary condition are always considered identical for the both domains, and the Neumann boundary condition (volume flux) is always proportionally scattered into the both domains. 3. Time step adaptivity It is well known that incorrect discretization of the mass term (time derivative term) yields mass balance errors. This topic, particularly for the Richards equa- tion, was already documented by Celia, et. al. [5] in the 1990s. The primary solution of the Richards equa- tion in h-based form is the capillary pressure head h. The Richards equation law originates from the law of mass balance, and thus the time derivative term can be referred to as a mass term. The original form of the Richards equation published by Richards [2] con- siders the mass term to be of the form ∂θ(h) ∂t , where θ(h) is the water content function (mass function). The so-called h-based form of the Richards equation considers the mass term in the form C(h)∂h∂t , where C(h) = θ(h)′ (if zero specific storage is assumed). Due to the non-linearities involved in the constitutive mass function, see Figure 1, an efficient and accurate time integration algorithm based on the Euler method should handle the selection of an adaptive time step. Time step adaptivity has been the subject of several works. Most of them focus on limiting the truncation error in the mass term. The original work presenting an implementation of the method of lines for limiting the truncation error of the mixed form of the Richards’ equation was published by Kavetski et al. [14]. Based on Celia’s work [5] the h-based form of the Richards’ equation, if numerically approximated by the implicit Euler method, is inappropriate due to the huge inac- curacies in the mass balance. However, it was later presented by Tocci et al. [15], that limiting the trunca- tion error in function h by adaptive time step selection preserves the mass balance and the h-based form of the Richards equation becomes a suitable governing equation for subsurface transport processes. In the works of both Kavetski et al. and Tocci et al., the trun- cation error was limited by applying the method of lines with a fixed order of approximation. Even more sophisticated work was later published by Miller et al. [16], where temporal adaption was accomplished by variable order and variable step size approximations. The proposed time integration method does not limit the truncation error of the time derivative term, as mentioned in the previous paragraph. It limits the error of the h-based approximation in a discrete analogue of the mass term (∂θ(h) ∂t ≈ ∆θ(h) ∆t ≈ C(h)∆h ∆t ). This error can be controlled by defining the ranges for the shift of the solution h in two sequential time steps. The method is depicted in Figure 2. The range of the accepted shift of the solution h in two sequential time steps is a function of h. It will be referred to as hhigh and hlow, and thus the following implicit functions are formed: |θ(h)− C(h)(h− hlow)− θ(hlow)| − ε = 0, |θ(h) + C(h)(hhigh − h)− θ(hhigh)| − ε = 0, (15) where ε is the error constant – an error of the first order Taylor series approximation of the constitutive mass function, see Figure 2. This implicit equation was solved by the bisection method. Because it is a resource-consuming computer operation, the equation is tabelized by the initialization procedures. Start- ing with some initial time step guess, the solution is checked always at the first and last Picard iteration level (the last iteration level of the Picard method is the level, where the distance between two consecutive Picard’s approximations is lower than a certain value (‖hn+1 k+1 − h n+1 k ‖2 < ε – see section 2.1). The range 〈hlow, hhigh〉 is evaluated for all mesh nodal values by solving equation (15), where value h is supplied by the previous time step function. It is apparent that evaluating this implicit function directly at each time 350 vol. 53 no. 4/2013 Algorithms for Solving Darcian Flow in Structured Porous Media Figure 3. Scheme of the one-dimensional vertical infiltration experiment. step for each nodal value would markedly influence the performance of the code, and therefore an efficient and simple algorithm that searches the values in a table was constructed. The solutions of (15) (hlow and hhigh) are precal- culated for several discrete values of the h function within some user defined ranges and density. Then the solution of (15) for any arbitrary h value is obtained by the following algorithm: • the h value is converted into an integer value that points the position of the nearest smallest discrete value of the tabelized function (15), • then the values of solutions hlow and hhigh are ob- tained by a linear interpolation between the nearest discrete values of the tabelized function (15). In order to speed up the code performance this technique is also offered to evaluate values for all ma- terial functions with real exponents. The user can switch between this approximate and direct mate- rial function evaluation method. This method was tested on a standard laptop (Intel CPU 1.6GHz) for a mesh of 106 elements with the nine points integration method. If constitutive functions were evaluated di- rectly, the finite element method matrix was built in approximately 20 min.; if constitutive functions were tabelized as described, the solution matrix was created in approximately one minute. The array initialization took not more than a few seconds. This experiment in- volved only the van Genuchten’s constitutive functions without evaluating the function (15) [24]. If some nodal value in a newly evaluated solution falls beyond the range 〈hlow, hhigh〉, the new time step solution is rejected, the time step is decreased, and the new solution is tested again – until the criterion is fulfilled. It can be assumed that the constant ε controls the error caused by discretization of the mass term in the h-based form of the Richards equation.The method will be further referred as to the Mass Curve Zone Approach. It should be emphasized that equation (15) has a unique solution only for monotone convex retention curves, e.g. the Gardner’s model [17]. For the van Genuchten function this algorithm uses only the solu- tion hlow, resp. hhigh, which is closer to the particular evaluated h value. 3.1. Case study A simulation of a vertical one-dimensional infiltra- tion experiment was conducted, see Figure 3. For this numerical example, the dual porosity approach was found to be useless for demonstration purposes. The problem was therefore evaluated on the classical Richards equation only. The selected flow domain Ω was a porous material block 5 m in length and of infinite width. Thus only vertical one-dimensional flow was taking place. Porous material properties (parameters of the van Genuchten and Mualem consti- tutive functions) were assumed for highly permeable sandstone (α = 0.075 cm−1, n = 1.89, m = 0.47, Ks = 4.42 cmh−1). The domain was discretized by a uniform mesh of grid size 5 cm. The non-linear Richards equation operator was lin- earized by the Picard method, the stopping criterion was defined as ‖hk+1 − hk‖2 ≤ 10−3, (16) where subscript k is an iteration level. In order to evaluate the method presented here, the simulation was performed three times, each time with a different temporal adaptivity mechanism. All selected adaptivity methods were purely empirical. (1.) Number of iterations of the Picard method (max. 10 iterations with a defined stopping criterion) — h-based form of the Richards equation. (2.) Number of iterations of the Picard method (max. 10 iterations with the defined stopping criterion) — mixed form of the Richards equation (the discrete analogue of the mass term was considered as ∆θ(h) ∆t ). (3.) Mass Curve Zone Approach (the maximal error of a first-order Taylor series, see Figure 2, was 5 ·10−6). The initial condition was defined as h(x, 0) = −8700.0 cm ∀x ∈ Ω. (17) The top boundary condition was defined as h(x, t) = 5.0 cm ∀(x, t) ∈ Γtop × (0, T ]. (18) And the bottom boundary was the homogeneous Neumann boundary condition (in technical practice this type of the bottom boundary setup is referred to as free drainage). 351 Michal Kuráž, Petr Mayer Acta Polytechnica Method Mass error Number [%] of iterations 1st 5.842 130 2nd 0.201 44382 3rd 0.193 1746 Table 1. Summary of code execution after the nu- merical infiltration experiment. 3.2. Results and Conclusions As stated in [5] the successful iteration criterion for adaptive time-step selection of the h-based Richards equation is not a sufficient condition to ensure ac- curacy. The mixed form of the Richards equation provides very accurate results, but the computational effort is excessive in the case evaluated here. The Mass Curve Zone Approach provides results of similar quality, but the total number of iterations (number of calls of the linear equation solver) is only approx- imately 4% of the amount required for solving the mixed-form Richards equation. This method offers an efficient condition for a nu- merical solution of the h-based equation based on the implicit backward Euler scheme, exactly reflecting the problem of the h-based approximation. Compared to the Tocci approach [15], our algorithm does not limit the truncation error of the time derivative term, it only limits the error of the discrete analogue of the mass term if the h-based form of the Richards’ equa- tion is considered. However, our algorithm behaves in a different way. It should be evaluated whether there is any relation between these two errors, or at least which criterion is more appropriate for achieving mass consistency. For particular results see Table 1. 4. Preconditioner Finite element approximation of the Richards equa- tion problem combined with the fully implicit Roth method typically forms very badly-conditioned ma- trices, especially if coarse-grain material is involved. The plot of the unsaturated hydraulic conductivity, Figure 1 right, covers a typical capillary pressure head range for a certain Richards equation problem, and thus it is not uncommon to solve problems where the non-linear model parameters can vary between eight or even ten orders of magnitude. The poor condition- ing of the Richards equation is caused by the huge range of values of the constitutive functions. Another property of the Richards equation also applies. Due to the small time steps that are typically required for an accurate solution of this problem, and due to lump- ing of the capacity matrices in order to preserve the discrete maximum principle, the matrices of its numer- ical approximation are diagonally dominant. Then it was observed that badly conditioned matrices of this type are formed out of clusters with a lower condition number. For simple infiltration problems, the clusters are formed at (1.) fracture and matrix flow domains, (2.) regions separated by the wetting front. And thus it can be summarized that, • due to the small time steps, matrices are diagonally dominant, • due to the small time steps and the continuity of the solution (it is in fact just a non-linear diffusion with convection), the matrices are "similar" in two sequential time steps, • the domain can be clustered into subdomains with similar values of the constitutive functions, and thus the matrix can be clustered into "submatrices" with a lower condition number (as the problem is nonlinear and often unsteady, the position of these subdomains is not constant in time). Based on these assumptions, the perfect precondi- tioner should adjust the material properties in the computational domain into such a state, that the ma- terial parameters (values of the constitutive functions K(h, x) and C(h, x))are transformed into certain ho- mogeneous values. Behavior of this kind is actually achieved by the standard Jacobi preconditioner (diag- onal scaling) assumed from the right hand side. The Jacobi preconditioner is one of the simplest forms of preconditioning, in which the preconditioner matrix is the diagonal of the original matrix A. The method can be described by the following scheme  a11 a12 . . . a1n a21 a22 . . . a2n ... ... . . . ... an1 an2 . . . ann  →  1.0 a12 a22 . . . a1n ann a21 a11 1.0 . . . a2n ann ... ... . . . ... an1 a11 an2 a22 . . . 1.0  , (19) and thus this method can be referred to as a right hand side diagonal preconditioner Aprcd = A×D, (20) where D =  1 a11 0 . . . 0 0 1 a22 . . . 0 ... ... . . . ... 0 0 . . . 1 ann  . (21) This method suffers from one notable disadvantage – it creates non-symmetric matrices even out of sym- metric matrices. However, transport processes are in 352 vol. 53 no. 4/2013 Algorithms for Solving Darcian Flow in Structured Porous Media general a convection-diffusion process, and thus they are non-symmetric. In the context of the conjugate gradient method for normal equations, this method of preconditioning can be formulated as DTATAD, which is actually a standard approach. 4.1. On Selecting an Algorithm for Solving the Systems of Linear Equations Transport processes in porous media are in general of a convection-diffusion-reaction form. In the case of the standard Richards equation the convection is caused by gravitational forces, and then the problem can be formulated symmetrically with respect to the total hydraulic head H instead of the pressure head h. This approach would be efficient for the case studies evaluated here, but the transport processes in porous media are more accurately described by a coupled process of moisture, heat and solute transport, as described by Noborio et al. [18]. Thus the search for an efficient solver should not be limited to symmetric problems only. Moreover, preserving the formulation of the Richards equation in the pressure head form h enables density dependent transport problems to be solved, e.g. sea water intrusion. This form was therefore preferred in implementing the DRUtES code. Secondly, the diagonal scaling method that is ap- plied forms non-symmetric matrices out of symmetric matrices. This property is obviously inconvenient, but based on the following assumption the cost for losing the symmetricity, either due to a less appropri- ate formulation of the Richards equation or due to the application of right hand side diagonal scaling, is insignificant. Let A be the matrix obtained by a numerical ap- proximation of the Richards equation problem into the system of linear equations Ax = b, as described in sec- tion 2.1. For each iteration of the Picard method, ma- trix A and vector b is assembled from equations (11)– (12) and from the boundary condition setup (4). Due to nonlinearities in the Richards equation, and due to preserving the mass balance, the problem is tem- porally discretized into small time steps. Together with lumping the capacity matrix as recommended by Rank et al. [19], in order to preserve the discrete max- imum principle, we could start with the hypothesis of a strong diagonal dominance of matrix A. Then matrix A can be decomposed into the sum of its diagonal and nondiagonal part A = D + B, (22) where D = diag(A), and obviously B = A−D. The right hand side diagonal scaling rearranges the original problem Ax = b into AD−1y = b, (23) where x = D−1y. Then the matrix AD−1 can be considered as AD−1 = I + BD−1, (24) where I is an identity matrix. If we apply the conju- gate gradient algorithm for normal equations to the matrix I + BD−1 from (24), then the following is obtained G = ( I + BD−1)T (I + BD−1) = = I + (BD−1)T + BD−1 +(BD−1)TBD−1, (25) where G is the matrix that is processed in the conju- gate gradient algorithm. Based on the hypothesis of the strong diagonal dom- inance of matrix A, the Frobenius norm of ‖BD−1‖ is small ‖BD−1‖ � 1.0, (26) and thus it could be concluded that Frobenius norm of the matrix G from (25) must be close to 1.0. This conclusion will be deeply explained. Let us consider matrix F as F = G− I, (27) and thus F = (BD−1)T + BD−1 + (BD−1)TBD−1. (28) Then the Frobenius norm of matrix F is referred to Fne = ‖F‖. (29) We have already assumed the Frobenius norm of ‖BD−1‖ to be small, see (26), due to the triangle in- equality and due the submultiplicativity of the Frobe- nius norm we can conclude that Fne ≤ 3‖BD−1‖. (30) In order to estimate the efficiency of the conjugate gradient solver applied to the matrix G from (25), we should estimate the condition numbers of matrix F from (28). It is well known that the absolute values of all eigenvalues are lower than its matrix norm (this applies for all commonly used norms including the Frobenius norm). And thus the following applies σ (F) ⊂ (−Fne,+Fne). (31) And thus matrix G from (25), which is related to matrix F by definition (27), has eigenvalues within the range σ (G) ⊂ (1− Fne, 1 + Fne). (32) The condition number is evaluated as κ = 1 + λF max 1− λF min , (33) where κ is the condition number of matrix G, λF min is the minimal eigen value of matrix F, and λF max is the maximal eigen value of matrix F, and thus 1 + λF max 1− λF min < 1 + Fne 1− Fne . (34) 353 Michal Kuráž, Petr Mayer Acta Polytechnica Figure 4. Scheme of the two-dimensional theoretical experiment. In order to fulfill the proposed prerequisite, the norm ‖BD−1‖ must be small, but as already men- tioned the value of this norm for our Richards equation problems was typically observed within the range ‖BD−1‖ ∈ (1 · 10−12, 1 · 10−8). (35) The conjugate gradient method for normal equa- tions is not a typical choice due to squaring the condi- tion number. The typical solution is based on GMRES or BiCGSTAB. GMRES is an optimal method on the Krylov subspace generated by the initial residual r0 and matrix A. But its great disadvantage is that long recurrences are required, so it is almost impossible to use the original form of GMRES, as the long re- currences may cause serious problems with computer memory allocation. GMRES is typically used in a vari- ety of modifications, as restarts or shorten recurrences. As presented by Greenbaum et al. [22] the convergence rate does not depend on the eigenvalues, and thus achieving an accurate solution in restarted GMRES is uncertain (at least on theoretical problems). The BiCGSTAB method uses short recurrences, but its great disadvantage is the possibility of a so-called break-down. The simple conjugate gradient method for normal equations does not suffer from any these disadvantages. The only disadvantage is squaring the condition num- ber, but on the basis of the given hypothesis for this particular class of problems worsening the condition number is negligible, and thus this only disadvantage is absent here. The conjugate gradient method for normal equa- tions is implementary simple, and as explained for diagonally scaled diagonally dominant problems, it is also an efficient algorithm. This method was im- plemented into DRUtES code, and bas been used successfully for a variety of theoretical and practical problems. 4.2. Theoretical Case Study This example is a purely theoretical case study. The matrix and fracture domain parameters are consid- 1 10 100 1000 10000 0 10 20 30 40 50 60 n u m b e r o f it e ra ti o n s time [hrs] preconditioning no preconditioning Figure 5. Conjugate gradient iterations in time. ered identical here, and thus the problem is described by the classical Richards equation model only. The material parameters are depicted in Table 3, 3rd line “gravel bar.”. The boundary condition setup is stated as: h(x, t) = −10.0 cm ∀(x, t) ∈ ΓD × (0, T ], ∂h ∂n (x, t) = 0 ∀(x, t) ∈ ΓN × (0, T ]. (36) The initial condition setup is stated as: h(x) = −200.0 cm ∀x ∈ Ω. (37) The domain scheme is depicted in Figure 4. The spatial discretization was a uniform triangular mesh with 10201 nodes and 20000 elements. 4.2.1. Conditioning Improvements The diagonal scaling algorithm greatly improved the problem conditioning. The overall number of Picard iterations — linear equation solver calls — was 49001. The matrix of dimension 10100 was solved by the con- jugate gradient method for a normal equation with an average of 33 iterations, and the maximal number of iterations required per CG solver call was 67. The stopping criterion was the residual size rk < 10−12. Table 2 depicts the results that were obtained. Fig- ure 5 depicts the plot of the conjugate gradient it- erations for preconditioned and non-preconditioned matrices. As expected, the CG method for a nor- mal equation without preconditioning converges so slowly that its application can be concluded as an inappropriateconsidered inappropriate. This is a well- known fact, but simple diagonal scaling enables fast convergence. 4.3. Practical Case study This example is motivated by an existing nuclear waste repository for low-level and intermediate-level waste. The name of the repository is Richard and it is situated near Litoměřice, Czech Republic. The repository is located within a deep rock formation 354 vol. 53 no. 4/2013 Algorithms for Solving Darcian Flow in Structured Porous Media Time No preconditioning Preconditioned λmax λmin κ λmax λmin κ 10−5 7.853 · 10−10 6.376 · 10−12 123.1612 0.99999999898 0.9999999919 1.00000000706 0.5 4.578 · 10−7 1.207 · 10−11 37926.33 0.9992155771 0.94875498673 1.053186113526 70.5 6.842 · 10−5 7.040 · 10−10 97199.057 2.135820803 0.28957144730 7.375799041142 154.5 0.000356 1.773 · 10−8 20109.385 2.9932628563 0.08109150218 36.91216435282 276.5 0.000103 3.553 · 10−9 29069.776 2.4120151933 0.20777359355 11.60886305218 339.5 0.000257 1.948 · 10−8 13234.228 2.869385334 0.1028734952 27.892367499 502.5 0.000410 3.371 · 10−8 12176.805 3.0701098732 0.06920022271 44.36560682 669.0 5.0563 · 10−5 1.818 · 10−9 27805.705 2.01150608 0.34435293431 5.841408287719 850.0 0.0005709 1.090 · 10−7 5236.54 3.196279963 0.0527000983 60.65036051186 950.0 0.0007849 2.026 · 10−7 3873.438 3.3012379801 0.04053514763 81.44137057972 Table 2. Conditioning improvements before and after right hand side diagonal scaling. If λmax < 1.0, then most likely more iterations are required to improve the accuracy of the eigenvalue estimation. Figure 6. Domain scheme of the repository problem. in the Bídnice hill, 40m below the surface and 70m above the ground water table, and it is surrounded by tectonically fractured claystone. Under normal circumstances, this is a very stable unsaturated en- vironment. This is therefore a typical problem for applying the Richards equation model. In order to protect its surroundings from radiotoxic pollution, the repository is shielded by a 30 cm thick gravel layer that acts as a capillary barrier. The capillary barrier effect is based on the unsaturated hydraulic properties of gravel. The unsaturated hydraulic conductivity of gravel reaches very high values at and near full satu- ration, but with decreasing saturation it falls steeply to almost zero values. The problem was conceptualized by model (1), where the tectonic fractures are represented by the fast domain (f), and the rock matrix is represented by the matrix domain (m), the material properties are depicted in Table 3, and the domain scheme is depicted in Figure 6. The coupling parameters of the model (1) were identified in [9], and are stated as: α = 8.89911 cm−1 d−1, and ωf = 0.0292884. The hy- draulic parameters (parameters of the van Genuchten and the Mualem function) are depicted in Table 3, and Figure 1 is a plot of the constitutive functions. The repository interior – concrete, and the repository barrier – gravel, were obviously formed out of identical material in both the matrix domain and the fracture domain. The different material conditions scattered into the matrix and fracture domains appeared in the rock material. The initial condition was defined as h(x, 0) = −200.0 cm ∀x ∈ Ω, (38) and the boundary conditions were defined as • Γ1 – Dirichlet boundary condition of pressure head +10.0 cm hm,f (x, t) = 10.0 cm x ∈ Γ1 × [0, T ) (39) • Γ2,3,4 – homogeneous Neumann boundary condi- tions that satisfy the initial condition (free drainage boundary, because in general the Neumann condi- tion for the Richards equation states that qn = −n(K(h)∇h + K(h)(:,2)), where K(h)(:,2) is the second column of the hydraulic conductivity tensor in the case of <2) The domain was discretized by a non-uniform tri- angular mesh that consists of 12027 nodes and 23771 elements. 4.4. Numerical Properties The material properties of this problem cover several orders of magnitude, and thus finite element approxi- mation forms ill-posed matrices. In order to evaluate the efficiency of the preconditioner that is applied here, several early states were solved both with and without preconditioned matrices. When the stopping criterion definition from the previous paragraph was used, it was observed that for non-preconditioned ma- trices the number of iterations reached an average 355 Michal Kuráž, Petr Mayer Acta Polytechnica Material Ks [cmd−1] θr [–] θs [–] αvg [cm−1] nvg [–] mvg [–] Ss [cm−1] Claystone – “m” 0.0864 0.064 0.14 0.01 2.5 0.864 0.1 Claystone – “f” 15.9464 0.0 0.91087 0.00682 5.4829 0.43369 0.1 Gravel bar. 712.8 0.0 0.43 1.0 2.0 0.626 0.2 Concrete fill 1.68 0.089 0.43 0.01 1.23 0.187 0.03 Table 3. The material properties table. 0 50 100 150 200 250 300 350 400 0 25 50 75 100 p c g i te ra ti o n s time [days] right hand s. precond. left hand s. precond. Figure 7. Number of iterations required by the conju- gate gradient algorithm to converge for the repository problem. value of 17687. In the case of a preconditioned state, the average number of iterations was 18. The matrix dimension was 23898. Thus, for this particular case, this simple preconditioner decreased the total number of iterations required to solve the unsteady problem by a multiple of approximately 10−3. As was mentioned above, this proposed method suffers from a noticeable disadvantage – for symmetric problems the matrix is converted into a non-symmetric matrix. Thus the diagonal scaling was also assumed in a different way – from the left hand side Aprcd = D×A, (40) and then the scaling can be described by the following scheme  a11 a12 . . . a1n a21 a22 . . . a2n ... ... . . . ... an1 an2 . . . ann  →  1.0 a12 a11 . . . a1n a11 a21 a22 1.0 . . . a2n a22 ... ... . . . ... an1 ann an2 ann . . . 1.0  (41) It is well known, see eg. [23], that if A were a symmetric positive definite matrix, then the left and right hand side preconditioners would be equivalent. It should be further noted, that the left hand side Z [cm] −3400 −3200 −3000 −2800 −2600 −2400 −2200 h [cm] 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800X −2e+02 −1.5e+02 −95 −43 10 Figure 8. Intense initial infiltration – the Dirichlet boundary condition does not satisfy the initial condi- tion, such discontinuity causes intense processes in the initial states. preconditioner in the context of the conjugate gradient method for normal equations forms ATDTDA, and this form is not commonly used. This preconditioner was again evaluated in our repository study. It was observed that the left hand and right hand side preconditioners both behave simi- larly if the solution is "smooth", while for less smooth solutions the right hand side preconditioner offers bet- ter performance, see Figure 7. The "smoothness" of the solution was defined by the criterion sm = n−1∑ i=1 |xi+1 − xi|, (42) where n is matrix dimension, x is vector of solution (Ax = b), and xi is some particular value at position i from the solution vector x. Constant sm is only a certain indicator of the "smoothness" of the values in vector x. Then, as observed (see Figure 7), at time e.g. t = 2 days → sm ≈ 104 – the behavior of the two methods was comparable. However, for some later simulation time levels the right hand side preconditioner offered better performance, e.g. at time t = 80 days → sm = 106. If the right hand side preconditioner is used, together with the small time steps involved, the "smoothness" of the solution remains hidden in the diagonal matrix D. This is because the solution xprcd of matrix Aprcd is related to solution x of matrix A by x = D−1 × xprcd, (43) 356 vol. 53 no. 4/2013 Algorithms for Solving Darcian Flow in Structured Porous Media Figure 9. Wetting front is getting closer to the capil- lary barrier shield, the repository is being bypassed. and thus for the small time steps involved, together with an assumption that the solution of the Richards equation is continuous, the sm constant from equa- tion (42) of the xprcd vector reaches a low value. 4.4.1. Conditioning Improvements The unsteady problem of infiltration around the repos- itory can be divided into the following periods (1.) intense initial infiltration – the Dirichlet boundary condition does not satisfy the initial condition. Such discontinuity causes intense processes in the initial states, see Figure 8; (2.) the wetting front is getting closer to the capillary barrier shield, the repository is being bypassed, see Figure 9; (3.) the wetting front has reached the repository shield, steep gradients appear at the wetting front/barrier interface, see Figure 10. The right hand side preconditioner efficiency was evaluated for each period, for a selected time and it- eration level. Minimal and maximal eigenvalues were evaluated for the matrix AT ×A, as the problem is non-symmetric, and thus a conjugate gradient algo- rithm for normal equations was applied. Maximal eigenvalues were evaluated by the power method and minimal eigenvalues were evaluated by inverse itera- tion with Cholesky decomposition. The results are depicted in Table 4. 5. Conclusions This paper has aimed to present several new algo- rithms that were implemented in the DRUtES open- source project [1] – a finite element solver for nonlin- ear coupled parabolic problems, namely the Richards equation with the dual porosity approach, in order to improve the accuracy and efficiency of its numerical treatment. The Mass Curve Zone Approach, an efficient and ac- curate method for the Euler time integration method, has been presented and tested. This method is based on an automatic time step selection. The length of the time step is selected in such a way that the error Figure 10. Wetting front has reached the reposi- tory shield, steep gradients appears at the wetting front/barrier interface. of the discrete approximation of the time derivative term (mass term) is always kept below a certain value. A different problem concerned the conditioning of the Richards equation problem. It is well known that finite element approximation of the Richards equa- tion forms badly-conditioned matrices. A very simple and efficient preconditioning method based on diago- nal scaling has been presented. Then the method of preconditioning was tested on a real problem – infil- tration into a capillary barrier layer that was used as a protective shield around a waste repository. The method presented here efficiently decreased the condi- tion number of the evaluated problem by more than 13 orders of magnitude. The best performance was observed for the initial time levels, where the condi- tion number of the original matrix was κ = 1.43 · 1013, while after preconditioning κ = 1.00215, and thus the preconditioned matrix was even close to an identity problem. Generally speaking, all matrices that ap- peared in this simulation experiment had a condition number around 1013–1015, and the proposed precondi- tioner decreased the conditioning to values lower than 102–103. The efficiency of the proposed preconditioner is highly dependent on the selected time step. All simulated cases demonstrate great effectiveness of the discussed automated time step selection method for Euler time integration, and for the proposed matrix preconditioner. Generally speaking, both the design and the implementation of the preconditioner are trivial, and it is therefore a recommended method for conjugate gradient based finite element solvers of transport processes in a porous medium. Acknowledgements Financial support from the Czech Science Foundation (Research Project GACR 13-11977P) and from the Czech Ministry of Education, Youth and Sports (project COST CZ LD13068) is gratefully acknowledged. References [1] M. Kuráž, P. Mayer: DRUtES – finite element solver for nonlinear coupled convection-diffusion-reaction 357 Michal Kuráž, Petr Mayer Acta Polytechnica Time No preconditioning Preconditioned [days] λmax λmin κ λmax λmin κ 1.0 2.547 · 1013 1.7823 1.43 · 1013 1.0011 0.9989 1.00215 15.0 9.286 · 105 2.499 · 10−9 3.72 · 1014 12.1142 0.03508 345.335 120.0 9.112 · 105 1.391 · 10−9 6.55 · 1014 12.7672 0.02345 544.338 Table 4. Conditioning improvements before and after right hand side diagonal scaling. problems, http://www.freecode.com/projects/drutes, WWW [2] L.A. Richards: Capillary conduction of liquids through porous media, Physics I, 1, 11, 1931, 318-333 [3] S. P. Neuman, P. A. Witherspoon: Finite Element Method Analyzing Steady Seepage with a Free Surface, Water Resources Research, 6, 3, 1970, 889-897 [4] H. W. Alt, S. Luckhaus: Quasilinear Elliptic-Parabolic Differential Equations, Mathematische Zeitschrift, 183, 3, 1983, 311-341 [5] M. A. Celia, E. T. Bouloutas, R. L. Zarba: A General Mass-Conservative Numerical Solution for the Unsaturated Flow Equation, Water Resources Research, 26, 7, 1990, 1483-1496 [6] J. Kačur: Solution of Nonlinear and Degenerate Convection-Diffusion Problems, Nonlinear Analysis, 47, 1, 2001, 123-134 [7] C.E. Kees, M.W. Farthing, C.N. Dawson: Locally Conservative, Stabilized Finite Element Methods for Variably Saturated Flow, Comput. Methods Appl. Mech. Engrg., 197, 51, 4610–4625, 2008 [8] P. Šolín, M. Kuráž: Solving the nonstationary Richards equation with adaptive hp-FEM, Advances in Water Resources, 34, 9, 2011, 1062-1081 [9] M. Kuráž, P. Mayer, M Lepš, D. Trpkošová: An adaptive time discretization of the classical and the dual porosity model of Richards’ equation, Journal of Computational and Applied Mathematics, 233, 12, 2010, 3167-3177 [10] M. Kuráž, P. Mayer, V. Havlíček, P. Pech, J. Pavlásek: Dual permeability variably saturated flow and contaminant transport modeling of a nuclear waste repository with capillary barrier protection, Applied Mathematics and Computation, 219, 13, 2013, 7127-7138 [11] H. H. Gerke, M. T. H. van Genuchten: A Dual-Porosity Model for Simulating the Preferential Movement of Water and Solutes in Structured Porous Media, Water Resources Research, 29, 2, 1993, 305-319 [12] M. T. H. van Genuchten: A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Journal of Soil Science, 44, 5, 1980, 892-898 [13] Y. Mualem: A New Model for Predicting the Hydraulic Conductivity of Unsaturated Porous Media, Water Resources Research, 12, 3, 1976, 513-522 [14] D. Kavetski, P. Binning, S.W. Sloan: Adaptive time stepping and error control in a mass conservative numerical solution of the mixedform of Richards’ equation, Advances in Water Resources. 24, 2001, 595-605 [15] M.D. Tocci, C. T. Kelley, C. T. Miller: Accurate and Economical Solution of the Pressure Head Form of the Richards Equation by the Method of Lines, Advances in Water Resources, 20, 1997, 1-14 [16] C. T. Miller, C. Abhishek, M. W. Farthing: A spatially and temporally adaptive solution of Richards’ equation, Advances in Water Resources, 29, 2006, 525-545 [17] W. R. Gardner: Some Steady State Solutions of the Unsaturated Moisture Flow Equation with Application to Evaporation from a Water Table, Soil Science, 85, 1958, 228-232 [18] K. Noborio, K.J. McInnes, J.L. Heilman: Two-dimensional model for water, heat, and solute transport in furrow-irrigated soil: I. Theory, soil Science Society of America Journal, 60, 4, 1996, 1001–1009 [19] E. Rank, C. Katz, H. Werner: On the Importance of the Discrete Maximum Principle in Transient Analysis Using Finite Element Methods, International Journal for Numerical Methods in Engineering, 19, 1983, 1771–1782 [20] G. Pini, G. Gambolati: Is Simple Diagonal Scaling the Best Preconditioner for Conjugate Gradients on Supercomputers?, Advances in Water Resources, 13, 1990, 147-153 [21] X. Xu, N. Qin, B.E. Richards: Alpha-GMRES: A New Parallelizable Iterative Solver for Large Sparse non-Symmetric Linear Systems Arising from CFD, International Journal for Numerical Methods in Fluids, 15, 1992, 613-623 [22] A. Greenbaum, V. Pták, Z. Strakoš: Any nonincreasing convergence curve is possible for GMRES, SIAM Journal on Matrix Analysis and Applications, 17, 1996, 465-469 [23] M. Benzi: Preconditioning techniques for large linear systems: A survey, Journal of Computational Physics, 182, 2002, 418-477 [24] M. Kuráž: Numerical solution of the flow and transport equations in a porous medium with the dual permeability conceptual approach, Verbum publishing, Zlín, 2011 358 http://www.freecode.com/projects/drutes Acta Polytechnica 53(4):347–358, 2013 1 Introduction 2 Mathematical Model 2.1 Strong and Weak Formulations 3 Time step adaptivity 3.1 Case study 3.2 Results and Conclusions 4 Preconditioner 4.1 On Selecting an Algorithm for Solving the Systems of Linear Equations 4.2 Theoretical Case Study 4.2.1 Conditioning Improvements 4.3 Practical Case study 4.4 Numerical Properties 4.4.1 Conditioning Improvements 5 Conclusions Acknowledgements References