300 likes | 411 Views
MA557/MA578/CS557 Lecture 32. Spring 2003 Prof. Tim Warburton timwar@math.unm.edu. Heat Equation. Recall Lecture 17 where we derived the advection-diffusion equation in 1D. The same type of derivation for the diffusion process of a concentrate C yields in 2D:.
E N D
MA557/MA578/CS557Lecture 32 Spring 2003 Prof. Tim Warburton timwar@math.unm.edu
Heat Equation • Recall Lecture 17 where we derived the advection-diffusion equation in 1D. • The same type of derivation for the diffusion process of a concentrate C yields in 2D:
Discretization of the Laplacian • The LDG, semidiscrete, equation satisfies: • Zero Dirichlet boundary conditions (concentration pinned to zero at boundary): • Zero Neumann boundary conditions (insulated boundaries):
Heat Equation • Our first approach will be an explicit time integration of the equation (using zero Neumann bc):
Explicit Numerical Heat Equation • The trouble is that the discrete diffusion matrix has a large spectral radius: • We do know that (using the integration by parts formula for the DG operator) that the operator is non-positive. • We also need to fit the eigenspectrum in the appropriate stability region:
We also know that all the eigenvalues are non-positive real numbers and dt*rho(A) must fit in the stability region..
Comment • We have shown that using an explicit time integration scheme may require very small time steps… • The root cause of this is that information diffuses through the data domain very quickly.. • We will investigate some schemes which solve globally coupled systems.
Backwards Euler Time Integration • Suppose we consider the backwards Euler scheme: • Iterator:
Backwards Euler Time Integration • Suppose we consider the backwards Euler scheme: • Thus the only conditions for stability are: • However, accuracy depends on dt. The discretization of the time derivative has a first order truncation error:
Crank-Nicholson • We replace the right hand side of the discrete scheme with:
Diagonalize • If we apply a change of basis so that the inv(M)*L is diagonal then we decouple the iterator into a sequence of independent scalar iterators:
Quick Result • We can bound the iteration multiplier: • As long as both the eigenvalues and dt are non-negative then the multiplier is bounded by one and the scheme is stable. • We have previously shown that L is a non-negative operator.
Crank-Nicholson Accuracy • We can briefly compute the temporal accuracy:
ESDIRK4 • A new implicit RK scheme with a lower triangle Butcher tableau has been proposed by Carpenter et al. http://fun3d.larc.nasa.gov/papers/carpenter.pdf • The linear system to invert at each sub stage is the same. • The RK scheme comes with an embedded scheme to provide s’th order and (s-1)’th order time approximations:
4th Order ESDIRK4 Implicit RK Schemeinitializer + 5 stages • Iterator: (Ctilde are intermediate values of C, Chat6 is the embedded 3rd order approximation of Cn+1)
4th Order ESDIRK4 Implicit RK Scheme • Iterator: (Ctilde are intermediate values of C, Chat5 is the embedded 3rd order approximation of Cn+1)
4th Order ESDIRK4 Implicit RK Scheme • If L is a linear operator then:
Coefficients Gamma = 1/4 a21 = 1/4 a31 = 8611/62500 a32 = -1743/31250 a41 = 5012029/34652500 a42 = -654441/2922500 a43 = 174375/388108 a51 = 15267082809/155376265600 a52 = -71443401/120774400 a53 = 730878875/902184768 a54 = 2285395/8070912 b1 = 82889/524892 b2 = 0 b3 = 15625/83664 b4 = 69875/102672 b5 = -2260/8211 bhat1 =4586570599/29645900160 bhat2 = 0 bhat3 = 178811875/945068544 bhat4 = 814220225/1159782912 bhat5 = -3700637/11593932 bhat6 = 61727/225920
Implementation Details • As they are computed we are required to store Ctilde1, Ctilde2, …, Ctilde6 until the end of the 6th ESDIRK stage. • Each stage 2,3,4,5,6 requires the solution of the same linear system, with different data on the right hand side. • Chat6 may be computed without solving a linear system. • The systems can be solved using an iterative (say conjugate gradient) method. i.e. the L matrix does not need to be computed and stored. • We need to be able to compute (M+dt*gamma*L)*x for a given vector x. • ||Ctilde6-Chat6|| will give an estimate of the difference between a 3rd order and the 4th order approximation in time. This can be used to estimate the error made in this time step and can suggest a change in dt.
Solving the System • We are going to use an iterative method to solve the following sequence of symmetric matrix problems:
Summary of Temporal Implicit Schemes • Backwards Euler is unconditionally stable for non-negative diffusion parameter D (i.e. any dt>=0) and first order in dt. • Crank-Nicholson is unconditionally stable for non-negative diffusion parameter D (i.e. any dt>=0) and second order in dt. • ESDIRK4 – generalizes to fourth order in dt.
Homework 10 Q1) Modify/use the umDIFFUSION.zip files to solve the diffusion equation with the following time integrators: a) Backwards Euler b) Crank-Nicholsonc) ESDIRK4 Q2) Create a domain and determine convergence rates in h for p=1,3,5,7 for a small dt (i.e. make sure the convergence does not bottom out above 1e-10). Repeat for two of the above time integrators. Q3) Choose p=7 and h small and verify rates of convergence in dt for integration to some fixed time T. Use the one time integrator you did not use in Q2.
Homework 10 • This homework can be completed in pairs or individually. • This homework is due Monday 04/21/03 • Remember – no more than 5 email questions will be replied to for this homework.
Driver for Backwards Euler Diffusion Scheme(umDIFFUSIONdemo.m)
Time Loop (umDIFFUSIONrun.m) 1-12) Set parameters 14-16) Set initial conditions 19) Set solver parameters 22-40) Time stepping loop 28) Solves linear system using conjugate gradient 32-40) Plot solution and compute integral of C
umDIFFUSIONop.m • Set up storage for computing: (M*C +dt*gamma*D* L*C)
umDIFFUSIONop.m 15-25) compute qx,qy using Neumann bcs 27-40) compute div(qx,qy) using Dirichlet bcs 42-48) compute (M*C+gamma*dt*D*L*C)
Discussion of Final Project • Now is the time to ready your final projects. • Submit a one page description of the PDE you intend to discretize with DG by 04/21/03 • This is a proposal and may be rejected – requiring a resubmission or assignment of a set of PDEs. • Include: • List of PDEs to be discretized • List boundary conditions and initial conditions to be discretized • Relevance to your own research (if any) • List of group members (max 3, with one proposal per group) • Interesting issues (application of PML, creation of specific PML, specific physical application/model…)