A novel VIP (Vector Image Polygon) slope limiter for scalar variables in multidimensional grids

1. “Numerical Methods for multi-Material Fluid Flows”, Arcachon, France,5-9 Sept 2011 A novel VIP (Vector Image Polygon) slope limiter for scalar variables in multidimensional grids Gabi Luttwak1 and Joseph Falcovitz2 1Rafael, P.O. Box 2250, Haifa 31021, Israel 2Institute of Mathematics, The Hebrew University of Jerusalem, Israel

2. VIP-Convex Hull based limiters • What is VIP? • the natural extension of the monotonicity from scalar to vector variables • VIP => for the Lagrangian phase of the SMG/Q scheme (Pavia 2009) • VIP => for the advection phase of ALE/Euler SMG calculations (Paris 2010) • VIP => for advection of scalar variables ? • In the advection phase of ALE/Euler SMG/Q calculations • Test cases: • Noh 2D, Sedov 2D/3D, Standing wave in Saltzmann mesh • New “Noh Corner” Test

3. In Finite Volume Schemes • The conservation laws solved over a control volume : • (e.g. the computational zone - for zone centered variables) • Fluxes evaluated on the control volume faces • Face-centered values of variables are required

4. Slope Limiters • Second order Godunov and some other high resolution schemes use the gradients of the variables to get face-value: • In hyperbolic PDE discontinuities (shock or slip lines) are present, or can be generated, during the flow. • At discontinuities gradients are unbounded. • Using them produces unphysical fluctuations in the solution • Thus either gradients or the related fluxes must be limitedto preserve a monotonic flow field

5. Slope Limiters in 1D • Limiters were first formulated in 1D: • Either the slope-extrapolated values of a variable at the zone faces must lie in the range spanned by the values of that variable in neighboring cells: • Or the limited gradient defined as some kind of monotonic averaging of the gradients in the zone c and its neighbors . e.g. van Leer’s double minmod:

6. Limiters for scalars in 2D/3D • In 1D the two approaches are equivalent • Equation (1) is the usual way to extend the limiters to 2D/3D • See Dukowicz, Barth&Jesperson • Equation (2) is van Leer’s double minmod • The gradient of scalars are vectors • What would be monotonic averaging of vectors in 2D/3D ?

7. How is (1) extended to 2D/3D • Look at zone centered scalar variables • Zone face in 1D=>either face/edgecenter or zone vertex in 2D/3D • Extrapolate using (1) from the zone center to each of its vertices • Find a limited gradient • Such that all vertex extrapolated values of the density remain within the density range spanned by values in the zones around that vertex

8. Extending Equation (1) to 2D/3D(a) Scalar Limiting • For each vertex v find the largest such that: • And take: • Straightforward but dissipative !! • This was our default scheme for scalar variables

9. Extending Equation (1) to 2D/3D(a) Non scalar α • Find obeying Eq. (1) at each vertex ʋ “closest” in some sense to the unlimited . • Limiting should be preferentially normal to shock direction, so that the gradient direction may change: • Berger et.al. - Linear programming • Ridder&Kothe – Least squarereconstruction

10. Let us extend Eq.(2) to 2D/3D ! • We apply the VIP convex hull based monotonicity criterion to the slopes: • In 1D: • In 2D/3D: • If a slope is outside the VIP generated by twice the slopes of its neighbors=>move into the hull • Factor 2 like in the double minmod • e.g. we move it to the nearest point on the VIP

11. What is the analog of ? • In 1D corresponds to an extremum • (a) In 2D/3D extremum is when . • In the vector space of the slope-vectors, extremum exists if the origin (0,0,0) is inside the VIP => set slope to zero • Difficulty: “Ridge configuration” a small slope component in some direction may prevent detecting a discontinuity normal to a major gradient

12. VIP slope limiter for scalars (b) At a surface discontinuity is the normal to the surface • changes sign along • Take in the direction of the unlimited gradient

13. VIP slope limiter for scalars • Let the unlimited slopes in cell c and its k neighbors. • Let and • with • If => zero the slope: • Else if lies outside the VIP (convex hull ) move it to the closest point on hull

14. 3 1 4 5 2 6 7 9 12 11 8 A Convex Hull (CH) 10

15. The VIP convex hull limiter V0 “Move to hull” V1 V5 Vc V2 V0 V6 V4 V3 • if V0 inside VIP => no limiting: • if V0 outside, VIPlimiter is applied and V0 => Vc

16. The Staggered Mesh Godunov (SMG/Q) Scheme SMG schemes for 3D Lagrangian and ALE hydrodynamics were presented at Oxford 2005, Prague 2007, Pavia 2009,APS 2003,2005. Use Riemann problem (RP) solutions to capture shocks. Vertex-centered velocities have jumps on the in-cell corner zone faces. These are simplified, “impact” RP or (IRP) with continuous . The IRP are solved in the normal to shock direction, i.e. along the velocity difference . This defines a uniaxial tensor pseudo-viscosity.

17. The Staggered Mesh Godunov (SMG/Q) Scheme The velocities on either side of the face, serving as data for the RP, are evaluated from the vertex velocities using a cell centered velocity gradient, limited to preserve a monotonic velocity distribution. A limited velocity gradient is also used to compute the momentum fluxes during the advection phase The resulting scheme captures shocks with sharp monotonic profiles. It also has a strong inherent mesh stabilizing effect. The rotation invariant VIP limiter is used to preserve the symmetries present Limited slope Original slope Velocity jump for RP L C R

18. Test Problems • Saltzman, Sedov, Noh problems test the propagation of 1D shocks over non-aligned or distorted meshes • Our corner Noh (a 2D/3D collision Riemann problem) tests 2D/3D sector-wise converging flows • VIP for scalars is used only in the advection phase • ALE could smooth the Saltzman mesh making it regular • Sedov test will lack resolution at the shock in an Eulerian or smoothed out non-adaptive ALE calculation • In the examples VIP for scalars is compared with using the convex hull limiter only for vectors and “scalar limiting” (Eqs.1,3,4) for scalar variables

19. The Standing Shock Saltzman Test • In the original test, a piston moves with a velocity of u=10 into an almost cold ideal gas. This 1D problem is solved on a 2D 100x10 mesh as shown. • Instead, we look at the same mesh fixed in a frame of reference in which the shock is stationary with the gas flowing in from the right and out to the left

20. Saltzman with standing wave T=0.4isodensity plot(a) VIP vectors & scalars (b) VIP vectors

21. Saltzman with standing wave T=0.4Profiles of density and pressure(a) VIP vectors & scalars (b) VIP vectors

22. 2D Noh Test • A disk of cold ideal gas collapses to its axis with • A 100x100x2 disk is meshed with 50x50x1 zones • The analytic solution is a shock expanding with • and density profile: • ALE motion is limited to preserve resolution at and behind the shock

23. Noh 2D T60 iso-density plot(a) VIP vectors & scalars (b) VIP vectors

24. Noh 2D T60 iso-density plot inside(a) VIP vectors & scalars (b) VIP vectors

25. “Noh Corner” Test • A disk of cold ideal gas moves toward the origin with the velocity directed in each quadrant along the diagonal as shown • Due to the symmetry only a quadrant has to be considered. We carry out this simulation in a frame in which the gas is initially stationary and the walls are moving • In a “hot gas” corner flow and

26. “Cold” Noh Corner Test Iso-densities at T=75.0100x100x1 zones(a) VIP vectors and scalars (b) VIP vectors only a b

27. “Cold” Noh Corner Test Iso-densities at T=75.0VIP vectors and scalars (a) 200x200x1 zones (b) 100x100x1 zones

28. “Hot” Noh Corner Test Iso-densities at T=70.0100x100x1 zones(a) VIP vectors and scalars (b) VIP vectors only

29. Sedov-Taylor blast wave problem • The initial data is: cold ideal gas with a “point source” in a single zone at the origin • The physical space is a 253 boxdivided into a uniform 903 mesh with 3D computation conducted in the octant with the symmetry planes at • The exact solution has a post-shock density and the front radius is at . • To keep the mesh fine near the shock, the ALE motion was restricted to start at densities

30. Sedov-Taylor 3D Test T=1 VIP vectors & scalars isodensity plot

31. 2D Sedov-Taylor test • The initial data is: • To remove any effect of boundary conditions we ran the simulations in all 4 quadrants: • We set at the 4 central zones • The stronger cylindrical shock reaches at • The ALE motion was again restricted to

32. Sedov-Taylor 2D Test T=0.24isodensity plot(a) VIP vectors and scalars (b) VIP vectors only

33. Sedov-Taylor 2D Test T=0.24isodensity plot inside (a) VIP vectors and scalars (b) VIP vectors only

34. Conclusions • We have devised a new VIP based limiter for the gradient of scalar zone centered variables • This limiter is applied if either the unlimited gradient lies outside the VIP spanned by (twice) the gradients in the neighbor zones or if changes sign • Near a shock the second criterion dominates, so that we do not have to compute the convex hull and we zero the limited gradient. Thus for this case the results are identical to the scalar limiter (Eq. (1),(3),(4) )

35. More Conclusions • In principle the new limiter should be less dissipative than the “scalar limiter” • While we see some differences in the finer details pointing to this, in the test cases shown the results using both limiters are quite similar, virtually identical. • On one hand the new limiter performs well in all the test cases, on the other hand we will consider more test cases where it may have a clear advantage over scalar limiters.

36. Future Perspectives • This limiter is applied once to the zone centered gradient. The same limited gradient is used for fluxes at all zone faces. • Other variants, e.g. a directional limiter, should also be considered

37. References • Luttwak G. , Falcovitz J., ”Slope Limiting for vectors: A Novel Vector Limiting Algorithm”, Int.J.Num.Meth.Fluids,65,p1365-1375,(2011), presented at the Conf. Num. Meth. for Multi-Material Flows , Pavia, Italy, (2009) • Luttwak G. , Falcovitz J., “Vector Image Polygon (VIP) limiters in ALE Hydrodynamics”, EPJ Web of Conferences,10,00020,(2010),presented at the 8th Int. Conf. on New Models and Hydro-codes, Paris, France, (2010) • M.Berger, M.J.Aftosmis, ”Analysis of slope limiters on Irregular Grids”, NAS Technical Report NAS-05-007,[2005] • W. J. Rider, D. B. Kothe, “Constrained Minimization for Monotonic Reconstruction,” AIAA-97–2036, 1997. • Dukowicz J. K. and Kodis J. W. “Accurate conservative remapping (rezoning) for ALE computations, SIAM J. Scient. Stat. Comp., (8),p305-321,1987 • Luttwak G. , Falcovitz J., ”Applying the SMG Scheme to Reactive Flow”, J. Energetic Materials, 28(1), p279 - 302,(2010), presented at the 7th Int. Conf. on New Models and Hydro-codes, Lisbon,Portugal,2008

38. References • Luttwak, G. ”Staggered Mesh Godunov (SMG) Schemes for Lagrangian Hydrodynamics”, p339-342, Shock Compression of Condensed Matter-2005, Furnish M. D. et al Eds, (2006), AIP, CP • Luttwak G., “ Sliding and Multifluid Velocities in Staggered Mesh (MMALE) Codes ", Conference/Workshop On Numerical Methods For Multi-Material Flows, Prague , Sept. 2007, www-troja.fjfi.cvut.cz/~multimat07 • Luttwak G. , Falcovitz J. , "Staggered Mesh Godunov (SMG) Schemes for ALE Hydrodynamics", Workshop On Numerical Methods For Multi-Material Flows , Oxford, UK, Sept. 2005 , www.extra.rdg.ac.uk/ifcd/Multimaterial-Workshop.htm • Luttwak, G., "Comparing Lagrangian Godunov and Pseudo-Viscosity Schemes for Multi-Dimensional Impact Simulations”, p255-258, Shock Compression of Condensed Matter-2001, Furnish M. D. et al Eds, (2002), AIP, CP620.

39. References • Christensen R.B., "Godunov Methods on a Staggered Mesh. An Improved Artificial Viscosity", L.L.N.L report UCRL-JC-105269, (1990). • Caramana E.J., Shaskov M.J., Whalen P.P., J. Comp. Phys. 144, p70, (1998). • D.J.Benson,”Computational Methods in Lagrangian and Eulerian Hydrocodes”, Comp.Meth.Appl.Mech. Engn.,99 (1992) ,p235-391

40. How does a slope limiter work? Familiar algorithms are formulated for scalar variables. The slope extrapolated values of a variable must lie in the range defined by the values of that variable in neighboring cells ν: For vectors, limiting is usually applied separately to each component. Such procedure is frame-dependent: Rotating the coordinate axes produces different results Component limiters do not preserve problem symmetry They do not transform like vector entitiesshould

41. 4 6 1 1 1 1 3 2 2 2 2 8 5 The SMG VIP Vector Limiter • Both for the Lagrangian and Advection phases of the SMG scheme (in 2D/3D). • The extension to a general connectivity (FEM) mesh is straightforward • The limiting is done separately along each edge of a cell (e.g. 1-2 in the figure) : • The outward extrapolated velocities must lie in the VIPs of edge nodes 1,2 • The limited gradient assumed to lie along the velocity difference 7

42. 4 6 1 1 1 1 3 2 2 2 2 8 5 SMG VIP vector limiter • The velocities of 1,2 are extrapolated outward along the edge. • e.g. velocity of 2 extrapolated toward middle of 23: • The cell centered velocity gradient is used, taking its component along the velocity difference : 7

43. 4 6 1 1 1 1 3 2 2 2 2 7 8 5 SMG VIP limiter • The value of chosen to keep inside the VIP of the edge neighbors (1,5,3,4) • The same way chosen to keep inside the VIP of the edge neighbors (2,6,7,8) • The VIP limited slope along 12 generates a velocity jump between the corner zones : where: • The above jump serves as the data for the IRP in the Lagrange phase

44. 2 3 1 O 4 8 5 6 7 Momentum Advection on Staggered mesh • Momentum integrated around corner zones around vertex O (12345678 in the figure)

45. The momentum advection phase • Momentum fluxes at the mid-faces between the corner zones found using the limited velocity profile • The vertex mass flux between these corner zones found from the mass flux through the zone faces (similar to Benson HIS) • The upstream weighted velocity at the middle of the fluxed volume is: here d,a are the donor and acceptor sides (1or 2) • And the momentum flux =>

46. The Staggered Mesh Godunov (SMG/Q) Scheme • Let p* be the RP solution pressure. We take as a uni-axial tensor pseudo-viscosity acting along the shock direction. • Its impulse is imparted to the two neighboring vertices. • Its work, which must be dissipative, is added to the zone internal energy. • Some of these ideas are related to Christensen’s split-Q, and to the edge viscosity and compatible hydro-scheme by Caramana et al.

47. When is a vector “monotonic” relative to its neighbors? • We seek a monotonicity preservation criterion for a vector, analogous to the scalar criterion: • We have defined VIP - Vector Image Polygon in 2D (Vector Image Polyhedron in 3D) as the convex hull of the vector-space points corresponding to the neighbor vectors. • If a slope-extrapolated vector lies inside the VIP, the slope is monotonicity preserving. Otherwise, slope limiting is required. • If the original vector lies entirely outside the VIP, we zero its slopes, in analogy to the scalar extremum case

48. 1 A 2 3 The VIP Monotonicity Criterion • A point inside a convex hull can be expressed as: • In 1D (n=2) this reduces to: which is equivalent to the scalar criterion (1) • A vector point lying on the line connecting two neighbor points should be monotonicity complying. Since VIP is convex any segment connecting two vertices lies inside the VIP (or on boundary). •  The VIP monotonicity criterion is a natural extension of (scalar) monotonicity to vectors.

49. Falcovitz: Cold Noh Corner Test1000x1000 Eulerian GRP calculation