1 / 30

MA557/MA578/CS557 Lecture 32

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:.

quin-serrano
Download Presentation

MA557/MA578/CS557 Lecture 32

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. MA557/MA578/CS557Lecture 32 Spring 2003 Prof. Tim Warburton timwar@math.unm.edu

  2. 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:

  3. 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):

  4. Heat Equation • Our first approach will be an explicit time integration of the equation (using zero Neumann bc):

  5. 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:

  6. We also know that all the eigenvalues are non-positive real numbers and dt*rho(A) must fit in the stability region..

  7. 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.

  8. Backwards Euler Time Integration • Suppose we consider the backwards Euler scheme: • Iterator:

  9. 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:

  10. Crank-Nicholson • We replace the right hand side of the discrete scheme with:

  11. 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:

  12. Quick Result • We can bound the iterator 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.

  13. Crank-Nicholson Accuracy • We can briefly compute the temporal accuracy:

  14. 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:

  15. 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)

  16. 4th Order ESDIRK4 Implicit RK Scheme • Iterator: (Ctilde are intermediate values of C, Chat5 is the embedded 3rd order approximation of Cn+1)

  17. 4th Order ESDIRK4 Implicit RK Scheme • If L is a linear operator then:

  18. 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

  19. 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.

  20. Solving the System • We are going to use an iterative method to solve the following sequence of symmetric matrix problems:

  21. 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.

  22. 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.

  23. 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.

  24. Driver for Backwards Euler Diffusion Scheme(umDIFFUSIONdemo.m)

  25. 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

  26. umDIFFUSIONop.m • Set up storage for computing: (M*C +dt*gamma*D* L*C)

  27. 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)

  28. Demo Run

  29. 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…)

More Related