1 / 63

Numerical Methods for Partial Differential Equations

Numerical Methods for Partial Differential Equations . CAAM 452 Spring 2005 Lecture 15 Instructor: Tim Warburton. Two-Dimensional Advection. Recall in one-dimensions we chose an arbitrary section of a pipe. We monitored the flux of fluid into and out of the ends of the pipe.

swain
Download Presentation

Numerical Methods for Partial Differential Equations

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. Numerical Methods for Partial Differential Equations CAAM 452 Spring 2005 Lecture 15 Instructor: Tim Warburton

  2. Two-Dimensional Advection • Recall in one-dimensions we chose an arbitrary section of a pipe. • We monitored the flux of fluid into and out of the ends of the pipe. • The conservation equation in 1D we derived was: • Now we start with a 2D domain and consider any simply connected sub-region:

  3. Our 2D Domain

  4. Our Sub-region

  5. n Outward Normal n

  6. n=outward pointing normal Conservation • Suppose the fluid has a mean velocity everywhere. • This time is a two-vector. • The flux of fluid through the boundary is the integral of the concentration C being advected normal to the surface of w:

  7. Conservation Law • Now we know how much of the concentrate is being advected through the boundary of our sub-region. • The conservation law is now: • We apply the divergence theorem (beware this relies on smoothness arguments):

  8. In Component Form • Assuming ubar is constant this becomes:

  9. Finalize • For all sub-regions: • Assuming enough continuity we obtain:

  10. Basic Technology For Meshes in 2D • Our target is to be able to use triangulations of a domain in 2D. • For example here we have a 5 triangle mesh of an L-shape domain:

  11. Global Vertex Numbering • In order to represent the topology of the triangles we create a global numbering of the unique vertices in the mesh: • Notice – the allocation of number to vertex is arbitrary. 3 5 7 4 6 1 2

  12. Triangle as Ordered Triplet • We next label each triangle in the mesh (numbers in triangles): • Again the labeling of the triangles is in an arbitrary order. 3 5 1 7 2 4 4 5 3 6 1 2

  13. 3 5 1 7 2 4 4 5 3 6 1 2 Triangle as Ordered Triplet cont We can represent each triangle as atriplet of integers which are the globalvertex numbers of the triangle in acounter-clockwise ordering: Triangle 1: 7,5,3 Triangle 2: 2,5,7 Triangle 3: 4,6,1 Triangle 4: 1,7,4 Triangle 5: 2,7,1

  14. Class Exercise Part 1 Given this element to vertex representation of the mesh,can you determine determine the topology of the triangle (i.e. what might the mesh look like)? Triangle 1: 2,3,6 Triangle 2: 5,3,2 Triangle 3: 6,3,1 Triangle 4: 2,6,4 What is the cost of your algorithm?

  15. Class Exercise Part 2 Given this element to vertex representation of the mesh,can you determine which triangles share an edge: Triangle 1: 1,2,9 Triangle 2: 2,3,9 Triangle 3: 3,4,5 Triangle 4: 3,5,6 Triangle 5: 3,6,9 Triangle 6: 9,6,7 Triangle 7: 8,9,7 What is the cost of your algorithm?

  16. Naïve Algorithm • We could naively try a brute force search. • Two triangles will share an edge if they both contain two vertices in common. • So the obvious algorithm is: a b

  17. Alogithm 1 Brute force – test to see which out ofall the edges connects to a given edge. Cost approx. 9*K^2 We can trim the number of testsby ½ easily (loop elmt2=elmt1+1:K and store result in (elmt1,edge1) and (elmt2,edge2). Very loopy so Matlab will be very slow.

  18. Alogithm 2 Set up the edge to node matrix using Matlab’s sparse matrix facility. Thinking about the connectivity matrices as sparse matrices allows us to easily find other connectivities: NodeToEdge = transpose(EdgeToNode) …

  19. Topology Determined • So with either approach we now know the topology of the triangles (i.e. which triangles share an edge). • Consider the advection equation: • We can use this as the foundation of an upwind finite volume solver:

  20. With being the k’th triangle we need to be able to compute and the outward facing unit normals on each edge Geometric Factors

  21. Computing the Area of a Triangle • We can use a simple geometric argument to determine the area of a triangle:

  22. cont • We make some approximations: • that q is in fact constant over the element • Euler-forward in time and we obtain:

  23. Notes on the Geometric Factors We use the formula: area of triangle = ½*(base*perpendicular height) The ratio of the edge length to triangle area gives 2 divided by the perpendicular heightfrom the edge to the vertex not on the edge.

  24. cont • Simplifying: • This looks very much like a one-dimensional scheme applied at each edge – which it is. • At each edge we use the component of the advection velocity in the direction of the outwards facing normal.

  25. Project 1- Overview • I grabbed a polygonal description of part of coast line of the UK and Eire: • I then used the Triangle program by Shewchukto generate a trianglemesh. • Your task: solve the linearizedEuler equations with a pressurepulse set off near the coast. http://www.mar.dfo-mpo.gc.ca/science/ocean/coastal_hydrodynamics/resolute/resolute.html

  26. Project 1 - Equations • The linearized Euler equations we will consider are: • where p is the pressure and u,v are the x and y components of the water velocity. • The boundary condition states that at the boundary the velocity field is tangential to the boundary.

  27. cont • We write a first order finite-volume scheme for the linearized Euler equations using upwind fluxes for: • as (after board demo of characteristic treatment):

  28. Project Tasks Q1) Write code to read in the mesh from a .neu mesh file (more about that next time) Q2) Write start up code to compute: • The area of each triangle • The length of each edge of each triangle • The outwards facing unit normal at each edge of each triangle Q3) Implement the first order finite volume scheme for the linearized acoustics. Use a dt which obeys:

  29. Project Tasks cont Q4) By generating a sequence of increasingly refined meshes for a square [-1,1]x[-1,1] find the solution order of accuracy of the scheme. A suitable initial condition is: The solution for all time is:

  30. Project Tasks cont Q5) Use the ocean mesh and model a pressure pulse (Gaussian profile) off the coast of England • Repeat with 3 increasingly refined meshes to visually display and compare results. • Plot snap shots at 4 time intervals showing the pulse spreading and interacting with the coast line. • I will now demo how to build meshes with some in house software based on Shewchuk’s triangle mesh generator.

  31. Project Details • You may work in groups of 1 or 2. • No groups of 3 or more will be permitted. • The project is due: 03/29/05 • Be prepared to present your results (i.e. make a 5 minute PowerPoint presentation). • Prepare a report – and I strongly suggest using the Latex stylesheet available from the website.

  32. More Advanced Than FV • As you will find out the finite volume method is robust, but not very accurate. • To increase the solution order of accuracy we are going to use the discontinuous Galerkin method (DGM or DG) • The idea: on each triangle we create a p’th order polynomial expansion local to the triangle. • We then use a flux type formulation to exchange information between triangles.

  33. Basics • First we need to discuss how to create local polynomial approximations on each triangle. • i.e. we need a robustly computable basis for • In finite volume and finite element methods it is common to perform calculus type operations (integration, differentiation and interpolation) on a standard element and map the results to a physical element.

  34. Reference Triangle • The following will be our reference triangle: • The k’th triangle is the image of this triangle under the map: s (-1,1) r (-1,-1) (1,-1)

  35. Reference Triangle Mapped To Physical Triangle • A typical map will have the following action: s (-1,1) r (-1,-1) (1,-1)

  36. Orthonormal Basis for the Triangle • Fortunately a well behaved, orthonormal basis for the triangle has been discovered (and rediscovered multiple times). • First we need to know some details about the Jacobi polynomials. These polynomials are parameterized by two reals: and their integer orders n,m such that they satisfy the orthogonality relationship (for integer alpha,beta):

  37. Jacobi Polynomial Recurrence Relation We can generate the n’th order Jacobi polynomials at a given x in [-1,1] through the following recurrence relationship: See: http://www.caam.rice.edu/~timwar/MA578S03/MatlabScripts/JACOBI1D.m

  38. Orthonormal Basis for the Triangle • The following basis is due to Koornwinder (later revived by Proriol, Dubiner, Owens,….) • It satisfies the following orthogonality condition:

  39. cont • Notice that the (n,m) basis member is an n+m’th order polynomial. • Moreover, the orthogonality property allows us to determine that the set of basis members with: • n+m<=p supports all bivariate polynomials with maximum total degree p.

  40. Basis Ordering In order to simplify the notation we can order the polynomials with a single index (this is not a unique ordering): where M=(p+1)(p+2)/2 and we have normalized each of the

  41. cont • Next class we will discuss how to use this basis to interpolate, differentiate and integrate data at points in the triangle. • We continue here to discuss the question of which points to use on the triangle for polynomial interpolation.

  42. Given a set of nodes lying in the triangle we use V to construct an interpolating polynomial for a function who’s values we know at the nodes: The interpolation condition yields: Interpolation Using Generalized Vandermonde Matrix

  43. Differentiation • Suppose we wish to find the derivative of a p’th order polynomial: • First we note that the approximation becomes equality: • And interpolation allows us to find the polynomial coefficients: • So differentiation requires us to compute:

  44. Differentiation cont • So we need to be able to compute: • Recall the definition of the basis functions: • R-derivative:

  45. Quick Jacobi Polynomial Identity • We will make extensive use of the following:

  46. r-Derivative • Ok we need to calculate: • We can compute these using the definition of the Jacobi polynomials. • Watch out for s=1 (top vertex) – the r-derivative of all the basis is functions is zero at r=1,s=-1

  47. s-Derivative We use the chain and product rule to obtain:

  48. s-Derivative From which:

  49. Special Cases Don’t worry about all those denominators having (1-s)since the functions are just polynomials and not singular functions…

More Related