MA557/MA578/CS557Lecture 19 Spring 2003 Prof. Tim Warburton email@example.com
DG Scheme for the Scalar ADE(Lax-Friedrichs flux – equivalent to upwind)
DG Derivative Operator • We are going to introduce a DG derivative operator to simplify the scheme definition. • The linear operator Dtilde is such that the following holds for all intervals Ij • With the choice of penalty terms tauL and tauR to be determined.
In Operator Notation Advectionterm(Lax-Friedrichs flux ~ upwind fluxfor scalar case) Diffusionterm
General Dtilde Operator • We previously defined Dtilde for the special case where the test function is compactly supported on the j’th cell. • We generalize by:
Skew-symmetry of Dtilde1,1 [integrate by parts]
Skew-symmetry • Note:the negative sign when the operator moves. • Trivially:
Recap • We already proved that the DG scheme, with Lax-Friedrichs fluxes, for advection (D=0) is stable in the sense that : • Adding the diffusion term (dropping boundary terms) gives: • In terms of norms:
Stability!! • Assuming D>=0 (reasonable since particles do not jump randomly with negative rates)
D operator • If we move to the Legendre basis we can write down a discrete description of Dtilde:
Implementation of Dtilde function [dfdx] = LEGdgderiv(f, nodex, tauL, tauR, D, F, G, H, J) [p,N] = size(f); % p = maximum polynomial order p=p-1; % N = number of nodes N=N+1; % space for derivative dfdx = zeros(p+1,N-1); % cells 2 to N-2 ids = (2:N-2); dfdx(:,ids) = (D*f(:,ids) + tauL*F*f(:,ids) + tauL*G*f(:,ids-1) + tauR*H*f(:,ids) + tauR*J*f(:,ids+1)); % cell 1 (periodicity assumed) dfdx(:,1) = (D*f(:,1) + tauL*F*f(:,1) + tauL*G*f(:,N-1) + tauR*H*f(:,1) + tauR*J*f(:,2)); % cell N-1 (periodicity assumed) dfdx(:,N-1) = (D*f(:,N-1) + tauL*F*f(:,N-1) + tauL*G*f(:,N-2) + tauR*H*f(:,N-1) + tauR*J*f(:,1)); % apply chain rule for physical cell width dx = nodex(2:N)-nodex(1:N-1); coeff = ones(p+1,1)*(2./dx); dfdx = dfdx.*coeff;
In Action (advection eqn) % TIME STEPPING tauL = (1/u)*(u+abs(u))/2; tauR = (1/u)*(u-abs(u))/2; for tstep = 1:Ntsteps sigma = rho; for rkstage = 3:-1:1 dsigmadx = LEGdgderiv(sigma, nodex, tauL, tauR, D, F, G, H, J); sigma = rho + (dt/rkstage)*(-u*dsigmadx); end rho = sigma; if ( ~mod(tstep, 100) ) q = V*rho; % plot(x(:),abs(q(:)-exp(-(x(:)-dt*tstep).^2))); plot(x,q); pause(0.1); tstep*dt end end
Discrete Scheme • Step 1: q = x DG derivative of C • Step 2: time rate of change of C Advectionterm(Lax-Friedrichs flux ~ upwind fluxfor scalar case) Diffusionterm
Dropping q • We can do away with q: • This should make the following observation clear. The dt dependence for stability is now: • The first term is due to the spectral radius of the advection operator. • The second term is due to the spectral radius of the diffusion operator.
In Action (ADE) % TIME STEPPING for tstep = 1:Ntsteps sigma = rho; for rkstage = p:-1:1 % advection: upwind bias on flux terms dsigmadx = LEGdgderiv(sigma, nodex, tauL, tauR, D, F, G, H, J); % diffusion: centered fluxes dsdx = LEGdgderiv(sigma, nodex, 1, 1, D, F, G, H, J); d2sdx2 = LEGdgderiv(dsdx, nodex, 1, 1, D, F, G, H, J); sigma = rho + (dt/rkstage)*(-u*dsigmadx+Dcoeff*d2sdx2); end rho = sigma; if ( ~mod(tstep, 400) ) q = V*rho; plot(x,q); axis([xmin xmax 0 1]) pause(0.1); end end
Downloads • You can download the Dtilde based scripts from: • http://www.math.unm.edu/~timwar/MA578S03/MatlabScripts/LEGadvectionLaxF.m • http://www.math.unm.edu/~timwar/MA578S03/MatlabScripts/LEGade.m • You will need the Dtilde routine (named LEGdgderiv.m): • http://www.math.unm.edu/~timwar/MA578S03/MatlabScripts/LEGdgderiv.m
Legendre Scheme To Lagrange Scheme(modal to nodal) • So far we have used Legendre expansions to represent the solution in each cell. • This representation causes some problems if we wish to evaluate non-linear function of the solution. • For instance solving the inviscid Burger’s equation:requires the evaluation of u2
Transforming the Modal Advection DG Scheme To a Nodal Advection DG Scheme • Recall that given a Legendre expansion we may evaluate it at a set of points by: • We substitute: into the scheme:
Nodal Scheme • Multiplying both sides on the left with V • This clearly shows that the nodal scheme is a change of basis away from the modal scheme.
Summary of Nodal Scheme Note: we do not need to transform the initial condition or solution for visualization…
Next Lecture • The lecture on Friday 03/07 is cancelled in favour of the Analysis Conference • The lecture on Monday 03/10 will be the first time we generalize the DG to two spatial dimensions. • Do not miss the Monday lecture • Run the LEGade.m script (you will need the LEGdgderiv.m script as well) • How to use nodal methods in general.