www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE Operator splitting and discontinuous Galerkin methods for advection-reaction-diffusion problem. Application to plant root growth Emilie Peynaud CIRAD, UMR AMAP, Yaoundé, Cameroun AMAP, University of Montpellier, CIRAD, CNRS, INRA, IRD, Montpellier, France University of Yaoundé 1, National Advanced School of Engineering, Yaoundé, Cameroon emilie.peynaud@cirad.fr Received: 19 September 2017, accepted: 3 December 2018, published: 18 December 2018 Abstract—Motivated by the need of developing numerical tools for the simulation of plant root growth, this article deals with the numerical res- olution of the C-Root model. This model describes the dynamics of plant root apices in the soil and it consists in a time dependent advection-reaction- diffusion equation whose unique unknown is the density of apices. The work is focused on the implementation and validation of a suitable nu- merical method for the resolution of the C-Root model on unstructured meshes. The model is solved using Discontinuous Galerkin (DG) finite elements combined with an operator splitting technique. After a brief presentation of the numerical method, the implementation of the algorithm is validated in a simple test case, for which an analytic expression of the solution is known. Then, the issue of the positivity preservation is discussed. Finally, the DG- splitting algorithm is applied to a more realistic root system and the results are discussed. Keywords-Time dependent advection-reaction- diffusion; Operator splitting; Discontinuous Galerkin method; Plant root growth simulation; I. INTRODUCTION The article is devoted to the numerical mod- eling of plant root growth. This work has been originally motivated by the need of developing numerical tools for the simulation of plant growth dynamics. Due to the difficulty of doing non- destructive observations of the underground part of plants (that allow to do long term studies of the dynamics of tree roots for example), mathematical models are achieving an essential role. Several theoretical and numerical challenges arise in the field of the simulation of the dynamics of plant roots [48], [47], [38], [2], [39]. The mathemat- ical description of plant root is not trivial, due to the presence of many interactions arising in the rhizosphere and also due to the diversity of plant root types. Mathematical models based on the use of partial differential equations are useful tools to simulate the evolution of root densities in Copyright: c© 2018 Peynaud. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Citation: Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction-diffusion problem. Application to plant root growth, Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 1 of 19 http://www.biomathforum.org/biomath/index.php/biomath https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... space and time [43], [44], [44], [45], [46], [41], [40], [1]. This formalism facilitates the coupling with physical models such as water and nutrient transports [42], [43], [44], [41], [49]. And the com- putational time for the simulation of such models is not dependent on the number of roots which is useful for applications at large scale. The C-Root model [1] is a generic model of the dynamics of root density growth. This model takes only one unknown which is related to root densities such as the density of apices, root length density or biomass density. It has only three parameters. The model is said to be generic in the sense that it can apply to a wide variety of root system types. The model consists in a single time-dependent advection-reaction-diffusion equation, and one of the challenge is to numerically solve the equation. In [1] and [2] the authors solved the problem with the finite difference method on Cartesian mesh grids combined with an operator splitting technique. Unfortunately, Cartesian mesh grids do not allow easily to mesh complex soil geometries. From the theoretical and computational point of view, Cartesian grids also lead to difficulties for a rigorous study and validation of the model. That is why this article focus on the development and implementation of a suitable numerical method for the resolution of the C-Root model on tri- angular mesh grids, that allow to mesh complex geometries. However, one of the main difficulties in the C-root model is that the advection and diffusion terms are not always of the same order of magnitude. It depends on the phase of the root system development [2]. As a result, the properties of the equation may vary along the simulation: it can be either close to a hyperbolic problem or close to a parabolic problem. In a previous work [3], the use of the Discon- tinuous Galerkin method has been implemented and validated. Indeed, the usual choice of the classical Lagrange finite element method suffers from a lack of stability when the advection term is dominant [4]. For this reason, we implemented a discontinuous Galerkin (DG) method for both the advection and diffusion terms. All the three operators where solved simultaneously using the same time approximation scheme (θ-scheme). However, as explained in [6], for multi- biophysic problems it is not efficient to use the same numerical scheme for the different operators of the system. For example, we may want to use the Euler explicit scheme for the advection term and an Euler implicit scheme for the diffusion. The operator splitting technique [7], [8] is a well known alternative for the resolution of equations having a multi-biophysic behaviour that allows the use of different time schemes for each operator of the equation. The idea of the splitting technique is to split the problem into smaller and simpler parts of the problem so that each part can be solved by an efficient and suitable time scheme. This methods has been used for a wide range of applica- tions dealing with the advection-reaction-diffusion equation [9]. Operator splitting techniques have been extensively used in combination with finite difference methods [10], [2], finite volume meth- ods [11], [12] but also with Continuous Galerkin methods [13], [14], [15], [16], [17]. To the best of my knowledge, only very few articles deal with the use of the operator splitting technique in combination with the discontinuous Galerkin approximation [18], [19], [20], [21]. In this paper, we present a new application of the operator splitting technique combined with discontinuous finite elements. The paper is structured as follows. In section II, the C-root growth model [1], [2], [3] is briefly described. An analysis is also provided, where I showed the existence and uniqueness of a positive real solution. In section III, the splitting operator technique is introduced and applied to the C-Root model, combined with the use of discontinuous Galerkin approximations. In section IV, the algo- ritm is implemented and validated using a simple test case for which an analytic expression of the solution is known. As an application, I provide simulations of the development of eucalyptus roots in section V. Finally, the paper ends with a con- clusion and further improvements. Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 2 of 19 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... II. THE MODEL A. Modelling root growth with PDE: the C-Root model The C-Root model [1] was developed to sim- ulate the growth of dense root networks, usually composed of fine roots, with negligible secondary thickening. As presented in [1], the unknown vari- able u is the number of apices per unit volume, but it can also stand for the density of fine root biomass. The soil is considered as a subdomain of Rd (with d = 1, 2 or 3). It is assumed that Ω has smooth boundaries (Lipschitz boundaries) de- noted ∂Ω. The C-Root model combines advection, diffusion and reaction, which aggregate the main biological processes involved in root growth, such as primary growth, ramification and root death. The reaction operator gives the quantity of apices (or root biomass) produced in time, whereas ad- vection and diffusion operators spatially distribute the whole apices (or biomass) in the domain. The reaction operator describes the evolution in time of the root biomass in a given domain. In the C-Root model it is a linear term characterized by the scalar parameter ρ which is the growth rate of the root system. The diffusion corresponds to the spread of the root biomass over space. It is described by the parameter σ which is a d × d matrix that characterizes the growth of the root biomass in any direction exploiting free space in the soil. The advection corresponds to the displacement of the root biomass in a direction and velocity given by v which is a vector in Rd. On the boundaries of Ω, what happens for the quantity being transported is different depending if the growth makes the roots to come inside Ω or to go outside of Ω. If v is going inside Ω (at the inlet boundary) the root biomass u will enter the domain and increase. On edges where v is going out of the domain (outlet boundary) the root biomass u is going to be pushed out of Ω. Since this phenomena is oriented (causality) and the behaviour of the solution is different on inlet and outlet boundaries, we need to specify in the model these parts of the boundaries. Mathematically, it is required to define the inlet boundary with respect to v as ∂Ω− = {x ∈ ∂Ω : (v ·n)(x) < 0} . (1) The outlet boundary Ω+ is given by ∂Ω+ = ∂Ω\∂Ω−. The dynamics of the root system is stud- ied between the time t0 and t1 with 0 ≤ t0 < t1. The problem reads as follow: find u such that  ∂tu + v ·∇u−∇· (σ∇u) + ρu = 0 in ]t0, t1[×Ω u(t0) = u0 at {t0}× Ω n ·σ∇u = g on ]t0, t1[×∂Ω (n ·v)u = gin on ]t0, t1[×∂Ω− (2) where g ∈ L2(∂Ω) and gin ∈ L2(∂Ω−) are given. And u0 is the given initial solution. Problem (2) is known as the time dependent advection- reaction-diffusion problem and belongs to the class of parabolic partial differential equa- tions. This equation is a model problem that often occurs in fluid mechanics but also in many other applications in life sciences (see for instance [22], [23], [24]). Depending on the boundary conditions, the problem has different meanings. To simplify the presentation we only consider the Neumann boundary condition combined with an inlet bound- ary condition at the inlet of the domain. The Neu- mann condition specifies the value of the normal derivative of the solution at the boundary of the domain. The inlet boundary condition specifies the quantity of u convected by v that enters in the domain. B. The weak problem Since the goal is to solve the problem on unstructured meshes, the spatial operators are ap- proximated using finite element methods. Within this framework, it is classical to write the problem in a variational form. Let us first introduce some functional spaces [50]. • The space H1(Ω) defined such that H1(Ω) = {v ∈ L2(Ω) : ∇v ∈ L2(Ω)} is a Hilbert space when equipped with the norm ‖ · ‖1,Ω. We recall that ∀v ∈ H1(Ω), ‖v‖1,Ω = (v,v)1,Ω Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 3 of 19 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... and the scalar product (·, ·)1,Ω is defined by ∀v ∈ H1(Ω), (u,v)1,Ω = ∫ Ω uv dx + ∫ Ω ∇u ·∇v dx. • We denote L2(]t0, t1[,H) the space of H- valued functions whose norm in H is in L2(]t0, t1[). This space is a Hilbert space for the norm ‖u‖L2(]t0,t1[,H) = (∫ t1 t0 ‖u(t)‖2H )1/2 . • Let B0 ⊂ B1 be two reflexive Hilbert spaces with continuous embedding, we de- note W(B0,B1) the space of functions v : ]t0, t1[−→ B0 such that v ∈ L2(]t0, t1[,B0) and dtv ∈ L2(]t0, t1[,B1). Equipped with the norm ‖u‖W(B0,B1) = ‖u‖L2(]t0,t1[,B0) +‖dtu‖L2(]t0,t1[,B1), the space W(B0,B1) is a Hilbert space [25]. Using the previous functional spaces, I now define the problem in the following weak form: Find u in W such that ∀v ∈ H 〈dtu(t),v(t)〉H′,H +a(t,u,v) =`(t,v) a.e. t∈]t0,t1[ u(t0) = u0, (3) where W = W(H1(Ω), (H1(Ω))′) and H = H1(Ω) and `(t,v) = ∫ ∂Ω g(t)vdγ (4) a(t,u,v) =aA(t,u,v)+aD(t,u,v)+aR(t,u,v) (5) with aA(t,u,v) = ∫ Ω v(v(t,x) ·∇u) dx, (6) aD(t,u,v) = ∫ Ω ∇v ·σ(t,x) ·∇udx, (7) aR(t,u,v) = ∫ Ω ρ(t)uv dx. (8) One can prove that problems (3) and (2) are equivalent almost everywhere in ]t0, t1[×Ω. Let us assume that there is a constant σ0 > 0 such that ∀ξ ∈ Rd, d∑ i,j=1 σijξiξj ≥ σ0‖ξ‖2d a.e. in Ω. (9) In addition, I assume that inf x∈Ω ( σ − 1 2 (∇·v) ) > 0 and inf x∈∂Ω (v ·n) ≥ 0. (10) Under assumption (9) and (10), one can prove that the problem is well-posed for sufficiently smooth v, σ and ρ (see for instance [25]). C. The positivity preserving property of the solu- tion In the framework of our applications to the simulation of root biomass densities one of the crucial property of the problem is the preservation of the positity of the solution along time. For a positive initial solution u0, the solution of (3) stays positive. Proposition II-C.1. Let u0 ∈ L2(Ω) and f ∈ L2(]t0, t1[,L 2(Ω)). We consider u the solution of (3) in W . We assume that u0(x) ≥ 0 a.e. in Ω and g(t,x) ≥ 0 a.e. in ]t0, t1[×∂Ω. Then u(t,x) ≥ 0 a.e in ]t0, t1[×Ω. Proof: I follow [25]. See also [26], [27]. We consider the function u− defined by u− = 1 2 (|u|−u). Let us note that u−= { 0 a.e in ]t0,t1[×Ω, if u≥0 a.e in ]t0,t1[×Ω, −u a.e in ]t0,t1[×Ω, if u<0 a.e in ]t0,t1[×Ω. That is we have u− ≥ 0 a.e. in ]t0, t1[×Ω. (11) We verify that u− is an admissible test function in W . Using the following obvious equations (∇|u|)2 = (∇u)2 ∇|u| ·∇|u| = ∇u ·∇u u∇|u| = |u|∇u u∇u = |u|∇|u| Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 4 of 19 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... that are valid a.e in ]t0, t1[×Ω we can verify that a(t,u−,u−) = −a(t,u,u−). By adding the same quantity on both sides of the equation we get 〈dtu−,u−〉+a(t,u−,u−) =〈dtu−,u−〉−a(t,u,u−). Since u satisfy (3) we have 〈dtu−,u−〉+a(t,u−,u−) =〈dtu−,u−〉+〈dtu,u−〉 − `(t,u−). One can notice that 〈dtu−,u−〉 + 〈dtu,u−〉 = 0. Then we have 1 2 dt‖u−‖20,Ω + a(t,u −,u−) = −`(t,u−) ≤ 0, with g(t,x) ≥ 0 a.e. in ]t0, t1[×∂Ω. Now from the coercivity of the bilinear form a we obtain 1 2 dt‖u−‖20,Ω + c‖u −‖20,Ω ≤ 1 2 dt‖u−‖20,Ω + a(t; u −,u−) ≤ 0, where c is a strictly positive constant. The estimate is then 1 2 dt‖u−‖20,Ω ≤−c‖u −‖20,Ω. By the Gronwall lemma we have that ∀t ∈ [t0, t1] × Ω,‖u−(t)‖20,Ω ≤ e −2ct‖u−(0)‖20,Ω. Since c > 0 and t ≥ t0 ≥ 0 , we have that e−2ct ≤ 1 , so we obtain ∀t ∈ [t0, t1] × Ω,‖u−(t)‖20,Ω ≤‖u −(0)‖20,Ω. Since u(0) = u0 ≥ 0 a.e in Ω we have u−(0) = 0 a.e in Ω. So we deduce that ∀t ∈ [t0, t1] × Ω,‖u−(t)‖20,Ω ≤ 0. But from the definition of u− we have u− ≥ 0 a.e in ]t0, t1[×Ω. So we deduce that ‖u−(t)‖20,Ω = 0 and thus u−(t) = 0 a.e in ]t0, t1[×Ω. It means that u ≥ 0 a.e in ]t0, t1[×Ω by definition of u−. III. APPROXIMATION OF THE MODEL A. The operator splitting technique Here we focus on the implementation of the operator splitting technique. The time inter- val [t0, t1] is divided in N subspaces of size δt such that [t0, t1] = ∪n=1,N ]tn, tn+1[ with ∩n=1,N ]tn, tn+1[= ∅. At each iteration step we solve the following problems • Find uA ∈ H such that ∀v ∈ H, for a.e t ∈ ]tn, tn+1[, 〈dtuA(t),v(t)〉H′,H + aA(t,uA,v) = 0 uA(tn) = u(tn). • Find uD ∈ H such that ∀v ∈ H, for a.e t ∈]tn, tn+1[, 〈dtuD(t),v(t)〉H′,H + aD(t,uD,v) = `(t,v) uD(tn) = uA(tn+1). • Find uR ∈ H such that ∀v ∈ H, for a.e t ∈ ]tn, tn+1[, 〈dtuR(t),v(t)〉H′,H + aR(t,uR,v) = 0 uR(tn) = uD(tn+1). • Set u(tn+1) = uR(tn+1). The bilinear forms aA(t,u,v), aD(t,u,v) and aR(t,u,v) are respectively given by (6), (7) and (8). And `(t,v) is the linear form (4). If the operators are commutative, then the splitting error vanishes. Otherwise, if the operators are not com- mutative, then the splitting error does not vanish and a second order splitting would be required (see [6]). In the following, I present the different schemes related to each operator. B. The advection step: DG upwind scheme The advection step consists in solving the fol- lowing transport problem : Find u such that ∀v ∈ H, for a.e t ∈]tn, tn+1[, 〈dtu(t),v(t)〉H′,H + aA(t,u,v) = 0 (12) uA(tn) = uR(tn) (13) where aA(t,u,v) is the bilinear form (6). For the space approximation of this problem, we imple- mented the DG upwind method presented below. Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 5 of 19 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... Let Th be a regular family of decomposition in triangles of the domain Ω such that Ω = N⋃ i=1 K̄i and Ki ∩Kj = ∅,∀i 6= j. The h subscript in Th denotes the size of the mesh cells and it is defined by h = max K∈Th hK where hK is the diameter of the element K. Let Eh be the set of edges of the elements of Th. Among the elements of Eh we denote by Ebh the set of edges belonging to ∂Ω. The sets Eb,−h and E b,+ h are the sets of edges belonging to ∂Ω− and ∂Ω+ respectively. And Eih is the set of interior edges. Let us consider an element of Eih. We denote by T + and T− the two mesh elements sharing the edge e so that e = ∂T + ∩∂T− where the minus and plus superscripts depend on the direction of the advection vector. By convention we suppose that v goes from T− to T + that is v ·n+e < 0 and v · n−e > 0 where n+e (resp. n−e ) is the outward normal vector of e in T + (resp. T−). When it is not necessary to distinguish the orientation of the normal vectors n+e and n − e we denote by n the unitary normal of e. Let us consider the advection problem on each element Ki of the domain : for all Ki, i = 1,N we look for u the solution of the equation (12) defined on Ki. Similarly to the problem defined on all the domain Ω, we look for a solution u that is in L2(Ki) and such that ∇u is in L2(Ki) for all Ki in Th. Let us introduce the following broken Sobolev space: H1(Th) = { v ∈ L2(Ω) : ∇v ∈ L2(Ki) and v ∈ H1/2+ε(Ki),∀Ki ∈Th } with ε a positive real number. The trace of the functions of H1(Th) are meaningful on e ⊂ Ki, ∀Ki ∈ Th. The functions v of H1(Th) have two traces along the edges e. We denote v+e the trace of v along e on the side of triangle T + and v−e the trace of v along e on the side of T−. On edges that are subsets of ∂Ω the trace is unique and we can note v+e = v if e ∈E b,− h and v − e = v if e ∈E b,+ h , and by convention, we set v−e = 0 if e ∈E b,− h and v + e = 0 if e ∈E b,+ h . The jump of functions of H1(Th) across the inter- nal edge e is defined by: JvK = v+e −v − e ,∀e ∈E i h. For edges belonging to the boundary of Ω we take JvK = ve,∀e ∈E b,− h and JvK = −ve,∀e ∈E b,+ h , with ve the trace of v along e. The mean value of u on e is defined by {{v}} = 1 2 (v+e + v − e ),∀e ∈E i h. Besides for edges on the boundaries we take {{v}} = ve,∀e ∈Ebh. Let us denote by X the functional space defined such that X = {v : ]t0, t1[−→ H1(Th) : v ∈ L2(]t0, t1[,H1(Th)); and dtv ∈ L2(]t0, t1[,H1(Th)′)}. This space is a Hilbert space equipped with the norm ‖v‖X = ‖v‖L2(]t0,t1[,H1(Th)) +‖v‖L2(]t0,t1[,H1(Th)′). The DG variational formulation of the advection step written on the broken Sobolev space takes the following form: Find u in X such that for a.e t ∈]t0, t1[, ∀v ∈ H1(Th) 〈dtu(t),v(t)〉H1(Th)′,H1(Th) +a up h (t;u,v) =` up h (t;v), u(t0) = u0, where the form auph (t; u,v) is the approximation of the advection term. It consists in the upwind formulation of the DG method [28]. It reads: a up h (t; u,v) = ∑ K∈Th ∫ K u(ρv−v ·∇v)dx − ∑ e∈Eb,±h ∪E b,+ h ∫ e |v ·n+e |u − e JvKds. (14) Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 6 of 19 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... The approximated linear form of the the right hand side reads ` up h (t; v) = − ∑ e∈Eb,−h ∫ e (v ·n+e )ginv + e ds. The DG-formulation (14) is consistent and stable, see for example [32]. The discontinuous Galerkin method consists in searching the solution in the approximation space Xh defined such that Xh = { v :]t0, t1[−→ Wkh ; v ∈ L 2(]t0, t1[,W k h ); and dtv ∈ L2(]t0, t1[, (Wkh ) ′) } , where Wkh is given by Wkh = { vh ∈ L2(Ω);∀K ∈Th,vh|K ∈ Pk } . Let us note that the functions of Wkh can be discontinuous from one element of the mesh to the other. Let us note that Wkh is embedded in H 1(Th) so that Xh ⊂ X . This problem can be written in a matrix form. Let us denote (λi)i=1,n the basis of the finite dimensional subspace Wkh where n = dim(Wkh ). In this basis the approximated solution takes the form: uh(t,x,y) = n∑ i=1 ξi(t)λi(x,y), where the ξi(t) are the degrees of freedom. Let us define X the vector of degrees of freedom: X(t) = (ξ1(t), . . . ,ξn(t)) T . The approximated problem then reduces to find X(t) ∈ [C2(0,T)]n such that M dX(t) dt + Aup(t)X(t) = L up h (t) MX(0) = MX0 where M and Aup(t) are two matrices defined such that M = (Mi,j)i,j and Mi,j = ∑ T∈Th ∫ K λiλjdx, (15) Aup = ( A up i,j ) i,j and Aupi,j = a up h (t; ,u,v), (16) and Luph (t) is the vector of size n defined such that( L up h (t) ) i = ` up h (t; λi) for i = 1,n. The problem reduces to a linear system of ordinary differential equations. The time approximation is based on a finite difference scheme. At each iteration step we solve the following problem: Find XN+1 ∈ Rn such that 1 δt M ( XN+1 −XN ) + (1 −θ)AupXN + θAupXN+1 (17) = (1 −θ)Lup,Nh + θL up,N+1 h and MX0 = MX0, where θ is a real parameter taken in [0, 1]. For θ = 0, we have the explicit Euler schema. For θ = 1, it is the implicit Euler schema. For θ = 1/2, it is the Crank-Nicolson schema. C. The diffusion step The diffusion step consists in solving the fol- lowing problem : Find u such that ∀v ∈ H, for a.e t ∈]tn, tn+1[, 〈dtu(t),v(t)〉H′,H + aD(t; u,v) = `(t; v) u(tn) = uA(tn) where aD(t; u,v) is the bilinear form (7) and `(t; v) is the linear form (4). In the setting intro- duced before, the DG variational formulation of the diffusion step written in the broken Sobolev space takes the following form: Find u in X such that ∀v ∈ H1(Th), for a.e. t ∈]t0, t1[ 〈dtu(t),v〉H1(Th)′,H1(Th) + a ip h (t; u,v) = ` ip h (t; v) u(t0) = u0. The form aiph (t; u,v) is the approximation of the diffusion term. It consists in the interior penalty Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 7 of 19 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... formulation (IP) that reads a ip h (t; u,v) = ∑ K∈Th ∫ K σ∇u ·∇v dx − ∑ e∈Eih ∫ e {{σ∇u}} ·n+e JvK ds + ∑ e∈Eih ∫ e {{σ∇v}} ·n+e JuK ds + ∑ e∈Eih η he ∫ e JuKJvK ds, where η is a positive penalization factor. The linear form `iph (t; v) is given by ` ip h (t; v) = ∑ e∈Ebh ∫ e gv ds. This formulation was introduced in [31] and is known as the non-symmetric interior penalty (NSIP) formulation, see [30], [32]. In matrix form the problem reduces to find X(t) ∈ [C2(0,T)]n such that M dX(t) dt + Aip(t)X(t) = L ip h (t) MX(0) = MX0 where M is defined by (15) and Aip is defined such that Aip = ( A ip i,j ) i,j and Aupi,j = a ip h (t; ,u,v). The vector Liph (t) is such that ( L ip h (t) ) i = ` ip h (t; λi) for i = 1,n. Similarly to the advection step, the time approximation of the problem is based on a finite difference scheme of the form (17). D. The reaction step The reaction step consists in solving the follow- ing problem : Find u such that ∀v ∈ H, for a.e. t ∈]tn, tn+1[ 〈dtu(t),v(t)〉H′,H + aR(t; u,v) = 0 u(tn) = uD(tn) where aR(t; u,v) is the bilinear form (8). This problem takes the following matrix form find X(t) ∈ [C2(0,T)]n such that dX(t) dt + ρX(t) = 0 X(0) = X0 where we recall that ρ is a constant real parameter. This problem can be solved by an exact scheme (a kind of schemes that provide exact solutions, i.e. a solution equal to the analytical solution). At each iteration we find XN+1 such that 1 Φ(δt) ( XN+1 −XN ) = −ρXN with Φ(δt) = 1 ρ (1 − exp(−ρδt)). This scheme is unconditionally stable, meaning that we can choose the time step independently from the space step. It is also positively stable, meaning that if XN ≥ 0 so is XN+1. IV. VALIDATION OF THE SPLITTING ALGORITHM WITH A SIMPLE TEST CASE Problem (3) has been already solved using dis- continuous Galerkin elements (DG) [3]. Advection and Diffusion operators were solved simultane- ously using the Crank-Nicolson scheme providing stable results. However, even for simple test cases some simulations did not always provide positive numerical solutions. One reason is that the same time approximation scheme is not necessarily suit- able for both the advection and for the diffusion. That is why a new operator splitting algorithm has been implemented with a different time scheme for each operator. The goal of this section is to validate the im- plementation of the code. To this end I compare the convergence of the approximation with and without the splitting technique. I briefly explore the question of the positivity of the approximated solution. A. Description of the simple test-case First let me introduce a simplified test-case for the validation of the splitting algorithm. Set L > 0, and Ω =] − L; L[2. Let v = (v1,v2) ∈ R2 and Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 8 of 19 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... d ∈ R be a constant and 0 ≤ t0 ≤ t1. Find u such that ∂u ∂t + v ·∇u + ρu = d∆u in ]t0, t1[×Ω, u(x,y, 0) = u0(x,y) on {t0}× Ω, (18) n ·∇u = g on ]t0, t1[×∂Ω. n ·vu = gin on ]t0, t1[×∂Ω−. The initial condition and the boundary condition are chosen such that the solution of problem (18) is explicitly given by ∀(x,y,t) ∈ Ω×]t0, t1[ u(x,y,t) = c0 ( a2 a2 + td ) κ(x,y,t)e−ρt. with κ(x,y,t) = c0 exp ( − (x−x0 − tv1)2 + (y −y0 − tv2)2 4(a2 + td) ) where c0 > 0, a > 0, x0 and y0 are real parameters and v1 and v2 are the two components of v. Notice that u(x,y,t) > 0 for all (x,y,t) in Ω×]t0, t1[. B. Numerical validation and convergence To validate the implementation of the splitting technique, I ran the previous test case with differ- ent mesh sizes and time steps and I computed the global L2-errors such that eh = ( δt N∑ k=1 ‖u(tk) −uh(tk)‖20,Ω )1/2 where tk = t0 + kδt, with k ∈ N+∗ and tN = t1. The flexibility of the splitting technique allows to choose different time schemes for each operator. I consider a θ-scheme with θ = 0 (explicit Euler), θ = 1 (implicit Euler), and θ = 1 2 (Crank- Nicolson) for both the advection step and the diffusion step, and I consider an exact scheme for the reaction step. For the simulations I took the parameters such that v = (0.1, 0)T , σ = 0.01 and ρ = −1. The triangular meshes used for the simu- lations are identified by h which is the size of the biggest triangle of the mesh. Table I, page 9, gives the number of triangles and the number of nodes of each mesh used for the simulations. Choosing h (≈) number of triangles number of nodes 2.63 × 10−1 68 45 1.31 × 10−1 272 157 6.57 × 10−2 1 088 585 3.29 × 10−2 4 352 2 257 1.64 × 10−2 17 408 8 865 8.22 × 10−3 69 632 35 137 4.11 × 10−3 278 528 139 905 TABLE I TRIANGULAR MESHES USED FOR THE SIMULATIONS. Fig. 1. Solution of the validation test case at t = t0 (left) and t = t1 (right) computed using the DG method with p1- finite elements and the Euler implicit scheme (θ = 1) and the operator splitting technique with h ≈ 8.2 × 10−3 and δt = 10−2. L = 1/2, the simulations are performed between t0 = 0 and t1 = 1 for different values of the time step δt. Fig. 1, page 9, shows the solution at t = t0 and t = t1. The code is implemented in Fortran 90 and it is run under a 64-bit Linux operating system on a 8-core processor Intel R©CoreTMi7- 7820HQ at a frequency of 2.9GHz and with 32 GB of RAM. The sparse matrices resulting from the finite element approximation are inverted using a solver provided by the library MUMPS [51], [52]. According to Fig. 2, page 10, all the three tem- poral schemes provide results with approximately the same level of accuracy with a spatial con- vergence rate of 2 computed with the global L2- Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 9 of 19 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... Fig. 2. Convergence of the solution with respect to the mesh size: plot of the total L2-error computed between t = 0 and t = 1 with and without the splitting technique for the explicit Euler scheme (θ = 0), the implicit Euler scheme (θ = 1) and the Crank-Nicolson scheme (θ = 1/2) for δt = 5 × 10−5. Fig. 3. Convergence of the solution with respect to the time step: plot of the total L2-error computed between t = 0 and t = 1 with and without the splitting technique for the implicit Euler scheme (θ = 1) and the Crank-Nicolson scheme (θ = 1/2) for h = 4.1 × 10−3. norm. The same order of convergence is obtained when the problem is solved without the splitting technique. Figure 3 on page 10 shows that the Crank- Nicolson scheme (θ = 1/2) converges in δt2 while the Euler Implicit scheme (θ = 1) converges in δt, with and without the splitting technique. The convergence rate in time has to be computed with a really refined mesh grid (here h ≈ 4.1 × 10−3). Fig. 4. Validation of the test case: plot of the CPU time against the mesh size (h) for the computations performed with a processor Intel R©CoreTMi7-7820HQ at 2.9 GHz and RAM 32 GB, between t = 0 and t = 1 with and without the splitting technique for the explicit Euler scheme (θ = 0), the implicit Euler scheme (θ = 1) and the Crank-Nicolson scheme (θ = 1/2) for δt = 5 × 10−5. Fig. 5. Validation of the test case: plot of the CPU time against the time step (δt) for the computations performed with a processor Intel R©CoreTMi7-7820HQ at 2.9 GHz and RAM 32 GB, between t = 0 and t = 1 with and without the splitting technique for the explicit Euler scheme (θ = 0), the implicit Euler scheme (θ = 1) and the Crank-Nicolson scheme (θ = 1/2) for h ≈ 4.1 × 10−3 and δt ranging from 5 × 10−1 to 5 × 10−4. Note that the computations performed here with θ = 0 gave unstable results. Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 10 of 19 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... δt global L2-errors min dof t+ CPU time 1 · 10−3 unstable unstable - 95 s 2 · 10−4 1.24 · 10−4 −1 · 10−4 0.221 475 s 1 · 10−4 1.22 · 10−4 −1 · 10−4 0.218 838 s 5 · 10−5 1.21 · 10−4 −1 · 10−4 0.216 1672 s 2.5 · 10−5 1.20 · 10−4 −9 · 10−5 0.215 3376 s 1 · 10−5 1.20 · 10−4 −9 · 10−5 0.215 10183 s TABLE II COMPUTATIONS PERFORMED WITH THE SPLITTING TECHNIQUE AND THE EXPLICIT EULER SCHEME (θ = 0) WITH h ≈ 1.6 · 10−2 (INTEL R©CORETM I7-7820HQ AT 2.9 GHZ, RAM 32 GB). It results an additional cost in term of CPU time, since it behaves like 1/h2, as shown on figure 4 page 10. For bigger values of h the plot of the errors gave convergence order in time less than 1 and 2 for the implicit Euler scheme and the Crank-Nicolson scheme respectively. As expected, the explicit Euler scheme is conditionally stable, such that, when the CFL condition is fulfilled, the computational time becomes prohibitive. Indeed, it behaves like 1/δt, as shown on figure 5. For instance, the computation with h ≈ 4.1 × 10−3 and δt = 10−5 takes more than 30 hours with the device specified above. That is why, in the rest of the paper, we will only focus on implicit Euler and Crank-Nicolson schemes. However I present here additional computations performed with a bigger mesh size (h ≈ 1.6×10−2) and smaller time steps chosen such that the CFL condition is fulfilled. The global L2-errors and the CPU time are shown on table II. Clearly, the mesh is not fine enough to recover the convergence order in δt, indeed decreasing the time step results only in an increase of the computational time but not in a significant decrease of the errors. C. Some comments on the positivity 1) Positivity of the full problem: Table III on page 11 and table V on page 12 give the min- imum values of the degrees of freedom (dof) obtained during the simulations performed respec- tively with and without the splitting technique. The minimum value of the dof is defined such that mintk (mini=1,n X k i ) where X k i is the i th dof at time tk. This quantity gives an idea about the δt h≈1.6·10−2 h≈8.2·10−3 h≈4.1·10−3 5 · 10−1 −7 · 10−1 −7 · 10−1 −7 · 10−1 2 · 10−1 −2 · 10−1 −2 · 10−1 −2 · 10−1 1 · 10−1 −8 · 10−4 −1 · 10−3 −2 · 10−3 4 · 10−2 −2 · 10−4 −6 · 10−4 −1 · 10−3 2 · 10−2 −3 · 10−13 −2 · 10−4 −6 · 10−4 1 · 10−2 −6 · 10−5 −2 · 10−20 −1 · 10−4 4 · 10−3 −3 · 10−4 4 · 10−65 5 · 10−66 2 · 10−3 −3 · 10−4 −9 · 10−11 5 · 10−88 1 · 10−3 −2 · 10−4 −8 · 10−10 8 · 10−114 1 · 10−4 −9 · 10−5 −1 · 10−11 −4 · 10−35 Crank-Nicolson scheme (θ = 1/2) δt h≈1.6·10−2 h≈8.2·10−3 h≈4.1·10−3 5 · 10−1 −5 · 10−8 −9 · 10−9 −2 · 10−9 2 · 10−1 4 · 10−9 1 · 10−9 4 · 10−10 1 · 10−1 3 · 10−11 3 · 10−11 1 · 10−11 4 · 10−2 6 · 10−17 5 · 10−17 5 · 10−17 2 · 10−2 3 · 10−23 2 · 10−23 2 · 10−23 1 · 10−2 1 · 10−31 4 · 10−32 3 · 10−32 4 · 10−3 −3 · 10−5 1 · 10−48 6 · 10−49 2 · 10−3 −5 · 10−5 2 · 10−65 2 · 10−66 1 · 10−3 −7 · 10−5 −1 · 10−13 2 · 10−88 1 · 10−4 −9 · 10−5 −7 · 10−12 −3 · 10−43 Implicit Euler scheme (θ = 1) TABLE III MINIMUM VALUE OF THE DOF (mini,k X k i ) COMPUTED WITH THE SPLITTING ALGORITHM WITH THE CRANK-NICOLSON SCHEME (TOP) AND THE IMPLICIT EULER SCHEME (BOTTOM). stability and the positivity preserving behaviour of the schemes. Tables III and V clearly show that the schemes are not always positivity preserving. In case where the approximated solution is not positive for all t > t0 I also check if it becomes non-negative for larger time ie. if there is t+ > t0 Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 11 of 19 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... δt h≈1.6·10−2 h≈8.2·10−3 h ≈4.1·10−3 5 · 10−1 - - - 2 · 10−1 - - - 1 · 10−1 - - - 4 · 10−2 - - - 2 · 10−2 0.24 - - 1 · 10−2 0.07 0.10 - 4 · 10−3 0.172 t0 t0 2 · 10−3 0.202 0.044 t0 1 · 10−3 0.210 0.079 t0 1 · 10−4 0.2140 0.0939 0.0297 Crank-Nicolson scheme (θ = 1/2) δt h≈1.6·10−2 h≈8.2·10−3 h≈4.1·10−3 5 · 10−1 1 1 1 2 · 10−1 t0 t0 t0 1 · 10−1 t0 t0 t0 4 · 10−2 t0 t0 t0 2 · 10−2 t0 t0 t0 1 · 10−2 t0 t0 t0 4 · 10−3 0.072 t0 t0 2 · 10−3 0.132 t0 t0 1 · 10−3 0.173 0.029 t0 1 · 10−4 0.2102 0.0874 0.0217 Implicit Euler scheme (θ = 1) TABLE IV POSITIVITY THRESHOLD (t+) COMPUTED WITH THE SPLITTING ALGORITHM AND THE CRANK-NICOLSON SCHEME (TOP) AND THE IMPLICIT EULER SCHEME (BOTTOM). such that Xki ≥ 0,∀i = 1,n for all tk > t+ > t0. The smallest such t+, if it exists, is referred as the positivity threshold, as defined in [36]. Table IV on page 12 and table VI on page 13 give the positivity thresholds computed with and without the splitting technique respectively. For the Crank-Nicolson scheme (θ = 1/2) and the implicit Euler scheme (θ = 1) the positivity is obtained under a specific condition on the time step and the mesh size. For a given mesh size, the time step δt must be bounded from above, but also from below to guarantee that the solution stays positive all along the simulation. In the case of the splitting technique those bounds are more restrictive than in the case of the resolution of the full problem without splitting. Those bounds are also more restrictive in the case of the Crank- Nicolson (θ = 1/2) scheme than in the case of the δt h≈1.6·10−2 h≈8.2·10−3 h≈4.1·10−3 5 · 10−1 −4 · 10−1 −4 · 10−1 −4 · 10−1 2 · 10−1 −1 · 10−1 −1 · 10−1 −1 · 10−1 1 · 10−1 1 · 10−13 1 · 10−13 1 · 10−13 4 · 10−2 1 · 10−21 1 · 10−21 9 · 10−22 2 · 10−2 3 · 10−30 1 · 10−30 1 · 10−30 1 · 10−2 −6 · 10−5 1 · 10−42 9 · 10−43 4 · 10−3 −3 · 10−4 3 · 10−64 5 · 10−65 2 · 10−3 −3 · 10−4 −8 · 10−11 3 · 10−87 1 · 10−3 −2 · 10−4 −7 · 10−10 3 · 10−113 1 · 10−4 −9 · 10−5 −1 · 10−11 −8 · 10−35 Crank-Nicolson scheme (θ = 1/2) δt h≈1.6·10−2 h≈8.2·10−3 h≈4.1·10−3 5 · 10−1 2 · 10−5 2 · 10−5 2 · 10−5 2 · 10−1 1 · 10−7 1 · 10−7 1 · 10−7 1 · 10−1 6 · 10−10 5 · 10−10 5 · 10−10 4 · 10−2 1 · 10−15 1 · 10−15 1 · 10−15 2 · 10−2 6 · 10−22 5 · 10−22 5 · 10−22 1 · 10−2 1 · 10−30 7 · 10−31 6 · 10−31 4 · 10−3 −2 · 10−5 2 · 10−47 8 · 10−48 2 · 10−3 −4 · 10−5 2 · 10−64 2 · 10−65 1 · 10−3 −7 · 10−5 −3 · 10−13 2 · 10−87 1 · 10−4 −9 · 10−5 −7 · 10−12 −1 · 10−42 Implicit Euler scheme (θ = 1) TABLE V MINIMUM VALUE OF THE DOF (mini,k X k i ) COMPUTED WITHOUT THE SPLITTING ALGORITHM AND THE CRANK-NICOLSON SCHEME (TOP) AND IMPLICIT EULER SCHEME (BOTTOM). implicit Euler scheme (θ = 1). Refining the mesh results in less restrictions on the time step but also lead to additional computational time. With the Crank-Nicolson scheme (θ = 1/2), for a given mesh size, if δt is too big, there is no threshold of positivity in tk ∈]t0, t1] and the computed solution is not non-negative all along the simulation. For θ = 1/2 and θ = 1, still with a given mesh size, if δt is too small, the simulations showed that there is a threshold of positivity t+ such that the approximated solution becomes non-negative for tk ≥ t+. The thresholds of positivity slightly depend on the time step and tend to increase when the time step δt is decreased. The computations clearly showed that the positivity thresholds diminish with the mesh size h (see for example [36]). Altogether, the positivity of the approximated Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 12 of 19 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... δt h≈1.6·10−2 h≈8.2·10−3 h≈4.1·10−3 5 · 10−1 - 1 1 2 · 10−1 0.8 0.8 0.8 1 · 10−1 t0 t0 t0 4 · 10−2 t0 t0 t0 2 · 10−2 t0 t0 t0 1 · 10−2 0.08 t0 t0 4 · 10−3 0.188 t0 t0 2 · 10−3 0.208 0.050 t0 1 · 10−3 0.213 0.083 t0 1 · 10−4 0.2143 0.0942 0.0300 Crank-Nicolson scheme (θ = 1/2) δt h≈1.6·10−2 h≈8.2·10−3 h≈4.1·10−3 5 · 10−1 t0 t0 t0 2 · 10−1 t0 t0 t0 1 · 10−1 t0 t0 t0 4 · 10−2 t0 t0 t0 2 · 10−2 t0 t0 t0 1 · 10−2 t0 t0 t0 4 · 10−3 0.0720 t0 t0 2 · 10−3 0.1380 t0 t0 1 · 10−3 0.1760 0.0290 t0 1 · 10−4 0.2104 0.0877 0.0221 Implicit Euler scheme (θ = 1) TABLE VI POSITIVITY THRESHOLD (t+) COMPUTED WITHOUT THE SPLITTING ALGORITHM AND THE CRANK-NICOLSON SCHEME (TOP) AND THE IMPLICIT EULER SCHEME (BOTTOM). solution is obtained at the expense of the compu- tational cost, but for a given mesh size h computa- tions performed with too small time step can also lead to a loss of positivity for small tk. In [36] (and references therein), Thomée showed that threshold values of tk > 0 may exist such that X(t) > 0 when t > tk. At this stage, one may wonder how each term of the splitting behaves in terms of positivity preser- vation. The reaction term is approximated using an exact scheme, so obviously the positivity of the solution is preserved. What about the diffusion and the advection term ? 2) Positivity of the pure diffusion problem: Here I set v = (0, 0) and ρ = 0, while keeping all others parameters to the same values as previously. Table VII clearly shows that the Crank-Nicolson scheme (θ = 1/2) is positivity preserving under a δt h≈1.6·10−2 h≈8.2·10−3 h ≈4.1·10−3 5 · 10−1 −5 · 10−1 −5 · 10−1 −5 · 10−1 2 · 10−1 −2 · 10−1 −2 · 10−1 −2 · 10−1 1 · 10−1 4 · 10−15 4 · 10−15 4 · 10−15 4 · 10−2 6 · 10−23 5 · 10−23 4 · 10−23 2 · 10−2 2 · 10−31 8 · 10−32 6 · 1032 1 · 10−2 −4 · 10−5 1 · 10−43 6 · 1043 4 · 10−3 −3 · 10−4 4 · 10−65 5 · 10−66 2 · 10−3 −3 · 10−4 −5 · 10−11 5 · 10−88 1 · 10−3 −2 · 10−4 −6 · 10−10 8 · 10−114 1 · 10−4 −9 · 10−5 −1 · 10−11 −1 · 10−34 Crank-Nicolson scheme (θ = 1/2) δt h≈1.6·10−2 h≈8.2·10−3 h≈4.1·10−3 5 · 10−1 4 · 10−6 4 · 10−6 4 · 10−5 2 · 10−1 2 · 10−8 2 · 10−8 2 · 10−8 1 · 10−1 2 · 10−11 2 · 10−11 2 · 10−11 4 · 10−2 5 · 10−17 5 · 10−17 5 · 10−17 2 · 10−2 3 · 10−23 2 · 10−23 2 · 10−23 1 · 10−2 1 · 10−31 4 · 10−32 3 · 10−32 4 · 10−3 −2 · 10−5 1 · 10−48 6 · 10−49 2 · 10−3 −4 · 10−5 2 · 10−65 2 · 10−66 1 · 10−3 −7 · 10−5 −3 · 10−13 2 · 10−88 1 · 10−4 −9 · 10−5 −7 · 10−12 −2 · 10−42 Implicit Euler scheme (θ = 1) TABLE VII MINIMUM VALUE OF THE DOF (mini,k X k i ) COMPUTED FOR THE PURE DIFFUSION PROBLEM WITH THE CRANK-NICOLSON SCHEME (TOP) AND THE IMPLICIT EULER SCHEME (BOTTOM). CFL-like condition with upper and lower bounds, like in the previous test. The implicit Euler scheme (θ = 1) seems to be more favorable, since it preserves the positivity even for big values of the time step. For both the Crank-Nicolson (θ = 1/2) and implicit Euler (θ = 1) schemes, the approx- imated solution suffers from a loss of positivity for small values of tk when the time step is too small. According to table VIII, there are positivity thresholds, like in [36] which indeed deals with the heat equation. 3) Positivity of the pure advection problem: Here I set σ = 0 and ρ = 0, while keeping all others parameters to the same values as in the first test. Table IX shows that none of the computations performed gave a non negative solutions, even though the minimum value of the dof can be really close to zero for small mesh sizes. Besides, I did Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 13 of 19 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... δt h≈1.6·10−2 h≈8.2·10−3 h≈4.1·10−3 5 · 10−1 - - - 2 · 10−1 0.8 0.8 0.8 1 · 10−1 t0 t0 t0 4 · 10−2 t0 t0 t0 2 · 10−2 t0 t0 t0 1 · 10−2 0.1 t0 t0 4 · 10−3 0.2 t0 t0 2 · 10−3 0.216 0.054 t0 1 · 10−3 0.219 0.085 t0 1 · 10−4 0.2203 0.0956 0.0303 Crank-Nicolson scheme (θ = 1/2) δt h≈1.6· 10−2 h≈8.2·10−3 h≈4.1·10−3 5 · 10−1 t0 t0 t0 2 · 10−1 t0 t0 t0 1 · 10−1 t0 t0 t0 4 · 10−2 t0 t0 t0 2 · 10−2 t0 t0 t0 1 · 10−2 t0 t0 t0 4 · 10−3 0.084 t0 t0 2 · 10−3 0.148 t0 t0 1 · 10−3 0.184 0.032 t0 1 · 10−4 0.2165 0.0893 0.0224 Implicit Euler scheme (θ = 1) TABLE VIII POSITIVITY THRESHOLD (t+) COMPUTED FOR THE PURE DIFFUSION PROBLEM WITH THE CRANK-NICOLSON SCHEME (TOP) AND THE IMPLICIT EULER SCHEME (BOTTOM). not observe any positivity threshold. The approx- imated solution stays non positive all along the simulation. However I run additional simulations with even smaller mesh size (h ≈ 2.0 × 10−3 and δt = 10−4). This time the computed solution was positive at the beginning of the simulation (before t− = 1.9 × 10−3), pointing the existence of a threshold of negativity, to finally reaching a negative minimum values of dof (around −10−44). Unfortunately, this threshold of negativity is really small compared to the ending time of the compu- tation (t1 = 1), while the computational time was reaching more than 14 hours (Intel R©CoreTMi7- 7820HQ at 2.9 GHz, RAM 32 GB) for both the Crank-Nicolson and the implicit Euler schemes. In fact it is well known that for the advection term the solution can be polluted by overshoot and undershoot oscillations near a discontinuity or a δt h≈1.6·10−2 h≈8.2·10−3 h≈4.1·10−3 5 · 10−1 −2 · 10−1 −2 · 10−1 −2 · 10−1 2 · 10−1 −4 · 10−2 −4 · 10−2 −4 · 10−2 1 · 10−1 −4 · 10−3 −2 · 10−3 −1 · 10−3 4 · 10−2 −1 · 10−3 −4 · 10−7 −5 · 10−7 2 · 10−2 −1 · 10−3 −6 · 10−8 −2 · 10−13 1 · 10−2 −1 · 10−3 −4 · 10−8 −3 · 10−28 4 · 10−3 −1 · 10−3 −3 · 10−8 −1 · 10−28 2 · 10−3 −1 · 10−3 −3 · 10−8 −1 · 10−28 1 · 10−3 −1 · 10−3 −3 · 10−8 −1 · 10−28 1 · 10−4 −1 · 10−3 −3 · 10−8 −1 · 10−28 Crank-Nicolson scheme (θ = 1/2) δt h≈1.6·10−2 h≈8.2·10−3 h≈4.1·10−3 5 · 10−1 −1 · 10−4 −1 · 10−10 −8 · 10−34 2 · 10−1 −2 · 10−4 −6 · 10−10 −3 · 10−32 1 · 10−1 −4 · 10−4 −1 · 10−9 −2 · 10−31 4 · 10−2 −6 · 10−4 −4 · 10−9 −9 · 10−31 2 · 10−2 −8 · 10−4 −7 · 10−9 −2 · 10−30 1 · 10−2 −9 · 10−4 −1 · 10−8 −5 · 10−30 4 · 10−3 −1 · 10−3 −1 · 10−8 −8 · 10−30 2 · 10−3 −1 · 10−3 −2 · 10−8 −2 · 10−29 1 · 10−3 −1 · 10−3 −2 · 10−8 −4 · 10−29 1 · 10−4 −1 · 10−3 −3 · 10−8 −9 · 10−29 Implicit Euler scheme (θ = 1) TABLE IX MINIMUM VALUE OF THE DOF (mini,k X k i ) COMPUTED FOR THE PURE ADVECTION PROBLEM WITH THE CRANK-NICOLSON SCHEME (TOP) AND THE IMPLICIT EULER SCHEME (BOTTOM). sharp layer, see [34], [33], [35], [30]. For low order accurate spacial approximations one can prove the positivity preserving property of the scheme [33]. But for high order schemes slopes limiters are often required to guarantee the positivity of the approximated solution. When slope limiters are used, explicit time schemes seem to be suitable for the advection [6]. However, in the next section we will only privilege a numerical scheme that is unconditionally stable, i.e. the Crank-Nicolson scheme, that is a two-order scheme. V. APPLICATION TO THE SIMULATION OF ROOT SYSTEM GROWTH In this section, I apply the previous DG-splitting approach to solve numerically the C-Root model. First, I detail the parameters used for the simula- tions, then, I present and validate the results of the Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 14 of 19 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... simulations. A. The C-Root parameters for Eucalyptus root growth The parameters and operators’ coefficients are chosen based on the previous calibration done in [2]. The diffusion coefficient, σ, is build using the following Gaussian function fα,µ(x,y) = α √ 2π exp ( − (r(x,y) −µ)2 2 ) where r(x,y) = √ (x−x0)2 + (y −y0)2 and (x0,y0) ∈ Ω =] − L,L[. The function fα,µ(x,y) depends on two real and positive parameters: α, related to the maximum amplitude of fα,µ, and µ, the distance from (x0,y0) to the point where the function fα,µ reaches its maximum. The diffusion tensor is taken such that σ(x,y) = fαd,µd (x,y) ( 1 0 0 1 ) , for all (x,y) ∈ Ω, and αd, µd ∈ R+ are given parameters. The advection vector is taken such that v(x,y) = (0,−v0)T , for all (x,y) ∈ Ω, with v0 a positive constant. The reaction term is constant in space and splited into two contributions: βr and µr, the branching and mortality rates, respectively. That is ρ = βr −µr ∈ R. The branching rate, βr, is estimated from biologi- cal knowledge: it is equal to zero before 9 months and equal to 1/3 after, since no roots die before 9 month. However, for the following simulations we will not distinguish the contribution of βr and µr, so that the reaction term will only be described by the parameter ρ. Fig. 6. Density of apices computed at t = 6, t = 12, t = 18 and t = 24 months (from the left to the right and from the top to the bottom). B. Some simulations For the simulation the initial solution is chosen equal to the following function: u0(x,y) = A [ exp(b(1 −x)) (exp(−b(1 −x)) + exp(b(1 −x))) − exp(b(−1 −x)) (exp(−b(−1 −x)) + exp(b(−1 −x))) ] × [ exp(b(1 −y)) (exp(−b(1 −y)) + exp(b(1 −y))) − exp(b(−1 −y)) (exp(−b(−1 −y)) + exp(b(−1 −y))) ] with A = 2 · 10−4 and b = 1. The parameters’ values µr, αd , µd are estimated using the code Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 15 of 19 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... Fig. 7. L2-error with respect to the solution obtained with the mesh of size h ≈ 9.87 × 10−2 and δt = 10−3 computed at t = 6, t = 12, t = 18 and t = 24 months and plotted against the mesh size. described in [2]. I run the simulations from t0 = 1 to t1 = 24 months, with L = 13. The simulations are performed for different values of the mesh size. Fig. 6, page 15, shows the solution computed at four different stages of the root system develop- ment. One can notice the diffusion of the apices in the soil and also the transport of the apices from the top to the bottom of the soil layer. Since there is no analytic solution, the convergence of the computation is evaluated by measuring the L2- errors with respect to the approximated solution computed with the finest mesh (h ≈ 8.97×10−2) and with δt = 10−3. The curves of the errors against the mesh size are plotted on figure 7 and clearly show that the DG-splitting algorithm converges with a convergence rate of almost two. However, one can note that the mesh sizes and the time steps chosen for the simulations presented here might not be small enough. The positivity of the solution is not preserved at all times and the full convergence might not be acheived. Un- fortunatly, refining the mesh sizes and the time steps can lead to prohibitive computational time as shown on table X. On top of that simulation of root system growth can last for a long period of time, particularly for trees. Finally, this application shows promising results for future simulations of h (≈) δt = 10−1 δt = 10−2 δt = 10−3 1.44 1 s. 7 s. 66 s. 7.18 × 10−1 16 s. 46 s. 6 min. 3.59 × 10−1 70 s. 3.5 min. 27 min. 1.79 × 10−1 7 min. 17 min. 2 h. 8.97 × 10−2 50 min. 2h30 9 h. TABLE X COMPUTATIONAL TIMES FOR THE SIMULATIONS OF A ROOT SYSTEM GROWTH PERFORMED (WITH THE PROCESSOR INTEL R©CORETM I7-7820HQ AT 2.9 GHZ, RAM 32 GB) BETWEEN t = 1 AND t = 24 MONTHS WITH THE DG-SPLITTING ALGORITHM AND THE CRANK-NICOLSON SCHEME (θ = 1/2). the root system growth, provided that the compu- tational cost is not limiting. Further simulations requiring much more computational power has to be done to check if the convergence is acheived. This application also point out the difficulties related to the rigorous simulation validation in realistic test-cases of root system growth. VI. CONCLUSION In this work, a discontinuous Galerkin approxi- mation method based on unstructured mesh com- bined with operator splitting has been described, implemented and tested, to solve an advection- diffusion-reaction equation used to model the growth of root systems. The code has been val- idated in a simple test case for which an analytic expression of the solution is known. The compu- tations showed that the method convergences with a convergence rate of two in space with P 1-finite elements. A convergence rate of one and two in time were obtained for respectively the implicit Euler scheme and the Crank-Nicolson scheme both with and without the splitting technique. The computations of those convergence rates required the use of fine mesh grids. For the explicit Euler scheme, such fine mesh computations were not performed since they require really small time steps to fulfill the CFL condition, resulting in huge additional computational cost. Indeed the computational time of the DG-splitting algorithm behaves like 1/δt and 1/h2 where δt and h are respectively the time step and the mesh size. Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 16 of 19 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... Similarly, the positivity of the approximated solution is obtained at the expense of the com- putational time since it requires meshes of small size and small time steps. In fact, there is a CFL- like condition for positivity that has to be fulfilled to guarantee the positivity of the approximated solution. But for a given mesh size computations performed with too small time step can also lead to a loss of positivity at the beginning of the computation [36]. In that cases, the computations showed that there is a positivity threshold in time after which the solution becomes positive. This positivity threshold clearly appeared to diminish with the mesh size. This behavior is specific to the diffusion term. For the advection term, the computations also showed that the positivity of the solution can be preserved, but only at the beginning of the simulation and it required a really small mesh size and time step leading to huge computational time. Further studies in terms of numerical analysis has to be done in that direction. I also performed a more realistic simulation of root system growth. The computations showed that the algorithm converged but additional simulations with smaller time steps and mesh sizes might be performed to recover the full convergence order and positivity. Validation of the computation, but above all the computational time appeared to be the major limitations of the root growth simulation based on the C-Root model, particularly when it comes to deal with trees for which the life span is rather a long period of time. Further improvements on the numerical method has to be done so that the scheme preserves the positivity of the approxi- mated solution under acceptable CFL conditions in terms of computational time. However, our work shows promising results for the simulation of the C-Root model which appears to be an appropriate methodology for future improvements, like root- soil coupling or nonlinear terms arising to handle competition phenomena. ACKNOWLEDGMENT The author would like to thank Y. Dumont (CIRAD, University of Pretoria) for discussions and valuable comments about the numerical schemes. REFERENCES [1] A. Bonneu, Y. Dumont, H. Rey, C. Jourdan and T. Four- caud, A minimal continuous model for simulating growth and development of plant root systems, Plant and Soil, Springer, 2012, 354, 211-22. https://doi.org/10.1007/ s11104-011-1057-7. [2] E. Tillier and A. Bonneu, Operator splitting for solving C-Root, a minimalist and continuous model of root system growth, Plant Growth Modeling, Simulation, Vi- sualization and Applications (PMA), 2012 IEEE Fourth International Symposium on, 2012, 396-402.https://doi. org/10.1109/PMA.2012.6524863. [3] E. Peynaud, T. Fourcaud and Y. Dumont Numerical resolution of the C-Root model using Discontinuous Galerkin methods on unstructured meshes: application to the simulation of roo system growth, 2016 IEEE International Conference on Functional-Structural Plant Growth Modeling, Simulation, Visualization and Appli- cations (FSPMA), IEEE, FSPMA. Qingdao: IEEE, 2016, 158-166. https://doi.org/10.1109/FSPMA.2016.7818302. [4] H.-G. Roos, M. Stynes and L. Tobiska, Robust numerical methods for singularly perturbed differential equations: convection-diffusion-reaction and flow problems Springer Science & Business Media, 2008, 24. https://doi.org/10. 1007/978-3-540-34467-4. [5] X. Zhang and CW. Shu, Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 2011, 467, 2752-2776. https://doi.org/10.1098/rspa.2011.0153. [6] W. Hundsdorfer and J. G. Verwer, Numerical solution of time-dependent advection-diffusion-reaction equations. Springer Science & Business Media, 2013, 33. DOI 10.1007/978-3-662-09017-6 [7] G. Strang, On the construction and comparison of difference schemes. SIAM Journal on Numerical Anal- ysis, SIAM. 5, 506-517. 1968. https://doi.org/10.1137/ 0705041. [8] N. Janenko, The method of fractional steps. Springer. 1971. https://doi.org/10.1007/978-3-642-65108-3. [9] A. Chertock and A. Kurganov, On splitting-based numerical methods for convection-diffusion equations Numerical methods for balance laws, Aracne Editrice Srl Rome, 2010, 24, 303-343. [10] D. Lanser and J.G. Verwer, Analysis of operator splitting for advection-diffusion-reaction problems from air pollution modelling, Journal of computational and applied mathematics, Elsevier, 111, 1-2, 1999, 201-216. https://doi.org/10.1016/S0377-0427(99)00143-0. [11] J. Kačur, B. Malengier, M. Remešı́ková, Convergence of an operator splitting method on a bounded domain for a convection-diffusion-reaction system, Journal of Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 17 of 19 https://doi.org/10.1007/s11104-011-1057-7 https://doi.org/10.1007/s11104-011-1057-7 https://doi.org/10.1109/PMA.2012.6524863 https://doi.org/10.1109/PMA.2012.6524863 https://doi.org/10.1109/FSPMA.2016.7818302 https://doi.org/10.1007/978-3-540-34467-4 https://doi.org/10.1007/978-3-540-34467-4 https://doi.org/10.1098/rspa.2011.0153 https://doi.org/10.1137/0705041 https://doi.org/10.1137/0705041 https://doi.org/10.1007/978-3-642-65108-3 https://doi.org/10.1016/S0377-0427(99)00143-0 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... Mathematical Analysis and Applications. 348, 894-914, 2008. https://doi.org/10.1016/j.jmaa.2008.08.017. [12] M. Remešı́ková, Numerical solution of two-dimensional convection-diffusion-adsorption problems using an operator splitting scheme, Applied mathematics and computation, 184, 116-130, 2007. https: //doi.org/10.1016/j.amc.2005.06.018. [13] J.C. Chrispell, V. Ervin V and E. Jenkins, A fractional step θ-method for convection-diffusion problems, Journal of Mathematical Analysis and Applications, 333, 204- 218, 2007.https://doi.org/10.1016/j.jmaa.2006.11.059. [14] S. Ganesan and L. Tobiska, Operator-splitting finite element algorithms for computations of high-dimensional parabolic problems, Applied Mathematics and Computa- tion, Elsevier, 219, 2013. 6182-6196. https://doi.org/10. 1016/j.amc.2012.12.027. [15] M. Wheeler and C. Dawson, An operator-splitting method for advection-diffusion-reaction problems, MAFELAP Proceedings, 6, 463-482, 1987. [16] R. Anguelov, C. Dufourd and Y. Dumont, Simulations and parameter estimation of a trap-insect model using a finite element approach, Mathematics and Computers in Simulation, 2017, 133, 47-75. https://doi.org/10.1016/ j.matcom.2015.06.014. [17] C. Dufourd and Y. Dumont, Impact of environmental factors on mosquito dispersal in the prospect of sterile insect technique control, Computers & Mathematics with Applications, Elsevier, 2013, 66, 1695-1715. https://doi. org/10.1016/j.camwa.2013.03.024. [18] V. Girault, B. Rivière and M.F. Wheeler, A splitting method using discontinuous Galerkin for the transient incompressible Navier-Stokes equations, ESAIM: Math- ematical Modelling and Numerical Analysis, EDP Sciences,39, 1115-1147, 2005. https://doi.org/10.1051/ m2an:2005048. [19] J. Zhu, Y.T. Zhang, S.A. Newman and M. Alber, Application of discontinuous Galerkin methods for reaction-diffusion systems in developmental biology, Journal of Scientific Computing, Springer, 2009, 40, 391- 418. https://doi.org/10.1007/s10915-008-9218-4. [20] N. Ahmed, G. Matthies and L. Tobiska, Finite element methods of an operator splitting applied to population balance equations Journal of Computational and Applied Mathematics, Elsevier. 236, 1604-1621, 2011. https://doi. org/10.1016/j.cam.2011.09.025. [21] R. Zhang, J. Zhu, A.F Loula and X. Yu, Operator splitting combined with positivity-preserving discontinuous Galerkin method for the chemotaxis model, Journal of Computational and Applied Mathematics. 302, 312 - 326, 2016. https: //doi.org/10.1016/j.cam.2016.02.018. [22] L. Edelstein-Keshet, Mathematical models in biology, SIAM, 46, 1988. ISBN 0-89871-554-7. [23] B. Perthame, Transport equations in biology. Springer Science & Business Media, 2006. https://doi.org/10.1007/ 978-3-7643-7842-4. [24] B. Perthame, Parabolic equations in biology. Growth, reaction, mouvement and diffusion, Springer, 2015. https: //doi.org/10.1007/978-3-319-19500-1 1 [25] A. Ern, J.L Guermond, Theory and practice of finite elements, Springer, 159, 2004. https://doi.org/10.1007/ 978-1-4757-4355-5. [26] V. Volpert, Elliptic Partial Differential Equations: Volume 2: Reaction-Diffusion Equations, Springer, 104, 2014. https://doi.org/10.1007/978-3-0348-0813-2. [27] D. Kuzmin, A guide to numerical methods for transport equations, Friedrich-Alexander-Universitt, Erlangen-Nrnberg, 2010. [28] F. Brezzi, L. D. Marini and E. Süli, Discontinuous Galerkin methods for first-order hyperbolic problems, Mathematical models and methods in applied sciences, World Scientific, 2004, 14, 1893-1903. https://doi.org/10. 1142/S0218202504003866. [29] W. H. Reed and T. R. Hill, Triangular mesh methods for the neutron transport equation, Los Alamos Scientific Lab., N. Mex.(USA), report LA-UR-73-479, 1973. [30] B. Rivière, Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation, Society for Industrial and Applied Math- ematics, 2008. https://doi.org/10.1137/1.9780898717440 [31] J. Oden, I. Babuŝka and C. E. Baumann, A discontinuous hp-finite element method for diffusion problems, Journal of computational physics, Elsevier 146, 491-519, 1998. https://doi.org/10.1006/jcph.1998.6032. [32] D. A. Di Pietro and A. Ern, Mathematical aspects of discontinuous Galerkin methods Springer, 69, 2011. https://doi.org/10.1007/978-3-642-22980-0. [33] X. Zhang and CW. Shu CW. Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 2011, 467, 2752-2776. https://doi.org/10.1098/rspa.2011.0153. [34] B. Cockburn and CW. Shu, Runge-Kutta discontinuous Galerkin methods for convection-dominated problems. Journal of scientific computing, Springer, 2001, 16, 173- 261. https://doi.org/10.1023/A:1012873910884. [35] JS. Hesthaven and T. Warburton, Nodal discontinuous Galerkin methods: algorithms, analysis, and applications. Springer, 2007, 54. https://doi.org/10.1007/978-0-387-72067-8. [36] V. Thomée, On positivity preservation in some finite element methods for the heat equation. International Con- ference on Numerical Methods and Applications, 2014, 13-24. https://doi.org/10.1007/978-3-319-15585-2 2. [37] J. Zhu, YT. Zhang, SA. Newman and M. Alber, Application of discontinuous Galerkin methods for reaction-diffusion systems in developmental biology. Journal of Scientific Computing, Springer, 2009, 40, 391- 418. https://doi.org/10.1007/s10915-008-9218-4. [38] J.-F. Barczi, H. Rey, S. Griffon and C. Jourdan, DigR: a generic model and its open source simulation software to mimic three-dimensional root-system architecture Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 18 of 19 https://doi.org/10.1016/j.jmaa.2008.08.017 https://doi.org/10.1016/j.amc.2005.06.018 https://doi.org/10.1016/j.amc.2005.06.018 https://doi.org/10.1016/j.jmaa.2006.11.059 https://doi.org/10.1016/j.amc.2012.12.027 https://doi.org/10.1016/j.amc.2012.12.027 https://doi.org/10.1016/j.matcom.2015.06.014 https://doi.org/10.1016/j.matcom.2015.06.014 https://doi.org/10.1016/j.camwa.2013.03.024 https://doi.org/10.1016/j.camwa.2013.03.024 https://doi.org/10.1051/m2an:2005048 https://doi.org/10.1051/m2an:2005048 https://doi.org/10.1007/s10915-008-9218-4 https://doi.org/10.1016/j.cam.2011.09.025 https://doi.org/10.1016/j.cam.2011.09.025 https://doi.org/10.1016/j.cam.2016.02.018 https://doi.org/10.1016/j.cam.2016.02.018 https://doi.org/10.1007/978-3-7643-7842-4 https://doi.org/10.1007/978-3-7643-7842-4 https://doi.org/10.1007/978-3-319-19500-1_1 https://doi.org/10.1007/978-3-319-19500-1_1 https://doi.org/10.1007/978-1-4757-4355-5 https://doi.org/10.1007/978-1-4757-4355-5 https://doi.org/10.1007/978-3-0348-0813-2 https://doi.org/10.1142/S0218202504003866 https://doi.org/10.1142/S0218202504003866 https://doi.org/10.1137/1.9780898717440 https://doi.org/10.1006/jcph.1998.6032 https://doi.org/10.1007/978-3-642-22980-0 https://doi.org/10.1098/rspa.2011.0153 https://doi.org/10.1023/A:1012873910884 https://doi.org/10.1007/978-0-387-72067-8 https://doi.org/10.1007/978-3-319-15585-2_2 https://doi.org/10.1007/s10915-008-9218-4 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Emilie Peynaud, Operator splitting and discontinuous Galerkin methods for advection-reaction- ... diversity. Annals of Botany, 2018, 121, 5, 1089-1104, https://doi.org/10.1093/aob/mcy018. [39] L. X. Dupuy, M. Vignes, An algorithm for the simulation of the growth of root systems on deformable domains. Journal of Theoretical Biology, 2012, 310, 164- 174. https://doi.org/10.1016/j.jtbi.2012.06.025. [40] L. Dupuy, P. J. Gregory, A. G. Bengough, Root growth models: towards a new generation of continuous approaches. Journal of experimental botany, Soc Experi- ment Biol, 2010. https://doi.org/10.1093/jxb/erp389. [41] P. Bastian, A. Chavarria-Krauser, C. Engwer, W. Jäger, S. Marnach, M. Ptashnyk, Modelling in vitro growth of dense root networks Journal of theoretical biology, Elsevier, 2008, 254, 99-109. https://doi.org/10.1016/j.jtbi. 2008.04.014. [42] T. Roose, A. Schnepf, Mathematical models of plant-soil interaction Philosophical Transactions of the Royal Soci- ety of London A: Mathematical, Physical and Engineer- ing Sciences, The Royal Society, 2008, 366, 4597-4611, https://doi.org/10.1098/rsta.2008.0198. [43] S. G. Adiku, R. D. Braddock, C. W. Rose, 1996, Simulating root growth dynamics, Environ- mental Software 11 : 99-103. https://doi.org/10.1016/ S0266-9838(96)00041-X. [44] H. Hayhoe, 1981, Analysis of a diffusion model for plant root growth and an application to plant soil-water uptake, Soil Science 131 : 334-343. [45] M. Heinen, A. Mollier, P. De Willigen, 2003 Growth of a root system described as diffusion numerical model and application, Plant and Soil 252 : 251-265. https://doi.org/ 10.1023/A:1024749022761. [46] V. R. Reddy, Ya. A. Pachepsky, 2001, Testing a convective dispersive model of two dimensional root growth and proliferation in a greenhouse experiment with mare plants, Annals of Botany 87 : 759-768. https://doi. org/10.1006/anbo.2001.1409. [47] P. -H. Tournier, F. Hecht, M. Comte, Finite element model of soil water and nutrient transport with root uptake: explicit geometry and unstructured adaptive meshing. Transp. Porous Media 106 (2), 487504 (2015). https://doi.org/10.1007/s11242-014-0411-7. [48] M. Comte, Analysis and Simulation of a Model of Phosphorus Uptake by Plant Roots in Current Research in Nonlinear Analysis: In Honor of Haim Brezis and Louis Nirenberg, Rassias, T. M. (Ed.), Springer Inter- national Publishing, 2018, 85-97. https://doi.org/10.1007/ 978-3-319-89800-1 4. [49] F. Gérard, Cé Blitz-Frayret, P. Hinsinger, L. Pagès, Modelling the interactions between root system architecture, root functions and reactive transport processes in soil Plant and Soil, 2017, 413, 161-180. https://doi.org/10.1007/s11104-016-3092-x. [50] H. Brezis, Functional analysis, Sobolev spaces and partial differential equations. Springer Science & Business Media, 2010. https://doi.org/10.1007/ 978-0-387-70914-7. [51] P. R. Amestoy, I. S. Duff, J. Koster and J.-Y. L’Excellent, A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM Journal of Matrix Analysis and Applications, Vol 23, No 1, pp 15-41 (2001). https: //doi.org/10.1137/S0895479899358194. [52] P. R. Amestoy, A. Guermouche, J.-Y. L’Excellent and S. Pralet, Hybrid scheduling for the parallel solution of linear systems. Parallel Computing Vol 32 (2), pp 136- 156 (2006). https://doi.org/10.1016/j.parco.2005.07.004. Biomath 7 (2018), 1812037, http://dx.doi.org/10.11145/j.biomath.2018.12.037 Page 19 of 19 https://doi.org/10.1093/aob/mcy018 https://doi.org/10.1016/j.jtbi.2012.06.025 https://doi.org/10.1093/jxb/erp389 https://doi.org/10.1016/j.jtbi.2008.04.014 https://doi.org/10.1016/j.jtbi.2008.04.014 https://doi.org/10.1098/rsta.2008.0198 https://doi.org/10.1016/S0266-9838(96)00041-X https://doi.org/10.1016/S0266-9838(96)00041-X https://doi.org/10.1023/A:1024749022761 https://doi.org/10.1023/A:1024749022761 https://doi.org/10.1006/anbo.2001.1409 https://doi.org/10.1006/anbo.2001.1409 https://doi.org/10.1007/s11242-014-0411-7 https://doi.org/10.1007/978-3-319-89800-1_4 https://doi.org/10.1007/978-3-319-89800-1_4 https://doi.org/10.1007/s11104-016-3092-x https://doi.org/10.1007/978-0-387-70914-7 https://doi.org/10.1007/978-0-387-70914-7 https://doi.org/10.1137/S0895479899358194 https://doi.org/10.1137/S0895479899358194 https://doi.org/10.1016/j.parco.2005.07.004 http://dx.doi.org/10.11145/j.biomath.2018.12.037 Introduction The model Modelling root growth with PDE: the C-Root model The weak problem The positivity preserving property of the solution Approximation of the model The operator splitting technique The advection step: DG upwind scheme The diffusion step The reaction step Validation of the splitting algorithm with a simple test case Description of the simple test-case Numerical validation and convergence Some comments on the positivity Positivity of the full problem Positivity of the pure diffusion problem Positivity of the pure advection problem Application to the simulation of root system growth The C-Root parameters for Eucalyptus root growth Some simulations Conclusion References