- Email: [email protected]
Solving the nonlinear Richards equation model with adaptive domain decomposition
Journal of Computational and Applied Mathematics 270 (2014) 2–11
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Solving the nonlinear Richards equation model with adaptive domain decomposition Michal Kuraz a,∗ , Petr Mayer b , Pavel Pech a a
Czech University of Life Sciences Prague, Faculty of Environmental Sciences, Department of Water Resources and Environmental Modeling, Czech Republic b
Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mathematics, Czech Republic
article
info
Article history: Received 10 October 2013 Received in revised form 16 February 2014 Keywords: Richards equation Highly heterogeneous material properties Additive Schwarz domain decomposition Nonlinear operator treatment
abstract Modeling the transport processes in a vadose zone plays an important role in predicting the reactions of soil biotopes to anthropogenic activity, e.g. modeling contaminant transport, the effect of the soil water regime on changes in soil structure and composition, etc. Water flow is governed by the Richards equation, while the constitutive laws are typically supplied by the van Genuchten model, which can be understood as a pore size distribution function. Certain materials with dominantly uniform pore sizes (e.g. coarse-grained materials) can exhibit ranges of constitutive function values within several orders of magnitude, possibly beyond the length of real numbers that computers can handle. Thus a numerical approximation of the Richards equation often requires the solution of systems of equations that cannot be solved on computer arithmetics. An appropriate domain decomposition into subdomains that cover only a limited range of constitutive function values, and that will change adaptively, reflecting the time progress of the model, will enable an effective and reliable solution of this problem. Parts of this problem have already been treated in our work Kuraz et al. (2014). This paper focuses on improving the performance of a nonlinear solver by locating the areas with abrupt changes in the solution. © 2014 Elsevier B.V. All rights reserved.
1. Introduction The problem of predicting fluid movement in an unsaturated/saturated zone is important in many fields, ranging from agriculture, and hydrology to technical applications of dangerous waste disposal in deep rock formations. The mathematical model of unsaturated flow was originally published by L.A. Richards [1]. 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 [2]. The existence and the uniqueness of its solution were discovered 10 years later, by Alt and Luckhaus [3]. An analysis of the Richards equation problem has been the topic of several works by Kačur, e.g. [4]. A variety of methods for the numerical treatment of the Richards equation have been published in recent times, see e.g. [5–8]. A distinction is often made between those domain decomposition algorithms that use overlapping subregions and those which do not [9]. The alternating Schwarz method is the oldest, dating back to the 19th century [10], and it belongs to the first category (overlapping subregions). This method, originally designed as a tool for analysis, has been the subject of many
∗
Corresponding author. Tel.: +420 777018743. E-mail address: [email protected] (M. Kuraz).
http://dx.doi.org/10.1016/j.cam.2014.03.010 0377-0427/© 2014 Elsevier B.V. All rights reserved.
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
3
more or less recent publications and computer applications, see e.g. [11]. An additive variant of this algorithm was well treated in the 1980s e.g. by Dryja and Widlund [12]. Cai and Dryja [13] have studied the convergence of nonlinear strongly elliptic equations discretized by the finite element method when the Schwarz additive algorithms are taken into account, together with the classical inexact Newton method. Cai and Dryja provided a proof of the convergence of this algorithm. The convergence rate is independent of the finite element mesh parameter and the number of subdomains used in the domain decomposition. It was also latter proved by Dryja and Hackbusch [14], that the nonlinear operator converges locally with the same asymptotic speed as the corresponding linear iteration applied to the linearized problem. A numerical solution of the Richards equation typically encounters two major obstacles—the poor conditioning of the arising system of linear equations, and the slow convergence or even divergence of the nonlinear operator. The aim of this paper is to investigate a method that adaptively splits the computational domain, so that the nonlinear operator always iterates only on small matrices, as the convergence problem of the Richards equation is in most cases only a local problem around the wetting front. This paper extends and improves the results from the authors’ previous work [15], which studied the adaptive decomposer and the conditioning of the subdomain matrices. We can already present the results obtained from the algorithm that iteratively solves the nonlinearities on local matrices only. In our implementation, a combination of the Schwarz–Picard algorithm is considered rather than the Schwarz–Newton algorithm, as presented in [13]. We enable a technique that controls the iteration updates and the local residuum value, in order to localize the solution disturbances, so that the computational power can be focused on local problems locally, rather than on solving local problems on the huge global matrix. The paper is organized as follows. The mathematical model is briefly presented in Section 2. The domain decomposition method and the construction of local matrices for the linear problem are described in Section 3, and the nonlinear problem with the activity switch is described in Section 4. The adaptive domain decomposition algorithm is similar to the algorithm from our paper [15], but it contains some updates and improvements, and it is therefore covered in Section 3. Section 5 provides a simple numerical example, and Section 6 evaluates the properties of the algorithm on a selected case study. 2. Mathematical model 2.1. Governing equations and their numerical treatment The problem of Darcian flow in a variably saturated porous medium is usually expressed by the Richards equation [1]. The governing partial differential equation and its finite element approximation will be described here in a brief only. A deep description of this model and its numerical treatment have already been given by the authors in a previous paper [15]. This section will maintain only the solution of the scattered system of linear equations after applying the domain decomposition on the entire computational domain Ω . Thus, in brief, the mathematical problem to be solved is the Richards equation in h-based form, where h is a so-called pore hydrostatic pressure head [L]. Several authors, e.g. [16], suggest avoiding this form due to inaccurate mass consistency in its numerical approximation. However, some later works by Tocci [17] and also by Kuraz et al. [7] suggest temporal adaptivities that preserve the mass balance. The Richards equation in the h-form is presumed as
∂h (x, t ) ∈ Ω × [0, T ), (1) ∂t where h is the capillary pressure head function [L], K(h) is the unsaturated hydraulic conductivity tensor of the second order θ (h) [L · T −1 ], C (h) is the water retention capacity [L−1 ], usually defined as C (h) = θ (h)′ + θ Ss , where θ (h) is the water content S function [−], Ss is the specific aquifer storage [L−1 ], θS is the saturated water content [−], and K(:,n) (h) is a vector formed out of the last column of the hydraulic conductivity tensor, where index n denotes the dimension of ℜn . A constitutive relation for function θ (h) was supplied by van Genuchten’s law [18], and for the K(h) function it was supplied by Mualem’s law [19]. ∇ · (K(h)∇ h) + ∇ · (K(:,n) (h)) = C (h)
Paper [15], mentioned above, explains in detail the finite element discretization of (1) and the Picard nonlinear iterations that lead to solving a system of linear equations Ax = b at each iteration. Because the problem (1) is nonlinear, we have to iteratively solve the problem A(x)x = b(x).
(2)
The Picard iterations perform a simple linearization of nonlinear product terms by moving all of them but one to the previous iteration level. And thus the method iteratively solves A(xk )xk+1 = b(xk ), where k is the iteration level, until
∥xk+1 − xk ∥2 < ε,
(3)
4
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
where ε is the desired iteration criterion. Matrix A is nonsymmetric, and we have already obtained results of good quality by applying conjugate gradient methods for normal equations combined with simple diagonal scaling. This was discussed by the authors in [20], and thus we solve the problem AT (xk )A(xk )xk+1 = AT (xk )b(xk ).
(4)
3. Domain decomposition The nonlinear problem (4) was solved here by a method that is close to alternating Schwarz domain decomposition. The method was applied in the following way: 1. consider index vectors Ii , where the subscript i = {1, . . . , nd}, where nd is the number of domains. Components of vector Ii define node indexes of subdomain Ωi ; 2. consider matrices Ni , where Ni ∈ ℜni ×n , where ni is the number of components of vector Ii (number of degrees of freedom of the subdomain i), and n is the number of degrees of freedom of the entire problem (3). Then the only nonzero elements are (Ni )j,jj = 1, where (Ii )j = jj ;
3. then local matrices on the subdomain are constructed out of the global matrix as Ai = A(:, Ii ), thus Ai = ANTi . And finally, in order to obtain symmetric local matrices, we have to perform the following: Bi = Ni AT ANTi .
(5)
The solution of (3) or matrix (4) is substituted by an algorithm solving the progression of matrices Bi , where i = {1, . . . , nd}. This measure, together with the adaptive domain decomposer, has significantly improved the performance of the nonlinear Picard solver, as explained in the following sections. 3.1. Algorithm The domain decomposition method was implemented in an algebraic way with Picard’s updates of the solution vector (deeply explained in Section 4). In our implementation, the Schwarz iterative procedure reduces the global residuum r = b − Ax. Calculating this global residuum is one of the most demanding computer operations. The iterative procedure is then as follows: 1. 2. 3. 4. 5. 6. 7.
calculate the global residuum; enter the loop across the subdomains; solve the problems Bi ci = ATi r (ATi Ai ci = ATi r); vector c is a vector of the updates for solution x, and thus xk+1 = xk + c, where k is the Schwarz iteration level; leave the subdomain loop; update the global matrix and go to step 1; the iterative loop is exited if the L2 norm of the global residuum r drops below the defined threshold.
In the case of the nonlinear Richards equation the problem becomes more complicated. The nonlinearities in certain parts of the domain require an excessive number of linear equation solver calls. The section below explains an application of this algorithm, and a method for the sequential deactivation of the solved subdomains. 4. Schwarz–Picard algorithm This section aims to provide a deep explanation of the local nonlinear operator treatment (5). The nonlinear problem (3) is solved in the following order. Roth time integration requires solving the nonlinear problem at each time level, and thus 1. at each time level, several Picard’s iterations must be performed—outer iterations. 2. for each Picard’s iteration we should perform several Schwarz iterations—inner iterations. 3. for each subdomain iterative linear equation solver (conjugate gradient method) has to perform several iterations— inner–inner iterations. In our implementation, only a single Schwarz iteration is performed within the Picard iteration level. The inner iteration corrects the solution on subdomains, the corrected solution is used to reconstruct matrix A, and the reconstructed matrix A is then used again for evaluating the global residuum. And thus the inner and outer iterations are combined together. The iteration algorithm can be partly considered as a block Jacobi method together with Picard’s updates of the solution vector. The following code demonstrates the implementation of the method.
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
1 3 5 7 9 11
5
time_integ : do ! prepare the data for solving the new time level call add_time_level () ! outer iteration schwarz_picard : do ! calculate the residuum r = b - matmul (A,x) ! inner iteration do i=1, ubound (sub ,1) ! check the subdomain activity and the local residuum if (sub(i)% active .or. norm2 (sub(i)%I) > r_err ) then
13
!add iteration level for each subdomain sub(i)% itcount = sub(i)% itcount (i) + 1
15
!get local matrix from the global matrix call get_submatrix (sub(i)%A, A)
17 19
! solve rectangular matrix (inner -- inner iteration ) call solve_matrix (sub(i)%A, sub(i)%c, r)
21
error = norm2 (sub(i)%c) 23
! update the solution x(sub(i)%I) = x(sub(i)%I) + sub(i)%c
25 27
if ( error < picard_error ) then sub(i)% active = . false . end if end if end do
29 31 33
! assemble the matrix A and vector b from the new results call assemble_matrix (A)
35 37 39
if ( all_domains_solved () .and. norm2 (r) < r_err ) then EXIT schwarz_picard end if end do schwarz_picard end do time_integ This algorithm will be more deeply described:
-line 1: -line 3: -line 5: -line 7:
enter the loop over time-steps prepare data for the new time level enter the Schwarz–Picard loop calculate the global residuum for the previous iteration level r = b(xk ) − A(xk )xk ,
the residuum is assumed for the global matrix, but can be also constructed out of the local matrices. -line 9: enter the subdomain loop and continue with the following: 1. (line 11) if the subdomain is active or if its local residuum is greater than the limit continue 2. (line 17) select the local matrix Ai from the global matrix A 3. (line 19) solve the local problem Bi ci = ATi r (ATi Ai ci = ATi r) by the conjugate gradient method for normal equations (inner–inner iterations). The input is the local rectangular matrix (Ai ), and the global residual vector r, and the output is the local vector of corrections ci 4. (line 21) the local error is considered as L2 norm of the local correction vector 5. (line 25) the solution in the next iteration level is updated in the following way: xk+1 = xk +
nd i=1
NTi wi oi ,
where wi is the diagonal matrix of the weight coefficients defined as (wi )α,α = 1, if
6
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
Fig. 1. Left: the domain is discretized both by the coarse mesh and by the fine mesh. Right: coarse mesh elements aggregate fine mesh elements into so-called clusters—the smallest subdomain units.
• the index jj = (Ii )α is unique in index sets (it is an internal node) • the index jj = (Ii )α is on the subdomain overlap and its subdomain index i is maximal, thus: jj ∈ Ii &(∀β)(jj ∈ Iβ : β ≤ i ) In other cases (wi )α,α = 0. 6. (line 27) if the Picard error in the subdomain is below the defined threshold, the subdomain is deactivated, but its local residuum is checked at each iteration (see line 9). If its value increases above the defined threshold, the subdomain is solved again. The local residuum is obtained by ri = Ni r. -line 33: update global matrix A -line 35: the Schwarz–Picard loop is exited in the case that the Picard error for all subdomains together with the L2 norm of the global residuum ∥r∥2 is below the defined threshold. -line 37: if the Schwarz–Picard loop is exited the current time level is solved, and the algorithm can move forward for the next time level solution (see line 3). The algorithm updates the residual vector outside the loop across the subdomains, and thus it is a type of additive method. 4.1. Richards equation application and domain decomposer The problem of the numerical solution of the Richards equation has already been summarized in our previous paper [15]. It concerns:
• the poor conditioning of the system of linear equations arising from the discretization of the system (1); • the slow convergence, or even divergence, of the nonlinear operator. Thus to solve system (2) one needs to resolve the slow convergence of the iterative linear equation solver and the slow convergence of the nonlinear operator. The moving wetting front plays an essential role in this problem. An optimal approach for the domain decomposition should then be time-adaptive. The domain decomposition algorithm should not be strictly algebraic, but should be linked with the Picard nonlinear solver, so that the Picard solver iterates on small matrices only. Unlike in our previous paper [15], we can now handle this capability. The domain decomposition adaptivity algorithm has been updated in order to ensure that the disturbances to the solution are covered by the small matrices. It was concluded that method described in our previous work was not optimal in terms of the treatment of the nonlinear operator. 4.1.1. Adaptive domain decomposer Based on the assumptions introduced here, the perfect and optimal domain decomposition algorithm should split the computational domain in such a way, that the wetting front is kept apart from the rest of the subdomains. In order to provide the subdomains with a proper shape (we should avoid needle-shaped subdomains) the domain decomposer algorithm works with two-level discretization—a coarse mesh and a fine mesh (see Fig. 1). At this point the algorithm is identical with the algorithm described in our previous work, but in order to keep our text consistent we should provide this explanation again. It is assumed that the mesh algorithm discretizes the domain in a geometrically optimal way. The mesh algorithm creates the fine mesh—a mesh that is used for the finite element method, and the coarse mesh—a mesh with optimal geometrical properties with a much lower density than the fine mesh. Elements of the coarse mesh aggregate the fine mesh elements into so-called clusters. The subdomain is formed out of clusters that can be merged or split on the basis of the solution properties. It is the base unit for constructing the subdomains. A single cluster cannot be split. The adaptivity algorithm was assumed here in the following way:
• the initial procedure creates a graph of the coarse mesh, so that the record of the mesh structure – neighbor elements – is kept;
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
7
Fig. 2. Sketch of the adaptive domain decomposition algorithm.
• the adaptivity works with the previous time and iteration level; • based on the previous time level solution the subdomains are formed in the following way; – the solution hk−1 on the clusters is evaluated as follows: 2 4 6
do i=1, no_of_clusters if (grad(h) > G .or. cluster (i)% it > 1) then cluster (i)% treat_alone = .true. else cluster (i)% treat_alone = . false . end if end do ∗ if the solution in the previous time level contains gradients steeper than G, or if the solution at this particular cluster required more than a single iteration to converge, then this cluster forms a single subdomain and cannot be joined with any other cluster; ∗ if the error of the Picard iteration at the first iteration is greater than the iteration criterion, the subdomains are re-evaluated, and this particular cluster again cannot be joined with any other cluster; – the Picard iteration error is measured for each cluster separately
1 3
error = 0.0 do ii=1, no_clusters_in_domain error = maxval (error , norm2 ( c_glob ( cluster (ii )%I))) end do ∗ where vector c is the vector of corrections for the next iteration level the vector component I of the cluster structure is the integer array of indexes of the fine mesh nodes inside the cluster ii, and the function norm2 is defined as a vector L2 norm (Fortran 2008 standard).
The adaptivity algorithm will be further explained by the following simple example—see Fig. 2, which demonstrates a source of water in the top left corner. The wetting front initially covers only a single cluster, and thus the domain decomposition algorithm decomposes the domain into a small domain of cluster size and into a large domain, which was created by merging the rest of the clusters outside of the wetting front – see Fig. 2 – left. Later, the wetting front intruded deeper into the domain, and thus the adaptive algorithm decomposed the domain into 6 subdomains: two big subdomains covering the area above and below the wetting front, and four subdomains of cluster size, covering the wetting front. 5. Numerical experiment This section presents a simple numerical example in order to evaluate the proposed algorithm. Because we would like to keep some continuity with our previous paper [15], an identical case study is presented here. The current results,
8
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
Fig. 3. Scheme of the numerical experiment setup. Table 1 Table of the unsaturated hydraulic parameters—parameters of the van Genuchten’s constitutive functions, see [18,19]. Ks (cm d−1 )
θr (–)
θs (–)
αvg (cm−1 )
nv g (–)
mv g (–)
Ss (cm−1 )
712.8
0.0
0.43
1.0
2.0
0.626
0.2
together with the previously published results, form a deep evaluation of the proposed algorithm. The scheme of this purely theoretical assignment is depicted in Fig. 3. The boundary condition setup is stated as: h(x, t ) = −10.0 cm, ∀(x, t ) ∈ ΓD × [0, T ), ∂h (x, t ) = 0, ∀(x, t ) ∈ ΓN × [0, T ). ∂ n⃗ The initial condition setup is stated as: h(x) = −200.0 cm,
∀x ∈ Ω .
The computational domain Ω was discretized by a uniform triangular fine mesh with 10 201 nodes and 20 000 elements, and by a uniform quadrilateral coarse mesh of 121 nodes and 100 elements. The material properties were assumed for a gravel material—a material with almost uniform pore sizes. Its constitutive functions are of a highly heterogeneous and steep nature. The parameters were obtained from Kuraz et al. [8], see Table 1. The computational time was discretized by a constant time step, despite the fact that the constant time step is inappropriate for optimal temporal discretization. Preserving the mass consistency is not crucial in this particular evaluation. It aims only to demonstrate the proposed domain decomposition adaptivity. The selected time step length was of constant size 1.33 × 10−2 days, the total simulation time was T = 786 days, and thus the uniform time step temporally discretized the solution into 59 098 time levels. 6. Results and conclusions The case study was solved both with a standard method, where each iteration of the Picard method required solving the entire system of linear equations, and also with our proposed algorithm. Obviously we obtained identical solutions. The standard method requires solving the system of linear equations with a dimension of the number of degrees of freedom for the entire domain discretized with a static mesh. The Picard method applied here only solves the nonlinearities (2), and thus due to the constant time step and the constant boundary condition setup it converged at an almost constant speed, and required three iterations to converge. The iteration count was within the range of ⟨3, 4⟩. However, the proposed domain decomposition technique is actually a combination of the block Jacobi method and the Picard nonlinear solver, and thus the algorithm obviously required more iterations than the standard approach to converge. The algorithm required on an average 4–5 iterations to converge, and the iteration count ranged in the interval of ⟨3, 8⟩, see Fig. 4. In the case of this proposed domain decomposition method, the iteration count is simply not global, it is scattered into subdomain iterations, and thus the iteration count mentioned here is the highest number of iterations for the most critical subdomain. The case study presented here was temporally discretized into 59 098 time levels, the domain decomposition algorithm changed the domain splitting 2685 times, and thus the subdomains were changed approximately after every 20th time level. As has been mentioned, the average number of iterations per time step was between 4 and 5, and thus on an average the mesh was constant for each 80–100 linear equation solver calls. The number of subdomains was initially two, exactly as expected in Fig. 2, and later the number of subdomains reached 34. As is shown in Fig. 5, the number of subdomains oscillated more or less around some optimal number. This actually means that the adaptive domain decomposer requires
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
9
Fig. 4. Number of iterations on critical subdomain for the first 200 days of the simulated time.
Fig. 5. Number of subdomains over the entire simulation time.
some further improvements, and that the subdomain distribution could be kept constant for a significantly greater period than 80–100 linear equation solver calls. See Fig. 6 for a solution and the subdomain distribution state at a simulation time of 550 days. Some of the most important results are the computational times for the proposed method compared to standard approaches. The computation was performed on Intel(R) Xeon(R) CPU E5520 2.27 GHz with memory speed 1333 MHz. The peak memory consumption of the code was 40 MB, and thus it easily fitted into the installed 8 GB of RAM. The code was compiled with gfortran gcc 4.7.2 with the full optimizations parameter. It is well known that nearly all Richards equation problems form ill-posed matrices in their numerical discretization. It was also discussed by the authors in [20] that due to the diagonal dominance, simple diagonal scaling greatly improves the properties of matrices of this type. Thus we solved the problem using two different approaches: 1. a standard approach with diagonal preconditioning; 2. the proposed domain decomposition approach, with no preconditioning applied to the local matrices. The first approach required 89 h of computer run-time, and the second approach required 34 h of computer run-time. It should be noted that the current version of the code is still highly non-efficient. In particular, the adaptivity requires a lot of memory allocations and deallocations that should be optimized together with the evaluation of the global residuum. All computations were performed serially on a single CPU core. It is assumed that further code optimizations will lead to huge savings in computer run time. The goal of this paper was to present the newly implemented algorithm for adaptive domain decomposition. The method is based on locating the nonlinearities in the solution, which requires a large number of linear equation solver calls to converge. The domain split adaptively sets up into such an order, that the problematic parts of the solution are always covered by matrices of small dimensions, so that the iterations are computationally inexpensive. Some partial results were already published in our paper [15], where the conditioning of the local matrices was studied. However, this current paper
10
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
Fig. 6. Solution, subdomain distribution, and number of iterations on subdomains at simulation time 550 days.
has provided a complete solution for the nonlinear problem, where the presented class of Richards equation problems on coarse porous media belongs. When this method is implemented in a stable manner, it will form a part of the DRUtES opensource FEM package [21]. Acknowledgment Financial support from the Czech Science Foundation (research project GACR 13-11977P) is gratefully acknowledged. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
L.A. Richards, Capillary conduction of liquids through porous media, Physics I 1 (1931) 318–333. S.P. Neuman, P.A. Witherspoon, Finite element method analyzing steady seepage with a free surface, Water Resour. Res. 6 (1970) 889–897. H.W. Alt, S. Luckhaus, Quasilinear elliptic–parabolic differential equations, Math. Z. 183 (1983) 311–341. J. Kačur, Solution of nonlinear and degenerate convection–diffusion problems, Nonlinear Anal. 47 (2001) 123–134. 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) (2008) 4610–4625. P. Šolín, M. Kuráž, Solving the nonstationary Richards equation with adaptive hp-FEM, Adv. Water Resour. 34 (2011) 1062–1081. M. Kuraz, P. Mayer, M. Leps, D. Trpkosova, An adaptive time discretization of the classical and the dual porosity model of Richards’ equation, J. Comput. Appl. Math. 233 (2010) 3167–3177. M. Kuraz, P. Mayer, V. Havlicek, P. Pech, J. Pavlasek, Dual permeability variably saturated flow and contaminant transport modeling of a nuclear waste repository with capillary barrier protection, Appl. Math. Comput. 219 (2013) 7127–7138. M. Dryja, O.B. Widlund, Towards a Unified Theory of Domain Decomposition Algorithms for Elliptic Problems, SIAM, Philadelphia, 1990, pp. 3–21. H.A. Schwarz, Gesammelete Mathematische Abhandlungen, Vol. 2, Springer, 1890, pp. 133–143. M. Benzi, A. Frommer, R. Nabben, D.B. Szyld, Algebraic theory of multiplicative Schwarz methods, Numer. Math. 89 (2001) 605–639. M. Dryja, O.B. Widlund, An Additive Variant of the Schwarz Alternating Method for the Case of Many Subregions, in: Ultracomputer Note, vol. 131, 1987. X.C. Cai, M. Dryja, Domain decomposition methods for monotone nonlinear elliptic problems, Contemp. Math. 180 (1994) 21–28. M. Dryja, W. Hackbusch, On the nonlinear domain decomposition method, BIT Numer. Math. 37 (1997) 296–311. M. Kuraz, P. Mayer, V. Havlicek, P. Pech, Domain decomposition adaptivityfor the Richards equation, Computing 95 (1) (2013) 501–519. M.A. Celia, E.T. Bouloutas, R.L. Zarba, A general mass-conservative numerical solution for the unsaturated flow equation, Water Resour. Res. 26 (1990) 1483–1496.
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
11
[17] 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, Adv. Water Resour. 20 (1997) 1–14. [18] M.T.H. van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, J. Soil Sci. 44 (1980) 892–898. [19] Y. Mualem, A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res. 12 (1976) 513–522. [20] M. Kuraz, P. Mayer, Algorithms for solving Darcian flow in structured porous media, Acta Polytech. 53 (2013) 347–358. [21] M. Kuraz, DRUtES—open source FEM library for coupled convection–difussion–reaction problems, http://freecode.com/projects/drutes, WWW.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Solving the nonlinear Richards equation model with adaptive domain decomposition Michal Kuraz a,∗ , Petr Mayer b , Pavel Pech a a
Czech University of Life Sciences Prague, Faculty of Environmental Sciences, Department of Water Resources and Environmental Modeling, Czech Republic b
Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mathematics, Czech Republic
article
info
Article history: Received 10 October 2013 Received in revised form 16 February 2014 Keywords: Richards equation Highly heterogeneous material properties Additive Schwarz domain decomposition Nonlinear operator treatment
abstract Modeling the transport processes in a vadose zone plays an important role in predicting the reactions of soil biotopes to anthropogenic activity, e.g. modeling contaminant transport, the effect of the soil water regime on changes in soil structure and composition, etc. Water flow is governed by the Richards equation, while the constitutive laws are typically supplied by the van Genuchten model, which can be understood as a pore size distribution function. Certain materials with dominantly uniform pore sizes (e.g. coarse-grained materials) can exhibit ranges of constitutive function values within several orders of magnitude, possibly beyond the length of real numbers that computers can handle. Thus a numerical approximation of the Richards equation often requires the solution of systems of equations that cannot be solved on computer arithmetics. An appropriate domain decomposition into subdomains that cover only a limited range of constitutive function values, and that will change adaptively, reflecting the time progress of the model, will enable an effective and reliable solution of this problem. Parts of this problem have already been treated in our work Kuraz et al. (2014). This paper focuses on improving the performance of a nonlinear solver by locating the areas with abrupt changes in the solution. © 2014 Elsevier B.V. All rights reserved.
1. Introduction The problem of predicting fluid movement in an unsaturated/saturated zone is important in many fields, ranging from agriculture, and hydrology to technical applications of dangerous waste disposal in deep rock formations. The mathematical model of unsaturated flow was originally published by L.A. Richards [1]. 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 [2]. The existence and the uniqueness of its solution were discovered 10 years later, by Alt and Luckhaus [3]. An analysis of the Richards equation problem has been the topic of several works by Kačur, e.g. [4]. A variety of methods for the numerical treatment of the Richards equation have been published in recent times, see e.g. [5–8]. A distinction is often made between those domain decomposition algorithms that use overlapping subregions and those which do not [9]. The alternating Schwarz method is the oldest, dating back to the 19th century [10], and it belongs to the first category (overlapping subregions). This method, originally designed as a tool for analysis, has been the subject of many
∗
Corresponding author. Tel.: +420 777018743. E-mail address: [email protected] (M. Kuraz).
http://dx.doi.org/10.1016/j.cam.2014.03.010 0377-0427/© 2014 Elsevier B.V. All rights reserved.
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
3
more or less recent publications and computer applications, see e.g. [11]. An additive variant of this algorithm was well treated in the 1980s e.g. by Dryja and Widlund [12]. Cai and Dryja [13] have studied the convergence of nonlinear strongly elliptic equations discretized by the finite element method when the Schwarz additive algorithms are taken into account, together with the classical inexact Newton method. Cai and Dryja provided a proof of the convergence of this algorithm. The convergence rate is independent of the finite element mesh parameter and the number of subdomains used in the domain decomposition. It was also latter proved by Dryja and Hackbusch [14], that the nonlinear operator converges locally with the same asymptotic speed as the corresponding linear iteration applied to the linearized problem. A numerical solution of the Richards equation typically encounters two major obstacles—the poor conditioning of the arising system of linear equations, and the slow convergence or even divergence of the nonlinear operator. The aim of this paper is to investigate a method that adaptively splits the computational domain, so that the nonlinear operator always iterates only on small matrices, as the convergence problem of the Richards equation is in most cases only a local problem around the wetting front. This paper extends and improves the results from the authors’ previous work [15], which studied the adaptive decomposer and the conditioning of the subdomain matrices. We can already present the results obtained from the algorithm that iteratively solves the nonlinearities on local matrices only. In our implementation, a combination of the Schwarz–Picard algorithm is considered rather than the Schwarz–Newton algorithm, as presented in [13]. We enable a technique that controls the iteration updates and the local residuum value, in order to localize the solution disturbances, so that the computational power can be focused on local problems locally, rather than on solving local problems on the huge global matrix. The paper is organized as follows. The mathematical model is briefly presented in Section 2. The domain decomposition method and the construction of local matrices for the linear problem are described in Section 3, and the nonlinear problem with the activity switch is described in Section 4. The adaptive domain decomposition algorithm is similar to the algorithm from our paper [15], but it contains some updates and improvements, and it is therefore covered in Section 3. Section 5 provides a simple numerical example, and Section 6 evaluates the properties of the algorithm on a selected case study. 2. Mathematical model 2.1. Governing equations and their numerical treatment The problem of Darcian flow in a variably saturated porous medium is usually expressed by the Richards equation [1]. The governing partial differential equation and its finite element approximation will be described here in a brief only. A deep description of this model and its numerical treatment have already been given by the authors in a previous paper [15]. This section will maintain only the solution of the scattered system of linear equations after applying the domain decomposition on the entire computational domain Ω . Thus, in brief, the mathematical problem to be solved is the Richards equation in h-based form, where h is a so-called pore hydrostatic pressure head [L]. Several authors, e.g. [16], suggest avoiding this form due to inaccurate mass consistency in its numerical approximation. However, some later works by Tocci [17] and also by Kuraz et al. [7] suggest temporal adaptivities that preserve the mass balance. The Richards equation in the h-form is presumed as
∂h (x, t ) ∈ Ω × [0, T ), (1) ∂t where h is the capillary pressure head function [L], K(h) is the unsaturated hydraulic conductivity tensor of the second order θ (h) [L · T −1 ], C (h) is the water retention capacity [L−1 ], usually defined as C (h) = θ (h)′ + θ Ss , where θ (h) is the water content S function [−], Ss is the specific aquifer storage [L−1 ], θS is the saturated water content [−], and K(:,n) (h) is a vector formed out of the last column of the hydraulic conductivity tensor, where index n denotes the dimension of ℜn . A constitutive relation for function θ (h) was supplied by van Genuchten’s law [18], and for the K(h) function it was supplied by Mualem’s law [19]. ∇ · (K(h)∇ h) + ∇ · (K(:,n) (h)) = C (h)
Paper [15], mentioned above, explains in detail the finite element discretization of (1) and the Picard nonlinear iterations that lead to solving a system of linear equations Ax = b at each iteration. Because the problem (1) is nonlinear, we have to iteratively solve the problem A(x)x = b(x).
(2)
The Picard iterations perform a simple linearization of nonlinear product terms by moving all of them but one to the previous iteration level. And thus the method iteratively solves A(xk )xk+1 = b(xk ), where k is the iteration level, until
∥xk+1 − xk ∥2 < ε,
(3)
4
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
where ε is the desired iteration criterion. Matrix A is nonsymmetric, and we have already obtained results of good quality by applying conjugate gradient methods for normal equations combined with simple diagonal scaling. This was discussed by the authors in [20], and thus we solve the problem AT (xk )A(xk )xk+1 = AT (xk )b(xk ).
(4)
3. Domain decomposition The nonlinear problem (4) was solved here by a method that is close to alternating Schwarz domain decomposition. The method was applied in the following way: 1. consider index vectors Ii , where the subscript i = {1, . . . , nd}, where nd is the number of domains. Components of vector Ii define node indexes of subdomain Ωi ; 2. consider matrices Ni , where Ni ∈ ℜni ×n , where ni is the number of components of vector Ii (number of degrees of freedom of the subdomain i), and n is the number of degrees of freedom of the entire problem (3). Then the only nonzero elements are (Ni )j,jj = 1, where (Ii )j = jj ;
3. then local matrices on the subdomain are constructed out of the global matrix as Ai = A(:, Ii ), thus Ai = ANTi . And finally, in order to obtain symmetric local matrices, we have to perform the following: Bi = Ni AT ANTi .
(5)
The solution of (3) or matrix (4) is substituted by an algorithm solving the progression of matrices Bi , where i = {1, . . . , nd}. This measure, together with the adaptive domain decomposer, has significantly improved the performance of the nonlinear Picard solver, as explained in the following sections. 3.1. Algorithm The domain decomposition method was implemented in an algebraic way with Picard’s updates of the solution vector (deeply explained in Section 4). In our implementation, the Schwarz iterative procedure reduces the global residuum r = b − Ax. Calculating this global residuum is one of the most demanding computer operations. The iterative procedure is then as follows: 1. 2. 3. 4. 5. 6. 7.
calculate the global residuum; enter the loop across the subdomains; solve the problems Bi ci = ATi r (ATi Ai ci = ATi r); vector c is a vector of the updates for solution x, and thus xk+1 = xk + c, where k is the Schwarz iteration level; leave the subdomain loop; update the global matrix and go to step 1; the iterative loop is exited if the L2 norm of the global residuum r drops below the defined threshold.
In the case of the nonlinear Richards equation the problem becomes more complicated. The nonlinearities in certain parts of the domain require an excessive number of linear equation solver calls. The section below explains an application of this algorithm, and a method for the sequential deactivation of the solved subdomains. 4. Schwarz–Picard algorithm This section aims to provide a deep explanation of the local nonlinear operator treatment (5). The nonlinear problem (3) is solved in the following order. Roth time integration requires solving the nonlinear problem at each time level, and thus 1. at each time level, several Picard’s iterations must be performed—outer iterations. 2. for each Picard’s iteration we should perform several Schwarz iterations—inner iterations. 3. for each subdomain iterative linear equation solver (conjugate gradient method) has to perform several iterations— inner–inner iterations. In our implementation, only a single Schwarz iteration is performed within the Picard iteration level. The inner iteration corrects the solution on subdomains, the corrected solution is used to reconstruct matrix A, and the reconstructed matrix A is then used again for evaluating the global residuum. And thus the inner and outer iterations are combined together. The iteration algorithm can be partly considered as a block Jacobi method together with Picard’s updates of the solution vector. The following code demonstrates the implementation of the method.
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
1 3 5 7 9 11
5
time_integ : do ! prepare the data for solving the new time level call add_time_level () ! outer iteration schwarz_picard : do ! calculate the residuum r = b - matmul (A,x) ! inner iteration do i=1, ubound (sub ,1) ! check the subdomain activity and the local residuum if (sub(i)% active .or. norm2 (sub(i)%I) > r_err ) then
13
!add iteration level for each subdomain sub(i)% itcount = sub(i)% itcount (i) + 1
15
!get local matrix from the global matrix call get_submatrix (sub(i)%A, A)
17 19
! solve rectangular matrix (inner -- inner iteration ) call solve_matrix (sub(i)%A, sub(i)%c, r)
21
error = norm2 (sub(i)%c) 23
! update the solution x(sub(i)%I) = x(sub(i)%I) + sub(i)%c
25 27
if ( error < picard_error ) then sub(i)% active = . false . end if end if end do
29 31 33
! assemble the matrix A and vector b from the new results call assemble_matrix (A)
35 37 39
if ( all_domains_solved () .and. norm2 (r) < r_err ) then EXIT schwarz_picard end if end do schwarz_picard end do time_integ This algorithm will be more deeply described:
-line 1: -line 3: -line 5: -line 7:
enter the loop over time-steps prepare data for the new time level enter the Schwarz–Picard loop calculate the global residuum for the previous iteration level r = b(xk ) − A(xk )xk ,
the residuum is assumed for the global matrix, but can be also constructed out of the local matrices. -line 9: enter the subdomain loop and continue with the following: 1. (line 11) if the subdomain is active or if its local residuum is greater than the limit continue 2. (line 17) select the local matrix Ai from the global matrix A 3. (line 19) solve the local problem Bi ci = ATi r (ATi Ai ci = ATi r) by the conjugate gradient method for normal equations (inner–inner iterations). The input is the local rectangular matrix (Ai ), and the global residual vector r, and the output is the local vector of corrections ci 4. (line 21) the local error is considered as L2 norm of the local correction vector 5. (line 25) the solution in the next iteration level is updated in the following way: xk+1 = xk +
nd i=1
NTi wi oi ,
where wi is the diagonal matrix of the weight coefficients defined as (wi )α,α = 1, if
6
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
Fig. 1. Left: the domain is discretized both by the coarse mesh and by the fine mesh. Right: coarse mesh elements aggregate fine mesh elements into so-called clusters—the smallest subdomain units.
• the index jj = (Ii )α is unique in index sets (it is an internal node) • the index jj = (Ii )α is on the subdomain overlap and its subdomain index i is maximal, thus: jj ∈ Ii &(∀β)(jj ∈ Iβ : β ≤ i ) In other cases (wi )α,α = 0. 6. (line 27) if the Picard error in the subdomain is below the defined threshold, the subdomain is deactivated, but its local residuum is checked at each iteration (see line 9). If its value increases above the defined threshold, the subdomain is solved again. The local residuum is obtained by ri = Ni r. -line 33: update global matrix A -line 35: the Schwarz–Picard loop is exited in the case that the Picard error for all subdomains together with the L2 norm of the global residuum ∥r∥2 is below the defined threshold. -line 37: if the Schwarz–Picard loop is exited the current time level is solved, and the algorithm can move forward for the next time level solution (see line 3). The algorithm updates the residual vector outside the loop across the subdomains, and thus it is a type of additive method. 4.1. Richards equation application and domain decomposer The problem of the numerical solution of the Richards equation has already been summarized in our previous paper [15]. It concerns:
• the poor conditioning of the system of linear equations arising from the discretization of the system (1); • the slow convergence, or even divergence, of the nonlinear operator. Thus to solve system (2) one needs to resolve the slow convergence of the iterative linear equation solver and the slow convergence of the nonlinear operator. The moving wetting front plays an essential role in this problem. An optimal approach for the domain decomposition should then be time-adaptive. The domain decomposition algorithm should not be strictly algebraic, but should be linked with the Picard nonlinear solver, so that the Picard solver iterates on small matrices only. Unlike in our previous paper [15], we can now handle this capability. The domain decomposition adaptivity algorithm has been updated in order to ensure that the disturbances to the solution are covered by the small matrices. It was concluded that method described in our previous work was not optimal in terms of the treatment of the nonlinear operator. 4.1.1. Adaptive domain decomposer Based on the assumptions introduced here, the perfect and optimal domain decomposition algorithm should split the computational domain in such a way, that the wetting front is kept apart from the rest of the subdomains. In order to provide the subdomains with a proper shape (we should avoid needle-shaped subdomains) the domain decomposer algorithm works with two-level discretization—a coarse mesh and a fine mesh (see Fig. 1). At this point the algorithm is identical with the algorithm described in our previous work, but in order to keep our text consistent we should provide this explanation again. It is assumed that the mesh algorithm discretizes the domain in a geometrically optimal way. The mesh algorithm creates the fine mesh—a mesh that is used for the finite element method, and the coarse mesh—a mesh with optimal geometrical properties with a much lower density than the fine mesh. Elements of the coarse mesh aggregate the fine mesh elements into so-called clusters. The subdomain is formed out of clusters that can be merged or split on the basis of the solution properties. It is the base unit for constructing the subdomains. A single cluster cannot be split. The adaptivity algorithm was assumed here in the following way:
• the initial procedure creates a graph of the coarse mesh, so that the record of the mesh structure – neighbor elements – is kept;
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
7
Fig. 2. Sketch of the adaptive domain decomposition algorithm.
• the adaptivity works with the previous time and iteration level; • based on the previous time level solution the subdomains are formed in the following way; – the solution hk−1 on the clusters is evaluated as follows: 2 4 6
do i=1, no_of_clusters if (grad(h) > G .or. cluster (i)% it > 1) then cluster (i)% treat_alone = .true. else cluster (i)% treat_alone = . false . end if end do ∗ if the solution in the previous time level contains gradients steeper than G, or if the solution at this particular cluster required more than a single iteration to converge, then this cluster forms a single subdomain and cannot be joined with any other cluster; ∗ if the error of the Picard iteration at the first iteration is greater than the iteration criterion, the subdomains are re-evaluated, and this particular cluster again cannot be joined with any other cluster; – the Picard iteration error is measured for each cluster separately
1 3
error = 0.0 do ii=1, no_clusters_in_domain error = maxval (error , norm2 ( c_glob ( cluster (ii )%I))) end do ∗ where vector c is the vector of corrections for the next iteration level the vector component I of the cluster structure is the integer array of indexes of the fine mesh nodes inside the cluster ii, and the function norm2 is defined as a vector L2 norm (Fortran 2008 standard).
The adaptivity algorithm will be further explained by the following simple example—see Fig. 2, which demonstrates a source of water in the top left corner. The wetting front initially covers only a single cluster, and thus the domain decomposition algorithm decomposes the domain into a small domain of cluster size and into a large domain, which was created by merging the rest of the clusters outside of the wetting front – see Fig. 2 – left. Later, the wetting front intruded deeper into the domain, and thus the adaptive algorithm decomposed the domain into 6 subdomains: two big subdomains covering the area above and below the wetting front, and four subdomains of cluster size, covering the wetting front. 5. Numerical experiment This section presents a simple numerical example in order to evaluate the proposed algorithm. Because we would like to keep some continuity with our previous paper [15], an identical case study is presented here. The current results,
8
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
Fig. 3. Scheme of the numerical experiment setup. Table 1 Table of the unsaturated hydraulic parameters—parameters of the van Genuchten’s constitutive functions, see [18,19]. Ks (cm d−1 )
θr (–)
θs (–)
αvg (cm−1 )
nv g (–)
mv g (–)
Ss (cm−1 )
712.8
0.0
0.43
1.0
2.0
0.626
0.2
together with the previously published results, form a deep evaluation of the proposed algorithm. The scheme of this purely theoretical assignment is depicted in Fig. 3. The boundary condition setup is stated as: h(x, t ) = −10.0 cm, ∀(x, t ) ∈ ΓD × [0, T ), ∂h (x, t ) = 0, ∀(x, t ) ∈ ΓN × [0, T ). ∂ n⃗ The initial condition setup is stated as: h(x) = −200.0 cm,
∀x ∈ Ω .
The computational domain Ω was discretized by a uniform triangular fine mesh with 10 201 nodes and 20 000 elements, and by a uniform quadrilateral coarse mesh of 121 nodes and 100 elements. The material properties were assumed for a gravel material—a material with almost uniform pore sizes. Its constitutive functions are of a highly heterogeneous and steep nature. The parameters were obtained from Kuraz et al. [8], see Table 1. The computational time was discretized by a constant time step, despite the fact that the constant time step is inappropriate for optimal temporal discretization. Preserving the mass consistency is not crucial in this particular evaluation. It aims only to demonstrate the proposed domain decomposition adaptivity. The selected time step length was of constant size 1.33 × 10−2 days, the total simulation time was T = 786 days, and thus the uniform time step temporally discretized the solution into 59 098 time levels. 6. Results and conclusions The case study was solved both with a standard method, where each iteration of the Picard method required solving the entire system of linear equations, and also with our proposed algorithm. Obviously we obtained identical solutions. The standard method requires solving the system of linear equations with a dimension of the number of degrees of freedom for the entire domain discretized with a static mesh. The Picard method applied here only solves the nonlinearities (2), and thus due to the constant time step and the constant boundary condition setup it converged at an almost constant speed, and required three iterations to converge. The iteration count was within the range of ⟨3, 4⟩. However, the proposed domain decomposition technique is actually a combination of the block Jacobi method and the Picard nonlinear solver, and thus the algorithm obviously required more iterations than the standard approach to converge. The algorithm required on an average 4–5 iterations to converge, and the iteration count ranged in the interval of ⟨3, 8⟩, see Fig. 4. In the case of this proposed domain decomposition method, the iteration count is simply not global, it is scattered into subdomain iterations, and thus the iteration count mentioned here is the highest number of iterations for the most critical subdomain. The case study presented here was temporally discretized into 59 098 time levels, the domain decomposition algorithm changed the domain splitting 2685 times, and thus the subdomains were changed approximately after every 20th time level. As has been mentioned, the average number of iterations per time step was between 4 and 5, and thus on an average the mesh was constant for each 80–100 linear equation solver calls. The number of subdomains was initially two, exactly as expected in Fig. 2, and later the number of subdomains reached 34. As is shown in Fig. 5, the number of subdomains oscillated more or less around some optimal number. This actually means that the adaptive domain decomposer requires
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
9
Fig. 4. Number of iterations on critical subdomain for the first 200 days of the simulated time.
Fig. 5. Number of subdomains over the entire simulation time.
some further improvements, and that the subdomain distribution could be kept constant for a significantly greater period than 80–100 linear equation solver calls. See Fig. 6 for a solution and the subdomain distribution state at a simulation time of 550 days. Some of the most important results are the computational times for the proposed method compared to standard approaches. The computation was performed on Intel(R) Xeon(R) CPU E5520 2.27 GHz with memory speed 1333 MHz. The peak memory consumption of the code was 40 MB, and thus it easily fitted into the installed 8 GB of RAM. The code was compiled with gfortran gcc 4.7.2 with the full optimizations parameter. It is well known that nearly all Richards equation problems form ill-posed matrices in their numerical discretization. It was also discussed by the authors in [20] that due to the diagonal dominance, simple diagonal scaling greatly improves the properties of matrices of this type. Thus we solved the problem using two different approaches: 1. a standard approach with diagonal preconditioning; 2. the proposed domain decomposition approach, with no preconditioning applied to the local matrices. The first approach required 89 h of computer run-time, and the second approach required 34 h of computer run-time. It should be noted that the current version of the code is still highly non-efficient. In particular, the adaptivity requires a lot of memory allocations and deallocations that should be optimized together with the evaluation of the global residuum. All computations were performed serially on a single CPU core. It is assumed that further code optimizations will lead to huge savings in computer run time. The goal of this paper was to present the newly implemented algorithm for adaptive domain decomposition. The method is based on locating the nonlinearities in the solution, which requires a large number of linear equation solver calls to converge. The domain split adaptively sets up into such an order, that the problematic parts of the solution are always covered by matrices of small dimensions, so that the iterations are computationally inexpensive. Some partial results were already published in our paper [15], where the conditioning of the local matrices was studied. However, this current paper
10
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
Fig. 6. Solution, subdomain distribution, and number of iterations on subdomains at simulation time 550 days.
has provided a complete solution for the nonlinear problem, where the presented class of Richards equation problems on coarse porous media belongs. When this method is implemented in a stable manner, it will form a part of the DRUtES opensource FEM package [21]. Acknowledgment Financial support from the Czech Science Foundation (research project GACR 13-11977P) is gratefully acknowledged. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
L.A. Richards, Capillary conduction of liquids through porous media, Physics I 1 (1931) 318–333. S.P. Neuman, P.A. Witherspoon, Finite element method analyzing steady seepage with a free surface, Water Resour. Res. 6 (1970) 889–897. H.W. Alt, S. Luckhaus, Quasilinear elliptic–parabolic differential equations, Math. Z. 183 (1983) 311–341. J. Kačur, Solution of nonlinear and degenerate convection–diffusion problems, Nonlinear Anal. 47 (2001) 123–134. 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) (2008) 4610–4625. P. Šolín, M. Kuráž, Solving the nonstationary Richards equation with adaptive hp-FEM, Adv. Water Resour. 34 (2011) 1062–1081. M. Kuraz, P. Mayer, M. Leps, D. Trpkosova, An adaptive time discretization of the classical and the dual porosity model of Richards’ equation, J. Comput. Appl. Math. 233 (2010) 3167–3177. M. Kuraz, P. Mayer, V. Havlicek, P. Pech, J. Pavlasek, Dual permeability variably saturated flow and contaminant transport modeling of a nuclear waste repository with capillary barrier protection, Appl. Math. Comput. 219 (2013) 7127–7138. M. Dryja, O.B. Widlund, Towards a Unified Theory of Domain Decomposition Algorithms for Elliptic Problems, SIAM, Philadelphia, 1990, pp. 3–21. H.A. Schwarz, Gesammelete Mathematische Abhandlungen, Vol. 2, Springer, 1890, pp. 133–143. M. Benzi, A. Frommer, R. Nabben, D.B. Szyld, Algebraic theory of multiplicative Schwarz methods, Numer. Math. 89 (2001) 605–639. M. Dryja, O.B. Widlund, An Additive Variant of the Schwarz Alternating Method for the Case of Many Subregions, in: Ultracomputer Note, vol. 131, 1987. X.C. Cai, M. Dryja, Domain decomposition methods for monotone nonlinear elliptic problems, Contemp. Math. 180 (1994) 21–28. M. Dryja, W. Hackbusch, On the nonlinear domain decomposition method, BIT Numer. Math. 37 (1997) 296–311. M. Kuraz, P. Mayer, V. Havlicek, P. Pech, Domain decomposition adaptivityfor the Richards equation, Computing 95 (1) (2013) 501–519. M.A. Celia, E.T. Bouloutas, R.L. Zarba, A general mass-conservative numerical solution for the unsaturated flow equation, Water Resour. Res. 26 (1990) 1483–1496.
M. Kuraz et al. / Journal of Computational and Applied Mathematics 270 (2014) 2–11
11
[17] 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, Adv. Water Resour. 20 (1997) 1–14. [18] M.T.H. van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, J. Soil Sci. 44 (1980) 892–898. [19] Y. Mualem, A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res. 12 (1976) 513–522. [20] M. Kuraz, P. Mayer, Algorithms for solving Darcian flow in structured porous media, Acta Polytech. 53 (2013) 347–358. [21] M. Kuraz, DRUtES—open source FEM library for coupled convection–difussion–reaction problems, http://freecode.com/projects/drutes, WWW.