Numerical Methods for Partial Differential Equations

1. Numerical Methods for Partial Differential Equations CAAM 452 Spring 2005 Lecture 15 Instructor: Tim Warburton

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

4. Our 2D Domain

5. Our Sub-region

6. n Outward Normal n

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

49. s-Derivative From which:

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