Benchmarking FEniCS for mantle convection simulations

Computers & Geosciences 50 (2013) 95–105

Contents lists available at SciVerse ScienceDirect

Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo

Benchmarking FEniCS for mantle convection simulations L. Vynnytska a,n, M.E. Rognes b, S.R. Clark a a b

Computational Geoscience Department, Simula Research Laboratory, P.O. Box 134, 1325 Lysaker, Norway Center for Biomedical Computing, Simula Research Laboratory, P.O. Box 134, 1325 Lysaker, Norway

a r t i c l e i n f o

abstract

Article history: Received 16 January 2012 Received in revised form 11 May 2012 Accepted 12 May 2012 Available online 22 May 2012

This paper evaluates the usability of the FEniCS Project for mantle convection simulations by numerical comparison to three established benchmarks. The benchmark problems all concern convection processes in an incompressible fluid induced by temperature or composition variations, and cover three cases: (i) steady-state convection with depth- and temperature-dependent viscosity, (ii) timedependent convection with constant viscosity and internal heating, and (iii) a Rayleigh–Taylor instability. These problems are modeled by the Stokes equations for the fluid and advection–diffusion equations for the temperature and composition. The FEniCS Project provides a novel platform for the automated solution of differential equations by finite element methods. In particular, it offers a significant flexibility with regard to modeling and numerical discretization choices; we have here used a discontinuous Galerkin method for the numerical solution of the advection–diffusion equations. Our numerical results are in agreement with the benchmarks, and demonstrate the applicability of both the discontinuous Galerkin method and FEniCS for such applications. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Finite element method Mantle convection Incompressibility FEniCS Advection–diffusion Discontinuous Galerkin

1. Introduction Simulations of Earth’s mantle are as important today as ever in understanding the possible and likely composition of the mantle and its governing physical laws. Since the 1980s, numerous software projects have arisen that can model the mantle. Some of these were purpose written for the problem; coordinate systems and numerical solvers are built into the code as in the CitcomS family (Zhong et al., 2000; Tan et al., 2006). Other projects have arisen with greater flexibility in mind (Gale and Underworld, Moresi et al., 2003). However, in all cases, modifying the underlying system of equations is not trivial. In this paper, we present benchmarks of a mantle simulation code for which the underlying physics are accessible and easily modified; this code utilizes the FEniCS Project framework as its base. The FEniCS Project is a modern collection of open source software components directed at the automated solution of differential equations by finite element methods (Logg et al., 2012; The FEniCS Team, 2011). The FEniCS components combine a high-level user interface (Alnæs, 2012), mimicking the mathematical formulation of a finite element discretization, with sophisticated computational techniques (Kirby, 2004; Kirby and

n

Corresponding author. E-mail addresses: [email protected] (L. Vynnytska), [email protected] (M.E. Rognes), [email protected] (S.R. Clark). 0098-3004/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cageo.2012.05.012

Logg, 2006), thus allowing a large class of finite element discretizations to be concisely defined while providing highly efficient code (Kirby and Logg, 2007; Ølgaard and Wells, 2010; Logg and Wells, 2010). In particular, the user interface allows for flexibly changing the equations and their method of discretization. This flexibility makes the project particularly interesting for fields in which computational experimentation with different constitutive relations, different auxiliary equations and different regimes are key. Its automation focus also allows for the automation of tedious and error-prone tasks such as the assembly of finite element matrices and vectors, thus encouraging rapid prototyping and development. Overall, the user-specified equations and methods of discretization remain explicit while the generic finite element implementation details are kept under the hood. The main aims of this paper are (a) to implement a set of mantle convection related benchmarks in FEniCS, (b) to compare the results with the references established in the literature and (c) to evaluate the suitability of FEniCS for geodynamical simulations. To this end, we consider three test examples: steady-state convection with temperature- and depth-dependent viscosity; time-dependent convection with constant viscosity and internal heating; and a Rayleigh–Taylor instability. The first two examples are established benchmarks, allowing us to estimate the accuracy of the numerical solution for the temperature-driven convection of an incompressible fluid (Blankenbach et al., 1989). The third example is widely used for the validation of numerical models that simulate compositionally heterogeneous flow (van Keken et al., 1997; Tackley and King, 2003).

96

L. Vynnytska et al. / Computers & Geosciences 50 (2013) 95–105

Mantle convection is typically modeled as a Stokes problem driven by thermal or compositional heterogeneities. The transfer of heat or composition, in turn, is modeled by advection–diffusion equations. Most existing software tools rely on stabilized upwind Petrov–Galerkin (SUPG) approaches utilizing continuous finite element spaces for solving these types of equations (cf. King et al., 2010 and references therein). Here, we instead utilize a discontinuous Galerkin (DG) method. DG methods have grown in popularity over the last decades (Cockburn et al., 1999; Arnold et al., 2002; Zienkiewicz et al., 2003). By construction, the methods fit the finite element framework, in particular the FEniCS components natively support DG methods (Ølgaard et al., 2008), while discontinuities in data and solutions are easily represented. Moreover, the resulting discretizations can often be related to finite volume methods and often demonstrate advantageous conservation properties (Cockburn and Shu, 2001). Part of the purpose of this paper is therefore also to evaluate this method for these convection problems, and in particular to compare the results obtained with previous results obtained using SUPG approaches. This paper is structured as follows. We begin with a brief overview of the governing equations in Section 2, followed by the specification of the benchmark problems in Section 3. Next, Section 4 deals with the discretization methods and solution algorithms, while some aspects of the implementation are commented upon in Section 5. The numerical results obtained are compared to established results and discussed in Section 6. Finally, some concluding remarks are provided in Section 7.

2. Governing equations for mantle convection

@T þu  rTr  ðkT rTÞ ¼ q, @t

ð4Þ

where q is the rate of internal heat production and kT is the coefficient of thermal diffusion. In the steady-state regime and under the assumption of no internal heating, the left-most and right-most terms of (4) vanish, yielding the steady energy equation u  rTr  ðkT rTÞ ¼ 0:

ð5Þ

Similarly, the evolution of the composition f with relative chemical diffusivity kC, advected by the velocity field u, is governed by the equation (van Keken et al., 1997; Hansen and Yuen, 1988) @f þ u  rfr  ðkC rfÞ ¼ 0: @t

ð6Þ

This equation is used in the simulation of so-called thermochemical convection in the lower mantle where seismological studies may indicate compositional variations (for details see e.g. Deschamps et al., 2007; Steinberger and Holme, 2008). Note that the hyperbolic case kC ¼0 is relevant (Tackley and King, 2003). Additionally, boundary conditions must be supplied for the velocity or the stress, for the temperature and for the composition. Moreover, in the time-dependent case, initial conditions must be provided for the temperature and the composition.

3. Specification of benchmark problems

Over long timescales, mantle rock behaves like a momentumless fluid, driven instantaneously by buoyancy differences created by temperature or compositional variations. The equations considered here have previously been presented as representative of mantle convection and are used, with small changes, pervasively throughout the geodynamic community. The flow of an incompressible Stokes fluid over a domain O with negligible Reynolds number induced by temperature and composition variations is modeled by the following nondimensional system of (van Keken et al., 1997; Hansen and Yuen, 1988; Schubert et al., 2001) for each t A ½0,t end , x A O r  srp ¼ ðRb fRa TÞe,

ð1Þ

r  u ¼ 0,

ð2Þ

here the principal unknowns are the velocity u ¼ uðx,tÞ, the pressure p ¼ pðx,tÞ, the compositional field f ¼ fðx,tÞ, and the temperature T ¼ Tðx,tÞ, while Ra and Rb are the thermal and compositional Rayleigh numbers, respectively, and e is a unit vector in the direction of gravity. Note that if there is no thermal contrast, Ra ¼ 0, while if there is no compositional contrast, Rb ¼ 0. The Rayleigh numbers Ra and Rb are defined in terms of other physical parameters as follows (see for instance van Keken 3 3 et al., 1997): Ra ¼ arg DThr =kT Zr , Rb ¼ Drghr =kT Zr , where g is the acceleration of gravity, hr is the thickness, kT is the thermal diffusivity, Zr is the reference viscosity, a is the coefficient of thermal expansion, r is the density, Dr is the density difference between layers, and DT is the thermal contrast. Moreover, s ¼ sðuÞ is the deviatoric stress tensor. In the Newtonian case, the stress–strain relation takes the form

s ¼ 2Ze_ ðuÞ, e_ ðuÞ ¼ 12ðru þ ðruÞT Þ,

The evolution of the temperature field T is linked to the flow of the fluid, in particular the velocity field u, by the thermal advection–diffusion equation

ð3Þ

where e_ ðuÞ is the strain rate tensor and Z denotes the dynamic viscosity.

This study considers three benchmark problems originating from Blankenbach et al. (1989) and van Keken et al. (1997). The three problems concern 1. Steady-state thermal convection. 2. Time-dependent isoviscous thermal convection. 3. Time-dependent isothermal compositional convection. Each problem is described in the respective subsections below, while the aforementioned references provide more details. All problems consider a fluid contained in a two-dimensional box of length l and height H: O ¼ ½0, l  ½0,H with coordinates x ¼ ðx1 ,x2 Þ. 3.1. Steady-state thermal convection The first benchmark considers a steady-state thermal convection problem in a compositionally homogeneous fluid, as presented by Blankenbach et al. (1989). The unknown fields are the velocity u, the pressure p and the temperature T only, satisfying (1), (2) and (5) with Rb ¼ 0 and s defined by (3). The viscosity in (3) is assumed to depend on the depth and the temperature in the following manner:

Z ¼ Z0 expðbT=DT þ cðHx2 Þ=HÞ,

ð7Þ

where Z0 is a reference viscosity, DT is the thermal contrast, and b and c are additional given parameters. Free slip boundary conditions are given as un 9@O ¼ snt 9@O ¼ 0,

ð8Þ

where n denotes the outward normal and t the tangent vector on the boundary @O. Further, the temperature is fixed on the top and bottom boundaries, while symmetry conditions apply on the left

L. Vynnytska et al. / Computers & Geosciences 50 (2013) 95–105

layers is defined by the curve

and right boundaries T9x2 ¼ 0 ¼ DT,

T9x2 ¼ H ¼ 0,

97

@x1 T9x1 ¼ 0 ¼ @x1 T9x1 ¼ l ¼ 0:

ð9Þ

Three computational quantities are of particular interest: the root-mean-square velocity  1=2 Z 1 urms ¼ JuJ2 dx , ð10Þ lH O where dx denotes integration with respect to both coordinates; the point-wise temperature gradients at the corners of the domain

f ðx1 Þ ¼ Hd þ k cosðpx1 =lÞ,

ð17Þ

where Hd is a given height, and where k is an amplitude perturbation parameter. At the first time step t 1 ¼ Dt 1 (defined by the numerical resolution), the growth rate of instability g is estimated by

g ¼ ðln urms ðt 1 Þln urms ð0ÞÞ=t1 :

ð18Þ n

Moreover, at each time step t , the root-mean-square velocity urms ðt n Þ is computed according to (10).

x1 ¼ @x2 Tð0; 0Þ, x2 ¼ @x2 Tðl,0Þ, 4. Numerical methods

x3 ¼ @x2 Tð0,HÞ, x4 ¼ @x2 Tðl,HÞ;

ð11Þ

and the Nusselt number (following Blankenbach et al., 1989) Rl 0 @x2 T9x2 ¼ H dx1 Nu ¼  R l : ð12Þ 0 T9x2 ¼ 0 dx1 3.2. Time-dependent isoviscous thermal convection The second benchmark problem considers the time-dependent thermal convection of a fluid with constant viscosity undergoing internal heating. The governing equations are (1) with Rb ¼ 0, (2) and (4). No slip boundary conditions apply for the velocity on the top and bottom of the domain, while symmetry conditions apply on the side walls u9x2 ¼ 0,x2 ¼ H ¼ 0,

u1 9x1 ¼ 0,x1 ¼ l ¼ s12 9x1 ¼ 0,x1 ¼ l ¼ 0:

ð13Þ

The temperature is fixed at the top of the box, while no exchange of heat is allowed on the remainder of the boundary G ¼ @O\fx : x2 ¼ Hg T9x2 ¼ H ¼ 0,

@n T9G ¼ 0:

ð14Þ

This section describes general algorithms for solving the listed benchmark problems as well as spatial and temporal discretization methods for the governing equations. 4.1. Steady-state thermal convection algorithm The stationary, nonlinear system described in Section 3.1 can be solved by a Picard algorithm: iterating between solving (1) and (2) for the velocity and pressure with a given temperature, and solving (4) with a given velocity until convergence. Algorithm 1 lists the precise scheme considered, while the discrete formulations of the subproblems are summarized below. Algorithm 1. A Picard iteration algorithm for steady-state convection. 1 Make an initial guess for T 0h . Set initial velocity u0h ¼ 0. 2 Compute u1h by solving (19) with f ¼ ðRa T 0h Þe. 3 Compute T 1h by solving (25) with w ¼ u1h . 4 Set k¼2, E ¼ 106 . k2 k1 5 while Juk1 h uh J=Juh J 4 E:

The energy Eq. (4) is augmented by the initial condition Tðx,0Þ ¼ T 0 ðxÞ where the initial temperature distribution is obtained as perturbation of the static solution (Lennie et al., 1988)

6 7

þ T k2 (1) Set T mid ¼ 0:5ðT k1 h h Þ. (2) Solve (19) with f ¼ ðRa T mid Þe to obtain ukh .

8

(3) Set umid ¼ 0:5ðukh þ uhk1 Þ.

T 0 ðx1 ,x2 Þ ¼ 0:5ð1x22 Þ þ0:01 cosðpx1 =lÞ sinðpx2 =HÞ:

9 (4) Solve (25) with w ¼ umid to obtain T kh . 10 (5) k ¼ k þ 1. 11 end while

ð15Þ

For each scenario, on a set of discrete times, the root-meansquare velocity urms and the Nusselt number Nu are computed according to (10) and (12). 3.3. Time-dependent isothermal compositional convection (Rayleigh–Taylor instability) The third benchmark problem considers a time-dependent isothermal Rayleigh–Taylor instability, as formulated in van Keken et al. (1997). The domain is assumed to contain a fluid with two compositionally distinct layers: a buoyant lower layer with density fbtm and viscosity Zbtm and an upper layer with density ftop and viscosity Ztop . The unknowns are the velocity u, the pressure p and the composition f satisfying (1), (2) and (6). The boundary conditions (13) are prescribed for the velocity and the stress. In addition, if kC 40, we include boundary conditions for the composition f: its values are fixed at the top and bottom boundaries, and no flux conditions apply on the remainder of the boundary

f9x2 ¼ 0 ¼ fbtm ¼ 1, f9x2 ¼ H ¼ ftop ¼ 0, @n f9x1 ¼ 0,x1 ¼ l ¼ 0:

ð16Þ

The initial distribution of the compositional field is modeled by a step function equal to ftop ¼ 0 in the upper layer and to fbtm ¼ 1 in the lower layer. The initial interface between the

The incompressible Stokes equations are discretized using a mixed variational formulation (Zienkiewicz and Taylor, 2000; Taylor and Hood, 1973). To this end, let Vh denote the space of continuous piecewise quadratic vector fields and Qh the space of continuous piecewise linears, both defined relative to a mesh T h of the domain O. Further, let V 0h ¼ fv A V h 9vn 9@O ¼ 0g. The standard mixed finite element formulation of (1) and (2), with the boundary conditions (8), then reads as follows: for a given body source f (specifically f ¼ ðRb fRa T Þe), find uh A V 0h and ph A Q h such that Z Z 2Ze_ ðuh Þ : e_ ðvh Þ þ ph r  vh þ qh r  uh dx ¼ f  vh dx ð19Þ O

O

vh A V 0h

and all qh A Q h . for all On the other hand, the numerical treatment of (advectiondominated) advection–diffusion equations, such as (4), (5), and (6) is nontrivial. There exists a large body of research on the development of efficient discretization schemes for such problems (Cockburn et al., 1999; Zienkiewicz and Taylor, 2000; Cockburn and Shu, 2001; Zienkiewicz et al., 2003; Li, 2006). Here, the discretization is performed by a DG scheme with an upwind numerical flux taken from di Pietro et al. (2006). This scheme is

98

L. Vynnytska et al. / Computers & Geosciences 50 (2013) 95–105

defined below using the following notation. For generic scalar fields bs and vector fields bv , define the jump 1  U, average fg and restriction operators on any internal edge, shared by elements K þ and K  with outward normals n þ and n , respectively, by þ



þ

1bs U ¼ bs n þ þ bs n , þ



þ

fbs g ¼ 12 ðbs þ bs Þ,



1bv U ¼ bv  n þ þ bv  n ,

ð20Þ



fbv g ¼ 12ðbv þ bv Þ,

ð21Þ

bsþ ¼ bs 9K þ , b bvþ ¼ bv 9K þ , b s ¼ bs 9K  , v ¼ bv 9K  :

ð22Þ

On boundary edges, the jump and average operators are defined as 1bs U ¼ bs n, fbs g ¼ bs ,

1bv U ¼ bv  n,

ð23Þ

fbv g ¼ bv :

ð24Þ

Let Gi denote the interior edges of the mesh and Ge denote the boundary edges. Let Ph be a finite element space of (discontinuous) piecewise linears over T h . The corresponding DG formulation of (5) reads: for a given velocity w, find T h A P h such that aA ðw; T h , ch Þ þaD ðT h , ch Þ ¼ 0

ð25Þ

for all ch A P h , where the discretization of the advection term reads XZ aA ðw; T h , ch Þ ¼  T h w  rch dx K AT h K

þ

XZ

e A Gi

þ

e

XZ

e A Ge e

ð26Þ

and the discretization of the diffusion term reads XZ kT rT h  rch dx aD ðT h , ch Þ ¼ X Z

e A Gi [Ge e

þ

XZ 

e A Gi e

fkT rT h g  1ch U ds

fkT rch g  1T h Uþ

akT hK

1T h U  1ch U

for all ch A P h , where aA and aD are as defined in (26) and (27), respectively. The corrector step can be viewed as applying a Crank–Nicolson scheme in time and the DG method in space to (4) and reads: with T k1 and uk1 given, and for given predicted velocity upr , find the h h h corrected temperature T hk1 satisfying Z k T h T k1 h ch dx þ 0:5ðaA ðupr ; T kh , ch Þ þ aD ðT kh , ch ÞÞ h Dtk O Z k1 k1 qch dx ð29Þ þ 0:5ðaA ðuk1 h ; T h , ch Þ þaD ðT h , ch ÞÞ ¼ O

for all ch A Ph . Each time step Dt k þ 1 , k ¼ 0; 1, . . . ,n1 is defined by the formula ð30Þ

where hmin ¼ minK A T h hK is the minimal cell size of the mesh and C CFL is a chosen positive number relating to a Courant–Friedrichs– Lewy (CFL)-type condition. 4.3. Rayleigh–Taylor instability algorithm

K AT h K

þ

The predictor step results from applying an implicit Euler scheme in time and the DG method in space to (4) and reads: with the temperature T k1 and velocity uk1 approximations at the h h previous time given, find the predicted temperature T pr satisfying h Z Z pr k1 T h T h pr pr ch dx þ aA ðuk1 qch dx ð28Þ h ; T h , ch Þ þ aD ðT h , ch Þ ¼ Dt k O O

Dtk þ 1 ¼ C CFL hmin =max9ukh 9,

ðw  1ch UfT h g þ 129w  n91ch U1T h UÞ ds w  1ch UfT h g ds,

11 (5) Compute new time step Dt k þ 1 by (30). 12 (6) Set t k þ 1 ¼ t k þ Dt k þ 1 and increase k ¼ k þ 1. 13 end while



ds, ð27Þ

where a is a chosen stabilization constant.

The Rayleigh–Taylor instability problem is computationally analogous to the time-dependent thermal convection problem. The equations are therefore also decoupled and solved using a splitting scheme. However, in contrast to the thermal problem, filtering (correcting for numerical dispersion) is applied to the compositional field directly instead of solving an additional variational problem. This algorithm and the filtering scheme are described in Lenardic and Kaula (1993). Since the viscosity varies with the composition field, Z is also mapped by linear interpolation

4.2. Time-dependent thermal convection algorithm

Zkh ¼ Zbtm þ fkh ðZtop Zbtm Þ:

The time-dependent benchmark problem described in Section 3.2 can be solved using a predictor–corrector based splitting scheme (van den Berg et al., 1993; Hansen and Ebel, 1988). The discrete formulations of the predictor and corrector steps are summarized below; the detailed algorithm is listed in Algorithm 2.

The complete numerical scheme is described in Algorithm 3.

Algorithm 2. A predictor–corrector algorithm for thermal convection. 1 Define the initial temperature T0 by (15). 2 Define V 1h ¼ fv A V h 9v satisfying (13)}.

Algorithm 3. A predictor–corrector algorithm for compositionaldriven convection. 0

1 Initialize f . 0

2 Compute initial velocity u0 by solving (19) with f . 3 Compute time step Dt 1 from u0 according to (30). 4 for k ¼ 1, . . . ,n. pr 5 (1) Solve (32) to obtain fh .

3 Compute an initial velocity u0h by solving (19) with

6

f ¼ ðRa T 0 Þe in V 1h  Q h . 4 Compute an initial time step Dt 1 by (30) with w ¼ u0h .

7

1

5 Set k¼1 and t ¼ Dt 1 . 6 while t k rt end . . 7 (1) Solve (28) for T pr h 8

(2) Solve (19) for upr with f ¼ ðRa T pr Þe in V 1h  Q h . h h

9

(3) Solve (29) for T kh . (4) Solve (19) for ukh

10

with f ¼ ðRa T kh Þe in V 1h  Q h .

ð31Þ

pr

k

(2) Filter predicted composition fh to obtain fh . (3) Map

fkh

to obtain viscosity Zkh . k

8 (4) Solve (19) with fh and Zkh as input to obtain ukh . 9 (5) Compute new Dt k þ 1 according to (30). 10 end for The discretization of the transport equation (6) is carried out by an implicit Euler scheme in time and the DG method in space, yielding the formulation: for given approximated composition fk1 and velocity uhk1 at the previous time, find the predicted h

L. Vynnytska et al. / Computers & Geosciences 50 (2013) 95–105 pr

composition fh satisfying Z

k1 fpr h fh

O

Dtk

pr ch dx þ aA ðuk1 ; fpr h , ch Þ þaD ðfh , ch Þ ¼ 0

ð32Þ

99

as the equations itself. Essential boundary conditions are applied by adding the constraints on the degrees of freedom in the resulting linear systems of equations.

for all ch A P h . In (32), aA and aD are as defined in (26) and (27), respectively, but with kT replaced by kC.

6. Results

5. Implementation

Here we present and analyze the results obtained for the benchmark problems posed in Section 3.

Each of Algorithms 1–3 has been implemented using DOLFIN 1.0 (Hake et al., 2011). The variational problems (19), (28), (29) and (32) are implemented directly by specifying the variational forms defining the equation. Each of the variational forms involves operators acting on basis functions or coefficients integrated over cells or edges. For (19), standard derivative operators such as gradients and divergences are invoked, while operators such as jumps and averages over internal edges are used for (26) and (27). Both the basis functions and the coefficients are defined on prescribed finite element spaces. For (19), the product space of continuous piecewise quadratic vector fields and continuous piecewise linear scalar fields are used. For (25), (28), (29) and (32), the element space of (discontinuous) piecewise linears is used for the temperature or composition approximations (at each discrete time). For the initial temperature guess T 0h in Algorithm 1, the L2 projection of (33) into Ph is used, similarly for (15) in Algorithm 2. The assembly of the resulting element matrices and vectors is handled automatically by DOLFIN. In particular, the integrals in the variational forms are evaluated using exact numerical quadrature. We remark that forms of arity 2 (bilinear forms resulting in element matrices), arity 1 (linear forms resulting in element vectors), and arity 0 (functionals resulting in scalars) are supported. Thus, functionals such as (10) and (12) are specified in the same language and automatically evaluated in the same manner

6.1. Steady-state thermal convection The objective of this test is to perform a quantitative analysis of our results for urms , xi , Nu with the results of Blankenbach et al. (1989). To this end, we consider five scenarios of the steady-state thermal convection problem by varying the parameters for the viscosity Z, the thermal Rayleigh number Ra and the length of the domain l. All scenarios share the common parameter values DT ¼ 1, kT ¼ 1, H¼1, Z0 ¼ 1:0, and q ¼0. The following initial guess is made for the temperature, which is a perturbation of a static purely diffusive solution Tðx1 ,x2 Þ ¼ 1x2 (Lennie et al., 1988) T 0 ðx1 ,x2 Þ ¼ ð1x2 Þ þ 0:05 cosðpx1 =lÞ sinðpx2 =HÞ:

ð33Þ

The meshes were here constructed by uniformly dividing the computational domain into m  n rectangles, and then dividing each rectangle into four triangles by the two diagonals. The resulting computed values for the Nusselt number, the rootmean-square velocity, and the values of the temperature gradient at the corners for the five different scenarios are presented in Table 1, along with the reference values from Blankenbach et al. (1989), and the relative percentage error (in parenthesis below the computed values) compared to these reference values. h-convergence implies that the difference between the computed and the reference values decreases toward zero as the mesh resolution is refined. For all quantities of interest and all five scenarios, we observe that this difference decreases as the mesh is

Table 1 Steady-state thermal convection results: computed quantities of interest for five different sets of parameter values, compared to reference values from Blankenbach et al. (1989) for each case. The number in parenthesis below each computed value gives the percentage relative difference; that is, the absolute difference of the computed and the reference value, divided by the reference value, and multiplied by 100. mn

Nu

urms

x1

x2

x3

x4

Case I: Ra ¼ 104 , l ¼ 1:0, b ¼ c¼ 0 80  80 4.8755 (0.182) 160  160 4.88216 (0.046) 200  200 4.88297 (0.029) Reference 4.88441

42.8664 (0.003) 42.8653 (0.001) 42.8652 (0.001) 42.8649

0.58953 (0.12) 0.589009 (0.03) 0.58894 (0.02) 0.58881

8.05631 (0.04) 8.05865 (0.01) 8.05892 (0.01) 8.05938

8.05635 (0.04) 8.05865 (0.01) 8.05892 (0.01) 8.05938

0.589551 (0.13) 0.589011 (0.03) 0.58894 (0.02) 0.58881

Case II: Ra ¼ 105 , l ¼ 1:0, b ¼c ¼0 80  80 10.4393 (0.900) 160  160 10.5099 (0.230) 200  200 10.5186 (0.147) Reference 10.5341

193.239 (0.013) 193.22 (0.003) 193.218 (0.002) 193.215

0.726027 (0.45) 0.724018 (0.18) 0.723613 (0.12) 0.722751

19.0312 (0.25) 19.0686 (0.06) 19.0726 (0.04) 19.0794

19.0322 (0.25) 19.0687 (0.06) 19.0727 (0.04) 19.0794

0.72637 (0.50) 0.724062 (0.18) 0.723635 (0.12) 0.722751

Case III: Ra ¼ 106 , l ¼ 1:0, b ¼ c ¼ 0 80  80 20.9926 (4.460) 160  160 21.7126 (1.183) 200  200 21.8056 (0.759) Reference 21.9725

834.692 (0.084) 834.111 (0.015) 834.062 (0.009) 833.99

0.861836 (1.75) 0.881698 (0.52) 0.881149 (0.45) 0.87717

45.0743 (1.94) 45.7802 (0.40) 45.8524 (0.24) 45.9642

45.0926 (1.90) 45.783 (0.39) 45.8539 (0.24) 45.9642

0.86787 (1.06) 0.882427 (0.60) 0.881528 (0.50) 0.87717

Case IV: Ra ¼ 104 , l ¼ 1:0, b ¼ lnð1000Þ, c ¼0 80  80 10.0203 (0.454) 160  160 10.0537 (0.122) 200  200 10.0581 (0.078) Reference 10.066

481.001 (0.118) 480.516 (0.017) 480.48 (0.010) 480.433

0.490599 (1.36) 0.500353 (0.60) 0.499969 (0.52) 0.49738

26.1727 (2.37) 26.6586 (0.56) 26.7162 (0.34) 26.8085

17.5007 (0.17) 17.5222 (0.05) 17.5255 (0.03) 17.5314

1.01203 (0.35) 1.00948 (0.10) 1.00914 (0.06) 1.00851

Case V: Ra ¼ 104 , l ¼ 2:5, b ¼ lnð16384Þ, c ¼ lnð64Þ 100  40 6.79347 (1.969) 171.528 (0.132) 200  80 6.88967 (0.581) 171.779 (0.014) Reference 6.9299 171.755

0.622902 (0.84) 0.620716 (0.49) 0.6177

14.0379 (0.92) 14.1392 (0.20) 14.1682

19.2556 (4.17) 18.4738 (0.06) 18.4842

0.175353 (1.17) 0.177719 (0.17) 0.17742

100

L. Vynnytska et al. / Computers & Geosciences 50 (2013) 95–105

Fig. 1. Time-dependent thermal convection: illustration of the P2 cycle in the Nusselt number (left) and the root-mean-square velocity (right) with Ra ¼ 216 000, mesh 144  96, a ¼ 4, C CFL ¼ 0:5. Within each cycle, we label the pairs of local extrema by S0 and S1: S0 corresponds to the largest local maximum and following minimum, S1 the next.

refined. On the finest mesh, the maximal relative error for all scenarios and all computed quantities is less than 0.76%, and is less than 0.1% for many quantities and scenarios. Thus, we conclude that the method presented here provides almost identical results to those of Blankenbach et al. (1989). For all scenarios, we observe that the difference in the rootmean-square velocity is generally small: the relative error is less than 0.15% even on the coarsest mesh for the most complex scenario with temperature and depth dependent viscosity (case V). On the finest meshes, this difference is in the range 0.001– 0.015% for the different scenarios. In contrast, the differences for the Nusselt number and the temperature gradients are higher, particularly for case III (constant viscosity, Ra ¼ 106 ) and case V (temperature and depth dependent viscosity). For instance, for case III, the relative error in the Nusselt number on the coarsest mesh is approximately 4.5% (but decreases to 0.8% on the finest mesh). For this problem, we tested the influence of the Rayleigh number Ra, the stability parameter a and mesh resolution on the magnitude of the difference. For the scenarios with constant viscosity, we observe that the differences in all the computed quantities increase as the Rayleigh number increases. Further numerical experiments reveal that the stability parameter a in (27) has a clear influence on the magnitude of the differences. The values reported here result from choosing a ¼ 3. Experiments with larger parameter (for instance a ¼ 10) yielded values following the same trends as reported above, however the magnitude of the differences were greater. Finally, we remark that the mesh configuration plays a certain role; in particular, the corner values of the temperature gradients are highly sensitive to the choice of mesh configuration. Additional numerical experiments with meshes formed by dividing the sub-rectangles into triangles by the right or left diagonals (instead of both diagonals) gave clearly differing values from the ones presented here, although difference between results obtained and the reference values still decrease as the meshes were refined. Moreover, for triangular meshes of a rectangular domain, corner points are typically shared between two cells. Since the approximated temperature gradients are discontinuous across cells, the point evaluation of such fields on cell boundaries is not immediately well-defined. Here, we projected the temperature gradient onto the space of continuous piecewise linear functions before performing the point-wise evaluation, and remark that such choices add an additional level of post-processing to the numerical approximations and the reported values.

6.2. Time-dependent convection with constant viscosity and internal heating Following Blankenbach et al. (1989), we examine this evolutionary problem with two particular values for Ra, namely, Ra ¼ 216 000 and Ra ¼ 218 000. The domain is defined by the parameters H ¼ 1:0, l ¼ 1:5. Further, q¼ 1.0, kT ¼1.0, and Z ¼ 1:0. Internal heating with an insulated lower boundary and a fixed temperature on the surface leads to the formation of a thermal instability in the top layer. The behavior of the system may be examined by changes in Nu which measures the heat flux on the surface, and urms . First, a layer of colder material is formed on the top resulting in a decreasing heat flux on the top boundary manifested by a decreasing Nu. This leads to the formation of a sinking blob and produces strong horizontal temperature gradients which in turn drive velocities to a maximum. When this material reaches the bottom boundary, urms approaches its minimum. At the same time the colder layer is thinning and Nu increases. Then a new blob is formed and the process repeats again. For lower Ra, every blob behaves in the same way giving a so-called P1 cycle. Higher values of Ra allows P2 and P4 cycles; that is, every second and fourth blobs produce the same pattern, respectively. Here we consider the particular values of Ra since a shift from P2 to P4 cycle lies between these Rayleigh numbers. For more details, see Lennie et al. (1988). The primary aim for the comparison is to examine whether a P2 cycle is obtained for the case Ra ¼ 216 000 and a P4 cycle is obtained for the case Ra ¼ 218 000. A secondary aim is to compare the output quantities with the benchmarking values presented in Blankenbach et al. (1989) for the case Ra ¼ 216 000. Definite benchmarking values are not provided for the case Ra ¼ 218 000, so we instead compare our values to those of the presumed most accurate discretization reported in Blankenbach et al. (1989): the results corresponding to acronym ‘‘Ha’’ on a 96  64 mesh therein. The simulations were run for different meshes of increasing sizes1 m  n, for different values of the DG stabilization parameter a (a ¼ 4; 10), and for different values of the adaptive time step parameter C CFL (C CFL ¼ 0:5,0:8). The expected cycles in the Nusselt number and the root-mean-square velocity were obtained for all combinations of these parameters, see Figs. 1 and 2 for a sample case.

1 Each mesh of size m  n was constructed by dividing the computational domain into m  n sub-rectangles and then dividing each rectangle into two triangles by the diagonal from the lower left to the upper right corner.

L. Vynnytska et al. / Computers & Geosciences 50 (2013) 95–105

101

Fig. 2. Time-dependent thermal convection: illustration of the P4 cycle in the Nusselt number (left) and the root-mean-square velocity (right) with Ra ¼ 218 000, mesh 144  96, a ¼ 4, C CFL ¼ 0:5. Within each cycle, we label the pairs of local extrema by S0, S1, S2, S3: S0 corresponds to the largest local maximum and following minimum, S1 the next pair, etc., and then repeating.

Table 2 Time-dependent thermal convection: computed values for different thermal Rayleigh number and different mesh sizes. Within each cycle, consisting of two pairs of local extrema in the case Ra ¼ 216 000 and four pairs of local extrema in the case Ra ¼ 218 000, the pairs are labeled by S0,S1, . . .: S0 corresponds to the largest local maximum and following minimum, S1 the next pair, etc., and then repeating. The period was computed by comparing the time at which the S0 maxima for the Nusselt number occurred for the last two computed complete cycles. The values labeled Reference/Comparison denote reference/comparison values from Blankenbach et al. (1989). Stage

Mesh

Ra ¼ 216 000, a ¼ 4, CFL ¼ 0:5 S0 48  32 96  64 144  96 Reference S1

Period

0.04821 0.04813 0.04807 0.04803

48  32 96  64 144  96 Reference

Ra ¼ 218 000, a ¼ 4, CFL ¼ 0:5 S0 48  32 96  64 144  96 Comparison

0.09583 0.09557 0.09551 0.09513

Nu

urms

max

min

max

min

7.3326 7.3598 7.3664 7.379

6.4149 6.4438 6.4526 6.468

60.5112 60.3981 60.3842 60.367

31.9962 31.9863 31.9812 31.981

7.1479 7.1749 7.1830 7.196

6.7310 6.7678 6.7778 6.796

57.5346 57.4577 57.4706 57.430

30.3306 30.3244 30.3210 30.320

7.3456 7.3735 7.3801 7.395

6.4412 6.4801 6.4885 6.503

61.0071 61.0095 60.9876 61.518

32.0921 32.0403 32.0377 32.195

S1

48  32 96  64 144  96 Comparison

7.1869 7.2278 7.2349 7.254

6.7701 6.8134 6.8230 6.839

58.4364 58.6746 58.6662 59.153

30.4249 30.4093 30.4066 30.512

S2

48  32 96  64 144  96 Comparison

7.3424 7.3690 7.3756 7.391

6.3960 6.4192 6.4278 6.439

60.3267 60.1055 60.0880 60.541

32.2782 32.2852 32.2808 32.449

S3

48  32 96  64 144  96 Comparison

7.1115 7.1290 7.1370 7.148

6.7271 6.7565 6.7665 6.777

56.6667 56.3569 56.3655 56.620

30.5027 30.5137 30.5109 30.624

The period of the cycle and the cycle extrema for the case a ¼ 4 and C CFL ¼ 0:5 are reported in Table 2 for the different meshes and compared to the reference values. Results obtained demonstrate that bifurcation point from P2 to P4 lies in the range between Ra ¼ 216 000 and Ra ¼ 218 000. For the case Ra ¼ 216 000, we see that all values are close to the reference values; that is, they agree to within at least two digits even on the coarsest meshes. Moreover, all values but one seem to converge, and toward values consistent with the error ranges reported by Blankenbach et al. (1989). The second local maximal value of urms displays less

conclusive behavior, but still agrees to within three digits. For the case Ra ¼ 218 000, we emphasize that the reference values are for qualitative comparison only and not established benchmark values. However, our results agree with these to within two digits on the finest meshes. Moreover, the values for the Nusselt number seem to converge as the mesh resolution is refined. The values for the root-mean-square velocity show less monotone behavior, but mainly only vary in the third significant digit for the maxima and the fourth digit for the minima as the meshes are refined.

102

L. Vynnytska et al. / Computers & Geosciences 50 (2013) 95–105

Finally, we note that the period and extreme values of the cycles are robust also with respect to the stabilization parameter a and the adaptive time step parameter C CFL . However, the different parameter choices can give significantly differing results in the early stages of the simulation, and also shift the onset of the cyclic behavior in time (cf. Fig. 3 and its caption).

(max urms ) and the point in time at which this value is achieved (tðmax urms Þ) converge as the mesh (and accordingly also the time step) is refined. Moreover, these values agree well with those

6.3. Rayleigh–Taylor instability Generally speaking the following techniques are used for the representation of chemically distinct layers: field-based methods, tracers and marker chains. The main disadvantage of the fieldbased methods is their continuous representation of the compositional field and, thus, smearing of the interface between distinct layers. Hence, tracers and marker chains were purposely designed to deal with this issue (van Keken et al., 1997; Tackley and King, 2003). In addition, finite element methods with continuous basis functions are typically unable to consider the case of zero chemical diffusivity: kC ¼ 0. However, DG discretizations allow this case, and may hence reduce artificial numerical diffusion. The results provided in this section demonstrate how the drawbacks of the field-based methods can be overcome by the usage of a DG scheme. The domain is described by l ¼ 0:9142, H¼1.0, k ¼ 0:02, and we let a ¼ 4:0. Meshes of sizes n  n (n ¼ 40; 80,160) were constructed as for the steady convection test case, but then adjusted to capture the initial interface given by (17) so to better represent the discontinuous initial composition. Following van Keken et al. (1997) and Tackley and King (2003), we investigate the case kC ¼0 corresponding to no chemical diffusion. Furthermore, for the sake of comparison with SUPG approaches in which some chemical diffusion must be prescribed (artificial or not), we also consider the cases kC ¼ 106 ,107 . We prescribe a 100-fold viscosity contrast between top and bottom layers. The resulting evolution of the root-mean-square velocity is plotted in Fig. 4 for the different choices of kC and C CFL . We observe that all three parameter choices give comparable results and that the common trend agrees with that of van Keken et al. (1997). The computed values and additionally the instability growth rate g are further examined in Table 3. In the primary case kC ¼0, C CFL ¼ 0:3, the peak root-mean-square velocity

Fig. 4. Rayleigh–Taylor instability: root-mean-square velocity urms versus time t for different choices of the chemical diffusivity kC and the time step parameter C CFL on meshes of size 160  160. Table 3 Rayleigh–Taylor instability: computed values for the instability growth rate g, peak value of the root-mean-square velocity (max urms ), and time at which the peak value was achieved (tðmax urms Þ) for different choices of the chemical diffusivity parameter kC, mesh sizes and C CFL . kC

C CFL

Mesh

g

max urms

tðmax urms Þ

0 0 0

0.30 0.30 0.30

40  40 80  80 160  160

0.0976 0.1018 0.1039

0.01140 0.01444 0.01458

57.21 52.11 51.55 53.71

106

0.50

40  40

0.0975

0.01133

106

0.50

80  80

0.1024

0.01401

50.86

106

0.50

160  160

0.1011

0.01362

49.70

107

0.50

40  40

0.1107

0.01145

54.53

107

0.50

80  80

0.1182

0.01430

51.64

107

0.50

160  160

0.1204

0.01451

51.14

Fig. 3. Time-dependent thermal convection: the evolution of the Nusselt number at different time ranges (early: left, later: right) with Ra ¼ 218 000, mesh 144  96 and varying a and C CFL . The different parameter choices all give comparable extrema values and same periodicity (P4) for the later stabilized cycle. However, the graphs differ considerably in the early stages of the simulation (left). Moreover, in this case, the difference in the stabilization parameter results in a small temporal shift of the cycle (red versus blue), while the change in the adaptive time step parameter (blue versus green) gives a shift by approximately a whole half period. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)

L. Vynnytska et al. / Computers & Geosciences 50 (2013) 95–105

reported in van Keken et al. (1997), which for the ‘‘best’’cases are in the range ð0:01385,0:01462Þ for max urms , and ð49:49,51:12Þ for tðmax urms Þ; in Tackley and King (2003) on the finest mesh, these values are tðmax urms Þ ¼ 51:0 and max urms ¼ 0:01398. As noted, the different parameter choices for kC and C CFL result in some differences in the computed max urms and tðmax urms Þ cf. Table 3. The difference in max urms is small (in the third digit) on the coarsest meshes 40  40; but, the differences increase to some extent as the discretization is refined. The differences are larger for tðmax urms Þ. We conclude that the chemical diffusion

103

coefficient kC plays a role: a lower kC seems to give a higher value for the peak velocity at a later time. Additionally, as we observe from Fig. 5, the chemical diffusivity has an influence on the results in the later stages of evolution. As expected, kC ¼ 106 causes diffusive processes, while for kC ¼ 0 sharp interfaces persist. The computed values for the instability growth rate g show greater variation, and are therefore studied more carefully in Fig. 6. For the case kC ¼0.0, it seems that g converges as the discretization is refined, and that it approaches a value close to 0.105. This is in agreement with the range ð0:1020,0:1051Þ of the

Fig. 5. Rayleigh–Taylor instability: compositional field with kC ¼ 0, C CFL ¼ 0:3 (left column) and kc ¼ 106 , C CFL ¼ 0:5 (right column) at nondimensional times t ¼ 500 (top row), t¼ 1000 (middle row) and t ¼1500 (bottom row); values for composition varies from 0 (blue) to 1 (red); mesh of size 160  160, a ¼ 4:0. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)

104

L. Vynnytska et al. / Computers & Geosciences 50 (2013) 95–105

significantly helped to improve the manuscript. The work of the first and the third authors is supported by a Statoil research grant. The work of the second author is supported by a Center of Excellence grant from the Research Council of Norway to the Center for Biomedical Computing at Simula Research Laboratory.

Appendix A. Supplementary data Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.cageo.2012.05.012.

References

Fig. 6. Rayleigh–Taylor benchmark problem: growth rate of instability versus mesh size n. We observe that the values do not in general behave uniformly with mesh refinement.

‘‘best’’values reported in van Keken et al. (1997). The case kC ¼ 106 , C CFL ¼ 0:5 was examined using one of the methods (labeled SK) by van Keken et al. (1997). That method is comparable to the one used here and this case is therefore also included here for further comparison. The main difference between the two methods lies in the use of a SUPG method rather than a DG method for the composition. For the range of mesh refinement reported there, the results obtained here are clearly comparable. However, g behaves in a non-monotone fashion as the discretization is further refined: the values first increase significantly, and then decrease. A full study of this issue is beyond the scope of this paper. However, we note that the computed g depends on the initial time step, and since an adaptive time-step is used, this value is sampled at different times for the different meshes and the different parameters. Moreover, one could conjecture that the filtering algorithm significantly influences this quantity. 7. Concluding remarks We conclude that the simulation results show agreement with published benchmarks. This paper has demonstrated that the discontinuous Galerkin scheme is an effective alternative for tracking compositional interfaces. The scheme can track such interfaces without significant numerical diffusion or dispersion errors. Overall, we have demonstrated the effectiveness of utilizing the FEniCS Project for the mantle convection simulations. The FEniCS components have indeed enabled easy implementation of the variational problems considered in this paper. Moreover, for benchmarking purposes in general, not only the direct output of the simulations is of interest, but also quantities such as integrals of gradients (12), norms (10) and arbitrary point evaluations. The rich form language provided by FEniCS is helpful in this regard as it allows these functionals to be quickly formulated and computed. Moreover, the use of an embedded finite element library has been particularly beneficial as this allows library functionality to be combined with pure Python (or Cþþ) functionality. In particular, it has allowed easy integration of finite element routines provided by FEniCS with special-purpose code such as for instance the filtering scheme used in Algorithm 3. Acknowledgments The authors are thankful to Dr. Bernhard Steinberger and anonymous reviewer for the constructive comments which

Alnæs, M.S., 2012. UFL: a finite element form language. In: Logg, A., Mardal, K.A., Wells, G. (Eds.), Automated Solution of Differential Equations by the Finite Element Method. Springer, Berlin, Heidelberg, pp. 303–338. Arnold, D.N., Brezzi, F., Cockburn, B., Marini, L.D., 2002. Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM Journal on Numerical Analysis 39, 1749–1779. van den Berg, A., van Keken, P., Yuen, D., 1993. The effects of a composite nonNewtonian and Newtonian rheology on mantle convection. Geophysical Journal International 115, 62–78. Blankenbach, B., Busse, F., Christensen, U., et al., 1989. A benchmark comparison for mantle convection codes. Geophysical Journal International 98, 23–28. Cockburn, B., Karniadakis, G.E., Shu, C.W., 1999. The development of discontinuous Galerkin methods. IMA Preprint Series #1662. /www.ima.umn.edu/preprints/ scanned-preprint/preprints5/1662.pdfS. Cockburn, B., Shu, C.W., 2001. Runge–Kutta discontinuous Galerkin methods for convection-dominated problems. Journal of Scientific Computing 16, 173–261, http://dx.doi.org/10.1023/A:1012873910884. Deschamps, F., Trampert, J., Tackley, P., 2007. Thermochemical structure of the lower mantle: seismological evidence and consequences for geodynamics. In: Yuen, D., Maruyama, S., Karato, S., Windley, B. (Eds.), Superplume: Beyond Plate Tectonics. Springer, pp. 293–320. Hake, J., Logg, A., Wells, G.N., 2011. DOLFIN: A Cþþ/Python Finite Element Library. Hansen, U., Ebel, A., 1988. Time-dependent thermal convection—a possible explanation for a multiscale flow in the Earth’s mantle. Geophysical Journal 94, 181–191. Hansen, U., Yuen, D., 1988. Numerical simulation of thermal–chemical instabilities at the core–mantle boundary. Nature 334, 237–240. van Keken, P., King, S., Schmeling, H., et al., 1997. A comparison of methods for the modeling of thermochemical convection. Journal of Geophysical Research 102, 22477–22495. King, S., Lee, C., van Keken, P., et al., 2010. A community benchmark for 2D Cartesian compressible convection for the Earth’s mantle. Geophysical Journal International 180, 73–87. Kirby, R.C., 2004. Algorithm 839: FIAT, a new paradigm for computing finite element basis functions. ACM Transactions on Mathematical Software 30, 502–516. Kirby, R.C., Logg, A., 2006. A compiler for variational forms. ACM Transactions on Mathematical Software 32, 417–444. Kirby, R.C., Logg, A., 2007. Efficient compilation of a class of variational forms. ACM Transactions on Mathematical Software 33. 17:1–17:20. Lenardic, A., Kaula, W., 1993. A numerical treatment of geodynamic viscous flow problems involving the advection of material interfaces. Journal of Geophysical Research 98, 8243–8260. Lennie, T., McKenzie, D., Moore, D., Weiss, N., 1988. The breakdown of steady convection. Journal of Fluid Mechanics 188, 47–85. Li, B.Q., 2006. Discontinuous finite elements in fluid dynamics and heat transfer. In: Computational Fluid and Solid Mechanics. Springer-Verlag London Limited. 595 pp. Logg, A., Mardal, K.A., Wells, G.N. (Eds.), 2012. Automated Solution of Differential Equation by the Finite Element Method. Springer, Berlin, Heidelberg. 723 pp. Logg, A., Wells, G.N., 2010. DOLFIN: automated finite element computing. ACM Transactions on Mathematical Software 32, 1–28. ¨ Moresi, L.N., Dufour, F., Muhlhaus, H.B., 2003. A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials. Journal of Computational Physics 184, 476–497. Ølgaard, K.B., Logg, A., Wells, G.N., 2008. Automated code generation for discontinuous Galerkin methods. SIAM Journal on Scientific Computing 31, 849–864. Ølgaard, K.B., Wells, G.N., 2010. Optimizations for quadrature representations of finite element tensors through automated code generation. ACM Transactions on Mathematical Software 37. 8:1–8:23. di Pietro, D., Forte, S.L., Parolini, N., 2006. Mass preserving finite element implementations of the level set method. Applied Numerical Mathematics 56, 1179–1195. Schubert, G., Turcotte, D., Olson, P., 2001. Mantle Convection in the Earth and Planets. Cambridge University Press 956 pp.

L. Vynnytska et al. / Computers & Geosciences 50 (2013) 95–105

Steinberger, B., Holme, R., 2008. Mantle flow models with core–mantle boundary constraints and chemical heterogeneities in the lowermost mantle. Journal of Geophysical Research 113. Tackley, P.J., King, S.D., 2003. Testing the tracer ratio method for modeling active compositional fields in mantle convection simulations. Geochemistry, Geophysics, Geosystems 4. Tan, E., Choi, E., Thoutireddy, P., Gurnis, M., Aivazis, M., 2006. GeoFramework: coupling multiple models of mantle convection within a computational framework. Geochemistry, Geophysics, Geosystems 7. Taylor, C., Hood, P., 1973. A numerical solution of the Navier–Stokes equations using the finite element technique. Computers & Fluids 1, 73–100.

105

The FEniCS Team, 2011. The FEniCS Project /http://www.fenicsproject.orgS. Zhong, S., Zuber, M., Moresi, L., Gurnis, M., 2000. The role of temperaturedependent viscosity and surface plates in spherical shell models of mantle convection. Journal of Geophysical Research 105, 11063–11082. Zienkiewicz, O., Taylor, R., 2000. The Finite Element Method, Fluid Dynamics, vol. 3, 5th ed. Butterworth-Heinemann, 352 pp. Zienkiewicz, O.C., Taylor, R.L., Sherwin, S.J., Peira, J., 2003. On discontinuous Galerkin methods. International Journal for Numerical Methods in Engineering 58, 1119–1148.