Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Moving Mesh Discontinuous Galerkin Methods for PDEs with Traveling Waves

Moving Mesh Discontinuous Galerkin Methods for PDEs with Traveling Waves In this paper, a moving mesh discontinuous Galerkin (dG) method is developed for nonlinear partial differential equations (PDEs) with traveling wave solutions. The moving mesh strategy for one dimensional PDEs is based on the rezoning approach which decouples the solution of the PDE from the moving mesh equa- tion. We show that the dG moving mesh method is able to resolve sharp wave fronts and wave speeds accurately for the optimal, arc-length and curvature moni- tor functions. Numerical results reveal the efficiency of the proposed moving mesh dG method for solving Burgers’, Burgers’-Fisher and Schlogl(Nagumo) ¨ equations. 1 Introduction The discontinuous Galerkin (dG) method is one of the most powerful discretization techniques for solving partial differential equations (PDEs) [2, 11], especially for con- vection dominated problems, exhibiting localized phenomena like sharp traveling wave fronts, internal and boundary layers [9, 13]. The dG method has been applied for this kind of singularly perturbed linear and nonlinear PDEs extensively using h-adaptive (refinement and coarsening in space), p-adaptive (enrichment of the local polynomial degree), hp-adaptive and space-time adaptive methods in the last two decades. An- other approach to deal with these kind of problems, is the r-method or moving mesh method. In the moving mesh method the grid points are relocated in the regions where the solution shows rapid variation, while keeping the number of the nodes fixed. The dG discretization is very flexible, since there is no continuity requirement between the inter-element boundaries, which makes it suitable as a moving mesh method on ir- regular meshes. Most of the studies with moving mesh methods are limited to finite difference and continuous finite element discretization [14]. There are only few pub- lications dealing with dG moving mesh method. They include the interior penalty dG method for preprocessing the solutions of steady state diffusion-convection-reaction equations [1], and the local dG moving mesh method for hyperbolic conservation laws [10]. In this paper we develop an adaptive dG moving mesh method for one dimensional arXiv:1602.05225v1 [math.NA] 31 Jan 2016 semi-linear differential equations with traveling wave solutions of the form u = e u f (u; u ); (x; t)2 W(t ; T ] (1a) t xx x 0 f u(x ; t) = u ; u(x ; t) = u ; t 2 (t ; T ] (1b) L L R R 0 f u(x; t ) = u ; x2 W; (1c) 0 0 where W = [x ; x ] R, t and T are the initial and final time instances, respectively, L R 0 and e denotes the diffusion coefficient. The model equation (1) becomes the Burgers equation with f (u; u ) = uu [12] , Burgers’-Fisher equation with f (u; u ) = a uu + x x x x b u(u 1) [12] and the Schlogl ¨ or Nagumo equation with f (u; u ) = u(1 u)(1b)=d [14]. A moving mesh method has three main components; the discretization of the phys- ical PDE, mesh strategy using monitor functions and discretization of the mesh equa- tion. The discretization of physical PDE is either coupled with the moving mesh equa- tion or separated. In the quasi-Lagrangian approach, a large system of the discretized PDE and moving mesh equation are solved simultaneously by the standard ordinary differential equation (ODE) solvers. Instead, we use the rezoning approach by solving alternately the PDE and mesh equation, which allows more flexibility; mesh genera- tion can be coded separately and embedded in the solution of the PDE. Since the mesh is updated at each time step, the physical PDE has to be discretized at the next time step on the new mesh. We use the static rezoning approach with the same number of points at each time step [7] in contrast to the dynamic rezoning [8] where the number of mesh points is changed at every time step. Therefore in the static rezoning approach the solutions from old to new mesh have to be interpolated. The paper is organized as follows. In the next section we describe briefly the dG method for the 1D model problem (1) on a uniform fixed mesh. Moving mesh adaption strategy and the adaptive moving mesh dG algorithm is presented in Section 3. Nu- merical results are given in Section 4 to demonstrate the effectiveness of the proposed method. 2 Discretization of the problem on a fixed mesh Before giving the moving mesh strategy in Section 3, in this section we describe the dG discretization of the model problem (1) on a fixed uniform mesh T : x = x + nh; n = 0; 1;:::; N ; (2) h n L I consisting of N elements (sub-intervals) I = [x ; x ], n = 1; 2;:::; N , and with the I n n1 n I fixed mesh size h = (x x )=N . R L I 2.1 Space discretization by discontinuous Galerkin method We use for the space discretization of the model problem (1) on a fixed mesh (2) the symmetric interior penalty Galerkin (SIPG) method [2, 11] which is a member of the 2 family of dG methods. The dG methods use the space of piecewise discontinuous polynomials of degree at most k: V =fv : vj 2 P (I ) ; 8n = 1;:::; N g; h I k n I where P (I ) is the space of polynomials of degree at most k on an interval I . Since the k n n functions in V are discontinuous at the inter-element nodes, we define the jump and average of a piecewise function v at the endpoints of I , n = 1;:::; N 1, respectively, n I as depicted in Figure 1, + + [v(x )] = v(x ) v(x ) ; fv(x )g = (v(x )+ v(x )); (3) n n n n n n with v(x ) = lim v(x) ; v(x ) = lim v(x): (4) n n x7!x x7!x n n On the boundary nodes, the jump and average are defined as + + + [v(x )] =v(x ); fv(x )g = v(x ); [v(x )] = v(x ); fv(x )g = v(x ): (5) 0 0 N N 0 0 I N I N I I v(x ) v| [v(x )] v| n+1 v(x ) x x x x x x N −1 N 0 1 n−1 n+1 I I Figure 1: Jump and limit terms of a piecewise discontinuous function v(x). The SIPG scheme is constructed by multiplying the continuous (the solution u is sufficiently smooth at the end points of I ) equation (1) by a test function v 2 V and n h integrating by parts on each element I , n = 1;:::; N : n I Z Z Z x x x n n n u vdx+ e u v dxe u (x )v(x )+e u (x )v(x )+ f (u; u )vdx = 0: t x x x n x n1 x n1 x x x n1 n1 n1 By adding all N equations, and using the definition of the jumps (3) and (5), we obtain Z Z Z N N I x x x I n n n u vdx+ e u v dx+ f (u; u )vdx [e u (x )v(x )] = 0: t x x x x n n å å x x x n1 n1 n1 n=1 n=0 One can verify that for 1 n N 1 [e u (x )v(x )] =fe u (x )g[v(x )]+[e u (x )]fv(x )g: (6) x n n x n n x n n Using the identity (6) and the fact that [e u (x )] = 0 for all 1  n  N 1 (u was x n I sufficiently smooth at the end points of I ), we obtain Z Z Z N N I x x x I n n n u vdx+ e u v dx+ f (u; u )vdx fe u (x )g[v(x )] = 0: (7) t x x x x n n å å x x x n1 n1 n1 n=1 n=0 | {z } Additionally, we have [u(x )] = 0 for all 1 n N 1. Then, adding the penalizing n I terms and the terms on the boundary nodes, n =f0; N g, to both sides of (7) by keeping them unknown on the left hand side and imposing the boundary conditions u and u L R on the right hand side, leads to the SIPG formulation: Z Z x x R R u vdx+ a(u; v)+ f (u; u )vdx = l(v): (8) t x x x L L In (8), a(u; v) and l(v) denote the symmetric bilinear form and the linear right hand side of the SIPG scheme x I a(u; v) = e u v dx+ fe u (x )g[v(x )]fe v (x )g[u(x )]+ [u(x )][v(x )] ; x x x n n x n n n n n=0 s s l (v ) = u e v (x ) v(x ) + u v(x )e v (x ) ; h h L x 0 0 R N x N I I h h where s is the penalty parameter which should be sufficiently large [11] in order to ensure the coercivity of the bilinear form. Hence, SIPG semi-discrete form of (1) reads as: a.e. t 2 (t ; T ], for all v 2 V , find u := u (x; t)2 V such that 0 f h h h h h Z Z x x R R u (x; t )v dx = u v dx; (9a) h 0 h 0 h x x L L Z Z x x R R (¶ u )v dx+ a(u ; v )+ f (u ; u )v dx = l(v ): (9b) t h h h h h h;x h h x x L L 2.2 Full discretization Let ff g denote the basis functions spanning the dG finite elements space V for i = 1;:::; N and n = 1;:::; N , where N stands for the local dimension depending on the k I k polynomial order k and is given by N = k + 1 in 1D. The local nature of dG methods leads to the basis functions and the approximate solution of the form n I k f (x) ; x2 I n n n f (x) = ; u (t) = u (t)f (x); (10) i h å å i i 0 ; x2 = I n=1 i=1 wherefu (t)g are the time-dependent unknown coefficients. We substitute the second identity in (10) into the system (9b) and we choose v = f for i = 1;:::; N and n = h k 1;:::; N , which leads to the N = N  N dimensional non-linear system of equations I k I of (9b) in matrix-vector form Mu + Su + h(u) d = 0; (11) NN NN where M2 R is the usual mass matrix and S2 R is the stiffness matrix related N N to the bilinear form a(u ; v ). The vectors h(u) 2 R and d 2 R are the non-linear h h vector of unknown coefficients u corresponding to the integral of the nonlinear term in (9b) and the right hand side vector related to the linear form l(v ), respectively. The initial vector u(0) is found by using the equation (9a) and the second identity in (10). 4 We solve the fully discrete system of (1), by applying the backward Euler scheme to the semi-discrete system (11). Let t < t < ::: < t = T be the uniform partition 0 1 J f of the time interval [t ; T ] into J time-steps (t ; t ], j = 1;:::; J, with the step size 0 f j1 j Dt = (T t )=J. Let us denote the approximate coefficient vector u(t) of (11) at the f 0 time t = t by u . Then, the fully discrete formulation of the model (1) is given as: for j N all j = 1;:::; J, find u 2 R such that j j1 u u j j j M + Su + h(u ) d = 0; (12) Dt which is solved by Newton’s method. 3 Adaptive moving mesh method In an adaptive moving mesh method the spatial mesh T in (2) varies with time T (t) : x = x(x ; t); n = 0; 1;:::; N ; (13) h n n I consisting of N elements I = [x ; x ], n = 1; 2;:::; N , of the mesh size h = x I n n1 n I n n x . In (13) the nodes x belong to the fixed and uniform mesh n1 n T : x = ; n = 0; 1;:::; N (14) n I on the auxiliary unit interval W = [0; 1], together with the boundary conditions x(0; t) = x and x(1; t) = x . Thus, in the adaptive moving mesh method, the mesh is adjusted as L R the time progresses in such a way that the mesh sizes h get smaller in the sub-intervals where the error is large, while the mesh sizes h are made coarser in the remaining part of the interval. The error indicator is chosen to relocate the mesh points where the solution shows large variations, based on the equidistribution principle, where a mesh T (t) with the mesh points x = x < x < ::: < x < x = x is determined by h L 0 1 N 1 N R I I satisfying the following relation x N 1 I Z Z r(x; t)dx = = r(x; t)dx: (15) x x 0 N 1 In this way a continuous function r(x; t) > 0 on the interval [x ; x ] can be distributed L R among the sub-intervals I = [x ; x ], n = 1;:::; N . In (15), the function r(x; t) is n n1 n I called the mesh density function, or monitor function, choice of which stands as the key point for an adaptive moving mesh method. The most popular choices for the monitor functions r := r(x; t) are [5] optimal 1=3 r = 1+ ju j ; (16) xx 5  arc-length 1=2 r = 1+ju j ; (17) curvature 1=4 r = 1+ju j ; (18) xx with the intensity parameter ( ) 2=3 a = max 1; ju j dx : (19) xx x x R L L Finding a proper mesh T (t) using the equidistribution condition (15) results in a system of so-called moving mesh partial differential equation (MMPDE) [6, 5] ¶ x 1 ¶ ¶ x = r ; (x; t)2 W (t ; T ]; (20a) c 0 f ¶ t tr ¶x ¶x x(0; t) = x ; x(1; t) = x ; t 2 (t ; T ]: (20b) L R 0 f The system (20) is solved through the central finite differences dx 1 r +r r +r n n+1 n n n1 = (x + x ) (x + x ) (21) n+1 n n n1 dt tr Dx 2 2 for n = 1; 2;:::; N 1, where the spatial nodes x 2 T (t) are the unknown solutions I n h of the nodes x 2T . The relaxation parameter t is specified by the user for adjusting the response time of mesh movement according to the changes of the density function r(x; t). The functions r , n = 0; 1;:::; N , are computed through the discrete form of n I the monitor function r . For instance, we have for the optimal monitor function (16) 1=3 r = 1+ ju j ; n = 0;:::; N ; (22) n xx;n I where the terms u are computed by the central difference approximations using xx;n the solutions at the mesh nodes x . At each time step, further, the discrete monitor functions r are smoothed by weighted averaging [5, Section 1.2]. There are different approaches [14] of the adaptive moving mesh method for solv- ing the fully discrete system (12), called the physical PDE, and the MMPDE (21). One common approach is solving both systems simultaneously using a quasi-Lagrange ap- proach. In this approach, the time derivative term requires a special attention since the mesh is assumed to move in a continuous manner by the time progresses. In such cases, there occur an extra convective term which may cause difficulties. Additionally, the solution of the both systems simultaneously needs the coupling of the systems, as a result the dimension of the system to be solved increases. Another choice, which we use in this paper, is the alternate solution using a rezoning approach. In this case, the physical PDE and the MMPDE are separated from each other, which allows more flex- ibility as the mesh generation can be coded separately and embedded in the solution 6 of the physical PDE. In the static rezoning approach, the change in the spatial mesh is derived in a discrete form similar to the solution [7, 8, 14]. In each time step, first the mesh adaptation is handled using the solution on the old mesh, then the solution is obtained by solving the physical PDE on the newly generated mesh. In the case of dG method, the unknown coefficient vectors have to satisfy u u(t ) in the physical PDE (12) with the discrete mesh density function (22) in the MMPDE (21). In general, it is difficult to use the coefficient vectors u directly to compute the discrete mesh density functions r , unless we use the Lagrange nodal basis. Let us choose the dG basis functions f as the Lagrange nodal basis functions, i = 1;:::; N , n = 1;:::; N . Then, recalling the definitions (4) of the limit terms at a node x , we have the relations n + n+1 u (x ; t) = u (t) ; u (x ; t) = u (t); n = 1;:::; N 1: (23) h h I n N n The relations (23) originate from the fact that dG methods use discontinuous basis functions at the inter-element nodes. As a result, a dG approximation u (x; t) has two traces at an inter-element node x from the neighboring sub-intervals I and I , as n n n+1 shown in Figure 1, which are not the same in general. Using the relations (23) and recalling again the average definition in (3), we can accept the value of the approximate solution u (x ; t) as h n u (x ; t) : =fu (x ; t)g = (u (x ; t)+ u (x ; t)) h n h n h h n n (24) n n+1 = (u (t)+u (t)); n = 1;:::; N 1: In this way, the solutions of the physical PDE (12) will be consistent with the MM- PDE (21), where the computation of the discrete mesh density functions r require the discrete approximations u (x; t) at the inter-element nodes x , n = 0;:::; N . h n I Algorithm 1 Moving mesh algorithm j1 j1 j1 Given the current spatial mesh T (t ) : x < x < ::: < x , co- h j1 0 1 N j1 efficient vector u , the parameter t and time-step size Dt , do for j = 1;:::; J, 1: Compute a temporary coefficient vector u ˜ by solving the physical PDE (12) on the current mesh T (t ), h j1 2: According to (24), calculate the consistent values of u ˜ (x ; t ) using u , h n j 3: Compute the discrete mesh density functions r ˜ using the values u ˜ (x ; t ), n n j j j j 4: Find the mesh T (t ) : x < x < ::: < x by solving the MMPDE (21) for the h j 0 1 N discrete mesh density functions r , j1 5: Interpolate the coefficient vector u to be used in the new mesh T (t ), h j 6: Compute the coefficient vector u by solving the physical PDE (12) on the new mesh T (t ), h j 7: Go to next time step. The general algorithm is summarized in Algorithm 1. We start with an initial mesh T (t ) (possibly a uniform mesh). Then on an arbitrary time-step (t ; t ], we, firstly, h 0 j1 j 7 j solve the physical PDE (12) on the mesh T (t ) for an auxiliary solution u ˜ . After, h j1 following the relation (24), we use the solution u ˜ in the MMPDE (21) to obtain the new mesh T (t ). Finally, we solve the physical PDE (12) on the mesh T (t ) for h j h j the solution u , and we proceed to the next time-step. On the other hand, the known j1 solution vector u will be no more consistent with the new mesh T (t ). This is h j a natural consequence of the alternate solution with rezoning approach. To make the j1 known solution u consistent with the updated spatial mesh T (t ), we interpolate it between the meshes T (t ) and T (t ). The interpolation procedure is as follows: let j1 j h h j1 j1 j1 j j j T (t ) : x < x < ::: < x be the current mesh andT (t ) : x < x < ::: < x j1 j h h 0 1 N 0 1 N I I j1 n ˆ ˆ denote the updated mesh. For the sake of simplicity, let also u := u with entries u being the coefficient vector defined on the current meshT (t ), and u with entries u h j1 denotes the interpolated coefficient vector of u into the new mesh T (t ), i = 1;:::; N , h j k j1 j1 n = 1;:::; N . The local nature of dG methods leads for any x2 I = [x ; x ] to I s s s1 s s u (x; t ) = u f (x); s = 1;:::; N : (25) h j1 I å i i i=1 On the other hand, using the Lagrange basis functions, on an arbitrary element we j j obtain I = [x ; x ] on the new mesh T (t ), uniformly distributed N nodal degrees n n h j k n1 of freedoms x 2 I such that n;d j j j x = x +(d 1)t; u (x ; t ) u ; (26) h j1 n;d n1 n;d j j1 for d = 1;:::; N , and with t = (x x )=(N 1). In other words, the entries u k n n k of the interpolated coefficient vector u are the approximate function values u (x; t ) h j1 at the nodal degrees of freedoms x , n = 1;:::; N , d = 1;:::; N . For any nodal I k n;d degrees of freedom x on the new mesh T (t ), we have to determine the intervals h j n;d j1 j1 j I = [x ; x ] on the current mesh T (t ) such that x 2 I . Then, using the s s h j1 s s1 n;d expansion (25), we will be able to obtain the entries of u as u = u (x ; t ). h j1 i n;i 4 Numerical results In this section we present several numerical examples demonstrating the effectiveness of the adaptive moving mesh dG method. In all of the examples, the traveling wave solutions are computed by the optimal mesh density function (16), but the correspond- ing mesh trajectories are given for the optimal (16), arc-length (17) and curvature (18) mesh density functions. 4.1 Burgers’ equation The first test example is the Burgers’ equation [14, 12] ¶ 1 u = e u u t xx ¶ x 2 8 with homogeneous Dirichlet boundary conditions on the space-time domain (x; t) 2 [0; 1] [0; 1] with the diffusion constant e = 10 . The initial condition is taken as u(x; 0) = sin(2p x) + 0:5 sin(p x), and linear dG basis functions are used. We set the time-step size Dt = 0:005, and we choose the relaxation parameter as t = 10 . Moving mesh solutions in Figure 3 are capable of resolving the sharp gradients of the moving fronts in contrast to the oscillatory solutions on the fixed mesh in Figure 2. In addition, in Figure 3, mesh trajectories for the curvature and arc-length monitor functions show large variations with respect to time, whereas for the optimal monitor function the mesh trajectories are smooth. 1.5 0.5 t=0 −0.5 t=0.2 t=0.4 −1 t=0.6 t=0.8 t=1 −1.5 0 0.2 0.4 0.6 0.8 1 Figure 2: Burgers’: Solutions at t = 0, 0.2, 0.4, 0.6, 0.8, 1 on the uniform fixed mesh with N = 120 elements. 4.2 Burgers’-Fisher equation Consider the Burgers’-Fisher equation [12] ¶ 1 u = u a u +b u(u 1); t xx ¶ x 2 with the exact solution 1 a u(x; t) = 1 tanh (x ct) ; 2 4 in the space-time domain (x; t)2 [1; 0](0:2; 0]. The parameters are a = 24, c = 8 and b = (2a c a )=4. We use quadratic dG basis functions and N = 40 spatial elements. The time step-size is taken as Dt = 0:001. In Figure 4, we give the solutions of Burgers’-Fisher equation at different time instances for the optimal mesh density function, and the moving mesh trajectories ob- tained by different monitor functions. In Table 1 the L -errors between the exact and numerical solutions are tabulated at different times. The results are very close to those in [12] computed with the same settings. u 1 0.9 0.9 0.8 0.8 0.7 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (a) Optimal (b) Arc-Length 1.5 0.9 0.8 0.7 0.6 0.5 0.5 0.4 0.3 t=0 t=0.2 0.2 −0.5 t=0.4 t=0.6 0.1 t=0.8 t=1 −1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x (c) Curvature (d) Solutions Figure 3: Burgers’ equation: (a-c) Moving mesh trajectories with N = 40 elements, (d) solutions at t = 0, 0.2, 0.4, 0.6, 0.8, 1 for the optimal monitor function. Table 1: Burgers’-Fisher equation: L -errors Monitor N t=-0.1 t=-0.05 t=-0.04 t=-0.035 t=-0.03 Optimal 40 4.6e-03 8.4e-03 1.1e-02 1.1e-02 1.3e-02 Arc-Length 40 2.4e-03 4.2e-03 5.2e-03 5.8e-03 6.5e-03 Curvature 40 2.4e-03 4.1e-03 5.1e-03 5.7e-03 6.4e-03 4.3 Schlogl ¨ equation The final example is Schlogl ¨ (Nagumo) equation [3, 14] u = e u f (u); t xx with the Ginzburg-Landau free energy E (u) = jÑuj + F(u) dx (27) t t t 0 −0.02 −0.02 −0.04 −0.04 −0.06 −0.06 −0.08 −0.08 −0.1 −0.1 −0.12 −0.12 −0.14 −0.14 −0.16 −0.16 −0.18 −0.18 −0.2 −0.2 −1 −0.8 −0.6 −0.4 −0.2 0 −1 −0.8 −0.6 −0.4 −0.2 0 x x (a) Optimal (b) Arc-Length 0 1 t=−0.1 −0.02 0.9 t=−0.05 −0.04 0.8 −0.06 0.7 t=−0.04 −0.08 0.6 t=−0.035 −0.1 0.5 t=−0.03 −0.12 0.4 −0.14 0.3 −0.16 0.2 −0.18 0.1 −0.2 −1 −0.8 −0.6 −0.4 −0.2 0 −1 −0.8 −0.6 −0.4 −0.2 0 x x (c) Curvature (d) Solutions Figure 4: Burgers’-Fisher equation: (a-c) Moving mesh trajectories with N = 40 el- ements, (d) solutions at t = -0.1, -0.05, -0.04, -0.035, -0.03 and corresponding exact solutions by solid lines. 2 2 with quartic potential function F(u) = u (3u 4(1+b)u+ 6b) and cubic bi-stable 12d nonlinearity f (u) = u(u 1)(ub) satisfying f (u) = F (u). We consider Schlogl ¨ equation in the space-time domain (x; t)2 (0; 1] [0; 1] with 3 3 the constant parameter values e = 10 , d = 10 and b = 0. The Dirichlet boundary conditions and the initial condition are taken according to the exact solution 1 x ct u(x; t) = 1 tanh p ; 8ed where c = e=2d is the speed of the traveling wave. The solutions are computed with linear and quadratic dG basis functions, and with the time-step size Dt = 0:001. As the number of spatial elements, we take N = 120 and N = 40 to construct a uniform and I I moving mesh, respectively. t t u Figure 5 shows that the steep wave fronts are captured well by the adaptive moving mesh method with linear and quadratic dG basis functions. On the other hand, Figure 6 shows that the numerical solutions with quadratic dG basis functions give the correct wave speed, whereas for the linear case the numerical solutions move faster than the exact solutions. In Table 2, we list the L -errors between the numerical and exact so- lutions at different time instances for different types of monitor functions. The correct wave speed is not captured by linear dG basis functions. Therefore the corresponding L -errors are larger than the quadratic case. Due the discontinuous nature of the dG discretization, higher order basis functions can be implemented in a more flexible way in contrast to continuous finite elements, providing continuity requirement between the inter-element boundaries. The free energy of Schlogl ¨ equation (27) decreases monoton- ically in time. The energy decreasing property is also captured numerically on uniform and moving meshes, as shown in Figure 6 for the conditionally energy stable backward Euler method [4]. 1 1 1 0.9 0.9 0.9 0.8 0.8 0.8 0.7 0.7 0.7 0.6 0.6 0.6 0.5 0.5 0.5 0.4 0.4 0.4 0.3 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.1 0 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x x (a) Optimal (b) Arc-length (c) Curvature 1 1 0.9 0.9 0.9 0.8 0.8 0.8 0.7 0.7 0.7 0.6 0.6 0.6 0.5 0.5 0.5 0.4 0.4 0.4 0.3 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.1 0 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x x (d) Optimal (e) Arc-length (f) Curvature Figure 5: Schlogl ¨ equation: Moving mesh trajectories with N = 40 elements; (a)-(c) linear dG basis functions, (d)-(f) quadratic dG basis functions. 5 Conclusions In this paper we have developed an adaptive discontinuous Galerkin moving mesh method for a class of one dimensional nonlinear PDEs. The moving mesh equations are solved using the static rezoning approach with the Lagrange dG basis functions t t t t t t t=0.5 t=0.5 t=0.5 1 1 1 0.8 0.8 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.2 0.2 0.2 Exact Exact Exact 0 0 0 DG DG DG −0.2 −0.2 −0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x x (a) (d) (g) t=1 t=1 t=1 1 1 1 0.8 0.8 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.2 0.2 0.2 Exact Exact Exact 0 0 0 DG DG DG −0.2 −0.2 −0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x x (b) (e) (h) 10 0 10 −20 −10 −10 −20 −40 −20 −30 −30 −60 −40 −40 −50 −80 −50 −60 −60 −70 −100 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 t t t (c) (f) (i) Figure 6: Schlogl ¨ equation: Traveling wave solutions and energy plots; (a)-(c) on a uniform mesh with N = 120 elements using quadratic dG basis functions, (d)-(f) on a moving mesh with N = 40 elements using linear dG basis functions, (g)-(i) on a moving mesh with N = 40 elements using quadratic dG basis functions. in the Algorithm 1. Numerical results for problems with different nature of traveling waves demonstrate the accuracy and effectiveness of the moving mesh dG method. As a future study we aim to extend the dG moving mesh to two dimensional problems in space. Energy u u Energy u u Energy u u 2 Table 2: Schlogl ¨ equation: L -errors with an adaptive moving mesh with N = 40 elements. Monitor degree t=0.001 t=0.01 t=0.25 t=0.5 t=0.75 t=1 Optimal 1 3.6e-03 1.6e-03 3.5e-01 4.3e-01 4.9e-01 5.3e-01 Arc-length 1 9.3e-03 2.1e-02 4.9e-01 5.1e-01 5.4e-01 5.3e-01 Curvature 1 3.9e-03 5.2e-03 1.2e-01 1.8e-01 2.3e-01 2.5e-01 Optimal 2 2.2e-04 3.3e-04 2.9e-03 4.4e-03 5.6e-03 6.7e-03 Arc-length 2 1.1e-03 1.5e-03 2.2e-03 3.8e-03 8.8e-03 1.4e-03 Curvature 2 3.4e-04 2.3e-04 2.9e-03 5.2e-03 7.4e-03 9.5e-03 References [1] P. Antonietti and P. Houston. A pre-processing moving mesh method for discon- tinuous Galerkin approximations of advection-diffusion-reaction problems. In- ternatioanl Journal of Numerical Analysis and Modelling, 5(4):704–728, 2008. [2] D. N. Arnold. An interior penalty finite element method with discontinuous ele- ments. SIAM J. Numer. Anal., 19:724–760, 1982. [3] Rico Buchholz, Harald Engel, Eileen Kammann, and Fredi Troltzsch. ¨ On the optimal control of the Schlogl-model. ¨ Comput. Optim. Appl., 56(1):153–185, [4] Ernst Hairer and Christian Lubich. Energy-diminishing integration of gradient systems. IMA Journal of Numerical Analysis, 2013. [5] W. Huang and R. Russell. Adaptive moving mesh methods. Applied Mathematical Sciences. Springer, 2011. [6] Weizhang Huang, Yuhe Ren, and Robert D. Russell. Moving mesh partial dif- ferential equations (MMPDES) based on the equidistribution principle. SIAM Journal on Numerical Analysis, 31(3):709–730, 1994. [7] J. M. Hyman and B. Larrouturou. Dynamic rezone methods for partial differential equations in one space dimension. Appl. Numer. Math., 5(5):435–450, 1989. [8] J.M. Hyman, Shengtai Li, and L.R. Petzold. An adaptive moving mesh method with static rezoning for partial differential equations. Computers & Mathematics with Applications, 46(1011):1511 – 1524, 2003. [9] B. Karasozen ¨ and M. Uzunca. Time-space adaptive discontinuous Galerkin method for advection-diffusion equations with non-linear reaction mechanism. International Journal on Geomathematics, 5:255–288, 2014. [10] Ruo Li and Tao Tang. Moving mesh discontinuous Galerkin method for hyper- bolic conservation laws. Journal of Scientific Computing, 27(1-3):347–363, 2006. 14 [11] B. Riviere. Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation. SIAM, 2008. [12] Ali R. Soheili, Asghar Kerayechian, and Noshin Davoodi. Adaptive numerical method for Burgers-type nonlinear equations. Applied Mathematics and Compu- tation, 219(8):3486–3495, 2012. [13] M. Uzunca, B. Karasozen, ¨ and M. Manguoglu. ˘ Adaptive discontinuous Galerkin methods for non-linear diffusion-convection-reaction equations. Computers and Chemical Engineering, 68:24–37, 2014. [14] Robert D. Russell Weizhang Huang. Adaptive Moving Mesh Methods. Spinger, http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Mathematics arXiv (Cornell University)

Moving Mesh Discontinuous Galerkin Methods for PDEs with Traveling Waves

Uzunca, Murat; Karasözen, Bülent; Küçükseyhan, Tuğba
Mathematics , Volume 2020 (1602) – Jan 31, 2016
Download PDF
15 pages

Loading next page...
 
/lp/arxiv-cornell-university/moving-mesh-discontinuous-galerkin-methods-for-pdes-with-traveling-CyGApvlv8q

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

ISSN
0096-3003
eISSN
ARCH-3343
DOI
10.1016/j.amc.2016.07.034
Publisher site
See Article on Publisher Site

Abstract

In this paper, a moving mesh discontinuous Galerkin (dG) method is developed for nonlinear partial differential equations (PDEs) with traveling wave solutions. The moving mesh strategy for one dimensional PDEs is based on the rezoning approach which decouples the solution of the PDE from the moving mesh equa- tion. We show that the dG moving mesh method is able to resolve sharp wave fronts and wave speeds accurately for the optimal, arc-length and curvature moni- tor functions. Numerical results reveal the efficiency of the proposed moving mesh dG method for solving Burgers’, Burgers’-Fisher and Schlogl(Nagumo) ¨ equations. 1 Introduction The discontinuous Galerkin (dG) method is one of the most powerful discretization techniques for solving partial differential equations (PDEs) [2, 11], especially for con- vection dominated problems, exhibiting localized phenomena like sharp traveling wave fronts, internal and boundary layers [9, 13]. The dG method has been applied for this kind of singularly perturbed linear and nonlinear PDEs extensively using h-adaptive (refinement and coarsening in space), p-adaptive (enrichment of the local polynomial degree), hp-adaptive and space-time adaptive methods in the last two decades. An- other approach to deal with these kind of problems, is the r-method or moving mesh method. In the moving mesh method the grid points are relocated in the regions where the solution shows rapid variation, while keeping the number of the nodes fixed. The dG discretization is very flexible, since there is no continuity requirement between the inter-element boundaries, which makes it suitable as a moving mesh method on ir- regular meshes. Most of the studies with moving mesh methods are limited to finite difference and continuous finite element discretization [14]. There are only few pub- lications dealing with dG moving mesh method. They include the interior penalty dG method for preprocessing the solutions of steady state diffusion-convection-reaction equations [1], and the local dG moving mesh method for hyperbolic conservation laws [10]. In this paper we develop an adaptive dG moving mesh method for one dimensional arXiv:1602.05225v1 [math.NA] 31 Jan 2016 semi-linear differential equations with traveling wave solutions of the form u = e u f (u; u ); (x; t)2 W(t ; T ] (1a) t xx x 0 f u(x ; t) = u ; u(x ; t) = u ; t 2 (t ; T ] (1b) L L R R 0 f u(x; t ) = u ; x2 W; (1c) 0 0 where W = [x ; x ] R, t and T are the initial and final time instances, respectively, L R 0 and e denotes the diffusion coefficient. The model equation (1) becomes the Burgers equation with f (u; u ) = uu [12] , Burgers’-Fisher equation with f (u; u ) = a uu + x x x x b u(u 1) [12] and the Schlogl ¨ or Nagumo equation with f (u; u ) = u(1 u)(1b)=d [14]. A moving mesh method has three main components; the discretization of the phys- ical PDE, mesh strategy using monitor functions and discretization of the mesh equa- tion. The discretization of physical PDE is either coupled with the moving mesh equa- tion or separated. In the quasi-Lagrangian approach, a large system of the discretized PDE and moving mesh equation are solved simultaneously by the standard ordinary differential equation (ODE) solvers. Instead, we use the rezoning approach by solving alternately the PDE and mesh equation, which allows more flexibility; mesh genera- tion can be coded separately and embedded in the solution of the PDE. Since the mesh is updated at each time step, the physical PDE has to be discretized at the next time step on the new mesh. We use the static rezoning approach with the same number of points at each time step [7] in contrast to the dynamic rezoning [8] where the number of mesh points is changed at every time step. Therefore in the static rezoning approach the solutions from old to new mesh have to be interpolated. The paper is organized as follows. In the next section we describe briefly the dG method for the 1D model problem (1) on a uniform fixed mesh. Moving mesh adaption strategy and the adaptive moving mesh dG algorithm is presented in Section 3. Nu- merical results are given in Section 4 to demonstrate the effectiveness of the proposed method. 2 Discretization of the problem on a fixed mesh Before giving the moving mesh strategy in Section 3, in this section we describe the dG discretization of the model problem (1) on a fixed uniform mesh T : x = x + nh; n = 0; 1;:::; N ; (2) h n L I consisting of N elements (sub-intervals) I = [x ; x ], n = 1; 2;:::; N , and with the I n n1 n I fixed mesh size h = (x x )=N . R L I 2.1 Space discretization by discontinuous Galerkin method We use for the space discretization of the model problem (1) on a fixed mesh (2) the symmetric interior penalty Galerkin (SIPG) method [2, 11] which is a member of the 2 family of dG methods. The dG methods use the space of piecewise discontinuous polynomials of degree at most k: V =fv : vj 2 P (I ) ; 8n = 1;:::; N g; h I k n I where P (I ) is the space of polynomials of degree at most k on an interval I . Since the k n n functions in V are discontinuous at the inter-element nodes, we define the jump and average of a piecewise function v at the endpoints of I , n = 1;:::; N 1, respectively, n I as depicted in Figure 1, + + [v(x )] = v(x ) v(x ) ; fv(x )g = (v(x )+ v(x )); (3) n n n n n n with v(x ) = lim v(x) ; v(x ) = lim v(x): (4) n n x7!x x7!x n n On the boundary nodes, the jump and average are defined as + + + [v(x )] =v(x ); fv(x )g = v(x ); [v(x )] = v(x ); fv(x )g = v(x ): (5) 0 0 N N 0 0 I N I N I I v(x ) v| [v(x )] v| n+1 v(x ) x x x x x x N −1 N 0 1 n−1 n+1 I I Figure 1: Jump and limit terms of a piecewise discontinuous function v(x). The SIPG scheme is constructed by multiplying the continuous (the solution u is sufficiently smooth at the end points of I ) equation (1) by a test function v 2 V and n h integrating by parts on each element I , n = 1;:::; N : n I Z Z Z x x x n n n u vdx+ e u v dxe u (x )v(x )+e u (x )v(x )+ f (u; u )vdx = 0: t x x x n x n1 x n1 x x x n1 n1 n1 By adding all N equations, and using the definition of the jumps (3) and (5), we obtain Z Z Z N N I x x x I n n n u vdx+ e u v dx+ f (u; u )vdx [e u (x )v(x )] = 0: t x x x x n n å å x x x n1 n1 n1 n=1 n=0 One can verify that for 1 n N 1 [e u (x )v(x )] =fe u (x )g[v(x )]+[e u (x )]fv(x )g: (6) x n n x n n x n n Using the identity (6) and the fact that [e u (x )] = 0 for all 1  n  N 1 (u was x n I sufficiently smooth at the end points of I ), we obtain Z Z Z N N I x x x I n n n u vdx+ e u v dx+ f (u; u )vdx fe u (x )g[v(x )] = 0: (7) t x x x x n n å å x x x n1 n1 n1 n=1 n=0 | {z } Additionally, we have [u(x )] = 0 for all 1 n N 1. Then, adding the penalizing n I terms and the terms on the boundary nodes, n =f0; N g, to both sides of (7) by keeping them unknown on the left hand side and imposing the boundary conditions u and u L R on the right hand side, leads to the SIPG formulation: Z Z x x R R u vdx+ a(u; v)+ f (u; u )vdx = l(v): (8) t x x x L L In (8), a(u; v) and l(v) denote the symmetric bilinear form and the linear right hand side of the SIPG scheme x I a(u; v) = e u v dx+ fe u (x )g[v(x )]fe v (x )g[u(x )]+ [u(x )][v(x )] ; x x x n n x n n n n n=0 s s l (v ) = u e v (x ) v(x ) + u v(x )e v (x ) ; h h L x 0 0 R N x N I I h h where s is the penalty parameter which should be sufficiently large [11] in order to ensure the coercivity of the bilinear form. Hence, SIPG semi-discrete form of (1) reads as: a.e. t 2 (t ; T ], for all v 2 V , find u := u (x; t)2 V such that 0 f h h h h h Z Z x x R R u (x; t )v dx = u v dx; (9a) h 0 h 0 h x x L L Z Z x x R R (¶ u )v dx+ a(u ; v )+ f (u ; u )v dx = l(v ): (9b) t h h h h h h;x h h x x L L 2.2 Full discretization Let ff g denote the basis functions spanning the dG finite elements space V for i = 1;:::; N and n = 1;:::; N , where N stands for the local dimension depending on the k I k polynomial order k and is given by N = k + 1 in 1D. The local nature of dG methods leads to the basis functions and the approximate solution of the form n I k f (x) ; x2 I n n n f (x) = ; u (t) = u (t)f (x); (10) i h å å i i 0 ; x2 = I n=1 i=1 wherefu (t)g are the time-dependent unknown coefficients. We substitute the second identity in (10) into the system (9b) and we choose v = f for i = 1;:::; N and n = h k 1;:::; N , which leads to the N = N  N dimensional non-linear system of equations I k I of (9b) in matrix-vector form Mu + Su + h(u) d = 0; (11) NN NN where M2 R is the usual mass matrix and S2 R is the stiffness matrix related N N to the bilinear form a(u ; v ). The vectors h(u) 2 R and d 2 R are the non-linear h h vector of unknown coefficients u corresponding to the integral of the nonlinear term in (9b) and the right hand side vector related to the linear form l(v ), respectively. The initial vector u(0) is found by using the equation (9a) and the second identity in (10). 4 We solve the fully discrete system of (1), by applying the backward Euler scheme to the semi-discrete system (11). Let t < t < ::: < t = T be the uniform partition 0 1 J f of the time interval [t ; T ] into J time-steps (t ; t ], j = 1;:::; J, with the step size 0 f j1 j Dt = (T t )=J. Let us denote the approximate coefficient vector u(t) of (11) at the f 0 time t = t by u . Then, the fully discrete formulation of the model (1) is given as: for j N all j = 1;:::; J, find u 2 R such that j j1 u u j j j M + Su + h(u ) d = 0; (12) Dt which is solved by Newton’s method. 3 Adaptive moving mesh method In an adaptive moving mesh method the spatial mesh T in (2) varies with time T (t) : x = x(x ; t); n = 0; 1;:::; N ; (13) h n n I consisting of N elements I = [x ; x ], n = 1; 2;:::; N , of the mesh size h = x I n n1 n I n n x . In (13) the nodes x belong to the fixed and uniform mesh n1 n T : x = ; n = 0; 1;:::; N (14) n I on the auxiliary unit interval W = [0; 1], together with the boundary conditions x(0; t) = x and x(1; t) = x . Thus, in the adaptive moving mesh method, the mesh is adjusted as L R the time progresses in such a way that the mesh sizes h get smaller in the sub-intervals where the error is large, while the mesh sizes h are made coarser in the remaining part of the interval. The error indicator is chosen to relocate the mesh points where the solution shows large variations, based on the equidistribution principle, where a mesh T (t) with the mesh points x = x < x < ::: < x < x = x is determined by h L 0 1 N 1 N R I I satisfying the following relation x N 1 I Z Z r(x; t)dx = = r(x; t)dx: (15) x x 0 N 1 In this way a continuous function r(x; t) > 0 on the interval [x ; x ] can be distributed L R among the sub-intervals I = [x ; x ], n = 1;:::; N . In (15), the function r(x; t) is n n1 n I called the mesh density function, or monitor function, choice of which stands as the key point for an adaptive moving mesh method. The most popular choices for the monitor functions r := r(x; t) are [5] optimal 1=3 r = 1+ ju j ; (16) xx 5  arc-length 1=2 r = 1+ju j ; (17) curvature 1=4 r = 1+ju j ; (18) xx with the intensity parameter ( ) 2=3 a = max 1; ju j dx : (19) xx x x R L L Finding a proper mesh T (t) using the equidistribution condition (15) results in a system of so-called moving mesh partial differential equation (MMPDE) [6, 5] ¶ x 1 ¶ ¶ x = r ; (x; t)2 W (t ; T ]; (20a) c 0 f ¶ t tr ¶x ¶x x(0; t) = x ; x(1; t) = x ; t 2 (t ; T ]: (20b) L R 0 f The system (20) is solved through the central finite differences dx 1 r +r r +r n n+1 n n n1 = (x + x ) (x + x ) (21) n+1 n n n1 dt tr Dx 2 2 for n = 1; 2;:::; N 1, where the spatial nodes x 2 T (t) are the unknown solutions I n h of the nodes x 2T . The relaxation parameter t is specified by the user for adjusting the response time of mesh movement according to the changes of the density function r(x; t). The functions r , n = 0; 1;:::; N , are computed through the discrete form of n I the monitor function r . For instance, we have for the optimal monitor function (16) 1=3 r = 1+ ju j ; n = 0;:::; N ; (22) n xx;n I where the terms u are computed by the central difference approximations using xx;n the solutions at the mesh nodes x . At each time step, further, the discrete monitor functions r are smoothed by weighted averaging [5, Section 1.2]. There are different approaches [14] of the adaptive moving mesh method for solv- ing the fully discrete system (12), called the physical PDE, and the MMPDE (21). One common approach is solving both systems simultaneously using a quasi-Lagrange ap- proach. In this approach, the time derivative term requires a special attention since the mesh is assumed to move in a continuous manner by the time progresses. In such cases, there occur an extra convective term which may cause difficulties. Additionally, the solution of the both systems simultaneously needs the coupling of the systems, as a result the dimension of the system to be solved increases. Another choice, which we use in this paper, is the alternate solution using a rezoning approach. In this case, the physical PDE and the MMPDE are separated from each other, which allows more flex- ibility as the mesh generation can be coded separately and embedded in the solution 6 of the physical PDE. In the static rezoning approach, the change in the spatial mesh is derived in a discrete form similar to the solution [7, 8, 14]. In each time step, first the mesh adaptation is handled using the solution on the old mesh, then the solution is obtained by solving the physical PDE on the newly generated mesh. In the case of dG method, the unknown coefficient vectors have to satisfy u u(t ) in the physical PDE (12) with the discrete mesh density function (22) in the MMPDE (21). In general, it is difficult to use the coefficient vectors u directly to compute the discrete mesh density functions r , unless we use the Lagrange nodal basis. Let us choose the dG basis functions f as the Lagrange nodal basis functions, i = 1;:::; N , n = 1;:::; N . Then, recalling the definitions (4) of the limit terms at a node x , we have the relations n + n+1 u (x ; t) = u (t) ; u (x ; t) = u (t); n = 1;:::; N 1: (23) h h I n N n The relations (23) originate from the fact that dG methods use discontinuous basis functions at the inter-element nodes. As a result, a dG approximation u (x; t) has two traces at an inter-element node x from the neighboring sub-intervals I and I , as n n n+1 shown in Figure 1, which are not the same in general. Using the relations (23) and recalling again the average definition in (3), we can accept the value of the approximate solution u (x ; t) as h n u (x ; t) : =fu (x ; t)g = (u (x ; t)+ u (x ; t)) h n h n h h n n (24) n n+1 = (u (t)+u (t)); n = 1;:::; N 1: In this way, the solutions of the physical PDE (12) will be consistent with the MM- PDE (21), where the computation of the discrete mesh density functions r require the discrete approximations u (x; t) at the inter-element nodes x , n = 0;:::; N . h n I Algorithm 1 Moving mesh algorithm j1 j1 j1 Given the current spatial mesh T (t ) : x < x < ::: < x , co- h j1 0 1 N j1 efficient vector u , the parameter t and time-step size Dt , do for j = 1;:::; J, 1: Compute a temporary coefficient vector u ˜ by solving the physical PDE (12) on the current mesh T (t ), h j1 2: According to (24), calculate the consistent values of u ˜ (x ; t ) using u , h n j 3: Compute the discrete mesh density functions r ˜ using the values u ˜ (x ; t ), n n j j j j 4: Find the mesh T (t ) : x < x < ::: < x by solving the MMPDE (21) for the h j 0 1 N discrete mesh density functions r , j1 5: Interpolate the coefficient vector u to be used in the new mesh T (t ), h j 6: Compute the coefficient vector u by solving the physical PDE (12) on the new mesh T (t ), h j 7: Go to next time step. The general algorithm is summarized in Algorithm 1. We start with an initial mesh T (t ) (possibly a uniform mesh). Then on an arbitrary time-step (t ; t ], we, firstly, h 0 j1 j 7 j solve the physical PDE (12) on the mesh T (t ) for an auxiliary solution u ˜ . After, h j1 following the relation (24), we use the solution u ˜ in the MMPDE (21) to obtain the new mesh T (t ). Finally, we solve the physical PDE (12) on the mesh T (t ) for h j h j the solution u , and we proceed to the next time-step. On the other hand, the known j1 solution vector u will be no more consistent with the new mesh T (t ). This is h j a natural consequence of the alternate solution with rezoning approach. To make the j1 known solution u consistent with the updated spatial mesh T (t ), we interpolate it between the meshes T (t ) and T (t ). The interpolation procedure is as follows: let j1 j h h j1 j1 j1 j j j T (t ) : x < x < ::: < x be the current mesh andT (t ) : x < x < ::: < x j1 j h h 0 1 N 0 1 N I I j1 n ˆ ˆ denote the updated mesh. For the sake of simplicity, let also u := u with entries u being the coefficient vector defined on the current meshT (t ), and u with entries u h j1 denotes the interpolated coefficient vector of u into the new mesh T (t ), i = 1;:::; N , h j k j1 j1 n = 1;:::; N . The local nature of dG methods leads for any x2 I = [x ; x ] to I s s s1 s s u (x; t ) = u f (x); s = 1;:::; N : (25) h j1 I å i i i=1 On the other hand, using the Lagrange basis functions, on an arbitrary element we j j obtain I = [x ; x ] on the new mesh T (t ), uniformly distributed N nodal degrees n n h j k n1 of freedoms x 2 I such that n;d j j j x = x +(d 1)t; u (x ; t ) u ; (26) h j1 n;d n1 n;d j j1 for d = 1;:::; N , and with t = (x x )=(N 1). In other words, the entries u k n n k of the interpolated coefficient vector u are the approximate function values u (x; t ) h j1 at the nodal degrees of freedoms x , n = 1;:::; N , d = 1;:::; N . For any nodal I k n;d degrees of freedom x on the new mesh T (t ), we have to determine the intervals h j n;d j1 j1 j I = [x ; x ] on the current mesh T (t ) such that x 2 I . Then, using the s s h j1 s s1 n;d expansion (25), we will be able to obtain the entries of u as u = u (x ; t ). h j1 i n;i 4 Numerical results In this section we present several numerical examples demonstrating the effectiveness of the adaptive moving mesh dG method. In all of the examples, the traveling wave solutions are computed by the optimal mesh density function (16), but the correspond- ing mesh trajectories are given for the optimal (16), arc-length (17) and curvature (18) mesh density functions. 4.1 Burgers’ equation The first test example is the Burgers’ equation [14, 12] ¶ 1 u = e u u t xx ¶ x 2 8 with homogeneous Dirichlet boundary conditions on the space-time domain (x; t) 2 [0; 1] [0; 1] with the diffusion constant e = 10 . The initial condition is taken as u(x; 0) = sin(2p x) + 0:5 sin(p x), and linear dG basis functions are used. We set the time-step size Dt = 0:005, and we choose the relaxation parameter as t = 10 . Moving mesh solutions in Figure 3 are capable of resolving the sharp gradients of the moving fronts in contrast to the oscillatory solutions on the fixed mesh in Figure 2. In addition, in Figure 3, mesh trajectories for the curvature and arc-length monitor functions show large variations with respect to time, whereas for the optimal monitor function the mesh trajectories are smooth. 1.5 0.5 t=0 −0.5 t=0.2 t=0.4 −1 t=0.6 t=0.8 t=1 −1.5 0 0.2 0.4 0.6 0.8 1 Figure 2: Burgers’: Solutions at t = 0, 0.2, 0.4, 0.6, 0.8, 1 on the uniform fixed mesh with N = 120 elements. 4.2 Burgers’-Fisher equation Consider the Burgers’-Fisher equation [12] ¶ 1 u = u a u +b u(u 1); t xx ¶ x 2 with the exact solution 1 a u(x; t) = 1 tanh (x ct) ; 2 4 in the space-time domain (x; t)2 [1; 0](0:2; 0]. The parameters are a = 24, c = 8 and b = (2a c a )=4. We use quadratic dG basis functions and N = 40 spatial elements. The time step-size is taken as Dt = 0:001. In Figure 4, we give the solutions of Burgers’-Fisher equation at different time instances for the optimal mesh density function, and the moving mesh trajectories ob- tained by different monitor functions. In Table 1 the L -errors between the exact and numerical solutions are tabulated at different times. The results are very close to those in [12] computed with the same settings. u 1 0.9 0.9 0.8 0.8 0.7 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (a) Optimal (b) Arc-Length 1.5 0.9 0.8 0.7 0.6 0.5 0.5 0.4 0.3 t=0 t=0.2 0.2 −0.5 t=0.4 t=0.6 0.1 t=0.8 t=1 −1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x (c) Curvature (d) Solutions Figure 3: Burgers’ equation: (a-c) Moving mesh trajectories with N = 40 elements, (d) solutions at t = 0, 0.2, 0.4, 0.6, 0.8, 1 for the optimal monitor function. Table 1: Burgers’-Fisher equation: L -errors Monitor N t=-0.1 t=-0.05 t=-0.04 t=-0.035 t=-0.03 Optimal 40 4.6e-03 8.4e-03 1.1e-02 1.1e-02 1.3e-02 Arc-Length 40 2.4e-03 4.2e-03 5.2e-03 5.8e-03 6.5e-03 Curvature 40 2.4e-03 4.1e-03 5.1e-03 5.7e-03 6.4e-03 4.3 Schlogl ¨ equation The final example is Schlogl ¨ (Nagumo) equation [3, 14] u = e u f (u); t xx with the Ginzburg-Landau free energy E (u) = jÑuj + F(u) dx (27) t t t 0 −0.02 −0.02 −0.04 −0.04 −0.06 −0.06 −0.08 −0.08 −0.1 −0.1 −0.12 −0.12 −0.14 −0.14 −0.16 −0.16 −0.18 −0.18 −0.2 −0.2 −1 −0.8 −0.6 −0.4 −0.2 0 −1 −0.8 −0.6 −0.4 −0.2 0 x x (a) Optimal (b) Arc-Length 0 1 t=−0.1 −0.02 0.9 t=−0.05 −0.04 0.8 −0.06 0.7 t=−0.04 −0.08 0.6 t=−0.035 −0.1 0.5 t=−0.03 −0.12 0.4 −0.14 0.3 −0.16 0.2 −0.18 0.1 −0.2 −1 −0.8 −0.6 −0.4 −0.2 0 −1 −0.8 −0.6 −0.4 −0.2 0 x x (c) Curvature (d) Solutions Figure 4: Burgers’-Fisher equation: (a-c) Moving mesh trajectories with N = 40 el- ements, (d) solutions at t = -0.1, -0.05, -0.04, -0.035, -0.03 and corresponding exact solutions by solid lines. 2 2 with quartic potential function F(u) = u (3u 4(1+b)u+ 6b) and cubic bi-stable 12d nonlinearity f (u) = u(u 1)(ub) satisfying f (u) = F (u). We consider Schlogl ¨ equation in the space-time domain (x; t)2 (0; 1] [0; 1] with 3 3 the constant parameter values e = 10 , d = 10 and b = 0. The Dirichlet boundary conditions and the initial condition are taken according to the exact solution 1 x ct u(x; t) = 1 tanh p ; 8ed where c = e=2d is the speed of the traveling wave. The solutions are computed with linear and quadratic dG basis functions, and with the time-step size Dt = 0:001. As the number of spatial elements, we take N = 120 and N = 40 to construct a uniform and I I moving mesh, respectively. t t u Figure 5 shows that the steep wave fronts are captured well by the adaptive moving mesh method with linear and quadratic dG basis functions. On the other hand, Figure 6 shows that the numerical solutions with quadratic dG basis functions give the correct wave speed, whereas for the linear case the numerical solutions move faster than the exact solutions. In Table 2, we list the L -errors between the numerical and exact so- lutions at different time instances for different types of monitor functions. The correct wave speed is not captured by linear dG basis functions. Therefore the corresponding L -errors are larger than the quadratic case. Due the discontinuous nature of the dG discretization, higher order basis functions can be implemented in a more flexible way in contrast to continuous finite elements, providing continuity requirement between the inter-element boundaries. The free energy of Schlogl ¨ equation (27) decreases monoton- ically in time. The energy decreasing property is also captured numerically on uniform and moving meshes, as shown in Figure 6 for the conditionally energy stable backward Euler method [4]. 1 1 1 0.9 0.9 0.9 0.8 0.8 0.8 0.7 0.7 0.7 0.6 0.6 0.6 0.5 0.5 0.5 0.4 0.4 0.4 0.3 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.1 0 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x x (a) Optimal (b) Arc-length (c) Curvature 1 1 0.9 0.9 0.9 0.8 0.8 0.8 0.7 0.7 0.7 0.6 0.6 0.6 0.5 0.5 0.5 0.4 0.4 0.4 0.3 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.1 0 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x x (d) Optimal (e) Arc-length (f) Curvature Figure 5: Schlogl ¨ equation: Moving mesh trajectories with N = 40 elements; (a)-(c) linear dG basis functions, (d)-(f) quadratic dG basis functions. 5 Conclusions In this paper we have developed an adaptive discontinuous Galerkin moving mesh method for a class of one dimensional nonlinear PDEs. The moving mesh equations are solved using the static rezoning approach with the Lagrange dG basis functions t t t t t t t=0.5 t=0.5 t=0.5 1 1 1 0.8 0.8 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.2 0.2 0.2 Exact Exact Exact 0 0 0 DG DG DG −0.2 −0.2 −0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x x (a) (d) (g) t=1 t=1 t=1 1 1 1 0.8 0.8 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.2 0.2 0.2 Exact Exact Exact 0 0 0 DG DG DG −0.2 −0.2 −0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x x (b) (e) (h) 10 0 10 −20 −10 −10 −20 −40 −20 −30 −30 −60 −40 −40 −50 −80 −50 −60 −60 −70 −100 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 t t t (c) (f) (i) Figure 6: Schlogl ¨ equation: Traveling wave solutions and energy plots; (a)-(c) on a uniform mesh with N = 120 elements using quadratic dG basis functions, (d)-(f) on a moving mesh with N = 40 elements using linear dG basis functions, (g)-(i) on a moving mesh with N = 40 elements using quadratic dG basis functions. in the Algorithm 1. Numerical results for problems with different nature of traveling waves demonstrate the accuracy and effectiveness of the moving mesh dG method. As a future study we aim to extend the dG moving mesh to two dimensional problems in space. Energy u u Energy u u Energy u u 2 Table 2: Schlogl ¨ equation: L -errors with an adaptive moving mesh with N = 40 elements. Monitor degree t=0.001 t=0.01 t=0.25 t=0.5 t=0.75 t=1 Optimal 1 3.6e-03 1.6e-03 3.5e-01 4.3e-01 4.9e-01 5.3e-01 Arc-length 1 9.3e-03 2.1e-02 4.9e-01 5.1e-01 5.4e-01 5.3e-01 Curvature 1 3.9e-03 5.2e-03 1.2e-01 1.8e-01 2.3e-01 2.5e-01 Optimal 2 2.2e-04 3.3e-04 2.9e-03 4.4e-03 5.6e-03 6.7e-03 Arc-length 2 1.1e-03 1.5e-03 2.2e-03 3.8e-03 8.8e-03 1.4e-03 Curvature 2 3.4e-04 2.3e-04 2.9e-03 5.2e-03 7.4e-03 9.5e-03 References [1] P. Antonietti and P. Houston. A pre-processing moving mesh method for discon- tinuous Galerkin approximations of advection-diffusion-reaction problems. In- ternatioanl Journal of Numerical Analysis and Modelling, 5(4):704–728, 2008. [2] D. N. Arnold. An interior penalty finite element method with discontinuous ele- ments. SIAM J. Numer. Anal., 19:724–760, 1982. [3] Rico Buchholz, Harald Engel, Eileen Kammann, and Fredi Troltzsch. ¨ On the optimal control of the Schlogl-model. ¨ Comput. Optim. Appl., 56(1):153–185, [4] Ernst Hairer and Christian Lubich. Energy-diminishing integration of gradient systems. IMA Journal of Numerical Analysis, 2013. [5] W. Huang and R. Russell. Adaptive moving mesh methods. Applied Mathematical Sciences. Springer, 2011. [6] Weizhang Huang, Yuhe Ren, and Robert D. Russell. Moving mesh partial dif- ferential equations (MMPDES) based on the equidistribution principle. SIAM Journal on Numerical Analysis, 31(3):709–730, 1994. [7] J. M. Hyman and B. Larrouturou. Dynamic rezone methods for partial differential equations in one space dimension. Appl. Numer. Math., 5(5):435–450, 1989. [8] J.M. Hyman, Shengtai Li, and L.R. Petzold. An adaptive moving mesh method with static rezoning for partial differential equations. Computers & Mathematics with Applications, 46(1011):1511 – 1524, 2003. [9] B. Karasozen ¨ and M. Uzunca. Time-space adaptive discontinuous Galerkin method for advection-diffusion equations with non-linear reaction mechanism. International Journal on Geomathematics, 5:255–288, 2014. [10] Ruo Li and Tao Tang. Moving mesh discontinuous Galerkin method for hyper- bolic conservation laws. Journal of Scientific Computing, 27(1-3):347–363, 2006. 14 [11] B. Riviere. Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation. SIAM, 2008. [12] Ali R. Soheili, Asghar Kerayechian, and Noshin Davoodi. Adaptive numerical method for Burgers-type nonlinear equations. Applied Mathematics and Compu- tation, 219(8):3486–3495, 2012. [13] M. Uzunca, B. Karasozen, ¨ and M. Manguoglu. ˘ Adaptive discontinuous Galerkin methods for non-linear diffusion-convection-reaction equations. Computers and Chemical Engineering, 68:24–37, 2014. [14] Robert D. Russell Weizhang Huang. Adaptive Moving Mesh Methods. Spinger,

Journal

Mathematics – arXiv (Cornell University)

Published: Jan 31, 2016

References

  • An interior penalty finite element method with discontinuous elements
    Arnold