Flux form semi-Lagrangian transport in ICON: construction and results of idealised test cases

1. Flux form semi-Lagrangian transport in ICON:construction and results of idealised test cases Daniel Reinert Deutscher Wetterdienst Workshop on the Solution of Partial Differential Equations on the Sphere Potsdam, 26.08.2010

2. Outline Flux form semi-Lagrangian (FFSL) scheme for horizontal transport A short recap of the basic ideas FFSL-Implementation approximations according to Miura (2007) Aim: higher order extension Results: Linear vs. quadratic reconstruction 2D solid body rotation 2D deformational flow Summary and outlook

3. Flux form semi-Lagrangian scheme (FFSL) for horizontal transport

4. FFSL: A short recap of the basic ideas • Scheme is based on Finite-Volume (cell integrated) version of the 2D continuity equation. • Assumption for derivation: 2D cartesian coordinate system • Starting point: 2D continuity equation in flux form • Control volume (CV): triangular cells • Discrete value at mass point is defined to be the average over the control volume • Problem: Given at time t0 we seek for a new set of at time t1= t0+Δt as an approximate solution after a short time of transport.

5. FFSL: A short recap of the basic ideas • In general, the solution can be derived by integrating the continuity equation over the Eulerian control volume Ai and the time interval [t0,t1]. • applying Gauss-theorem on the rhs, assuming a triangular CV, … • … we can derive the following FV version of the continuity equation. • No approximation as long as we know the subgrid distribution ρq and the velocity field (or aie) analytically.

6. Physical/graphical interpretation ? Let this be our Eulerian control volume (area Ai), with area averages stored at the mass point control volume

7. Physical/graphicalinterpretation ? Assume that we know all the trajectories terminating at the CV edges at n+1 trajectories control volume

8. Physical/graphical interpretation ? Now we can construct the Lagrangian CV (known as „departure cell“) In a ‚real‘ semi-Lagrangian scheme we would integrate over the departure cell control volume

9. Physical/graphicalinterpretation ? • We apply the Eulerian viewpoint and do the integration just the other way around: • Compute tracer mass that crosses each CV edge during Δt. • Material present in the area aie which is swept across corresponding CV edge

10. Basic algorithm The numerical algorithm to solve for consists of three major steps: example for edge 1 Determinethedepartureregionaiefortheethedge DetermineapproximationtounknowntracersubgriddistributionforeachEuleriancontrolvolume Integratethesubgriddistributionoverthe (yellow) areaaie. Note: For tracer-mass consistency reasons we do not integrate/reconstruct ρ(x,y)q(x,y) but only q(x,y). Mass flux is provided by the dynamical core.

11. FFSL-implementation

12. ApproximationsaccordingtoMiura (2007) • Departure region aie: • aproximated by rhomboidally shaped area. • Assumption: v=const on a given edge • Reconstruction of qn(x,y): • SGS tracer distribution approximated by 2D first order (linear) polynomial. • conservative weighted least squares reconstruction • Integration: • Gauss-Legendre quadrature • No additional splitting of the departure region. Polynomial of upwind cell is applied for the entire departure region. unknown mass point of control volume i

13. Possible improvement - Higher order reconstruction • improvements to the departure regions appear to be too costly for operational NWP • we investigated possible advantages of a higher order reconstruction. Default: first order (linear) polynomial (3 unknowns) Test: second order (quadratic) polynomial (6 unknowns) Stencil for least squares reconstruction CV CV 3-point stencil 9-point stencil

14. How to deal with spherical geometry ? • All computations are performed in local 2D cartesian coordinate systems • We define tangent planes at each edge midpoint and cell center • Neighboring points are projected onto these planes using a gnomonic projection. • Great-circle arcs project as straight lines Gnomonic projection 3-point stencil Plane of projection Equator

15. Results: Linear vs. quadratic reconstruction

16. Solid body rotation test case Example of flow in northeastern direction (=45°) Error norms (l1, l2, l∞) are calculated after one complete revolution for different resolutions • Uniform, non-deformational and constant in time flow on the sphere. • Initial scalar field is a cosine bell centered at the equator • After 12 days of model integration, cosine bell reaches its initial position • Analytic solution at every time step = initial condition

17. Shape preservation • Setup • R2B4 (≈ 140km) • CFL≈0.5 • =45° • flux limiter • conservative reconstruction quadratic L1 = 0.392E-01 L2 = 0.329E-01 linear Errors are more symmetrically distributed for the quadratic reconstruction. L1 = 0.887E-01 L2 = 0.715E-01

18. Non-conservative reconstruction quadratic Setup • R2B4 (≈ 140km) • CFL≈0.5 • =45° • flux limiter • non-conservative reconstruction L1 = 1.293E-01 L2 = 1.133E-01 linear Conservative reconstruction: important when using a quadratic polynomial L1 = 0.854E-01 L2 = 0.699E-01

19. Convergence rates (solid body) Setup: • CFL≈0.25 • =45° (i.e. advection in northeasterndirection) • fluxlimiter quadratic, conservative linear, non-conservative • Quadratic reconstruction shows improved convergence rates and reduced absolute errors.

20. Deformational flow test case • based on Nair, D. and P. H. Lauritzen (2010):A classofDeformational Flow Test-CasesfortheAdvection Problems on theSphere, JCP • Time-varying, analyticalflowfield • Tracer undergoesseveredeformationduringthesimulation • Flow reversesitscourseat half time andthetracerfieldreturnstotheinitialpositionandshape • Test suiteconsistsof 4 casesofinitialconditions, threefor non-divergent andonefor divergent flows. Example: Tracer field for case 1 t=0.5 T t=0 T t=T

21. Convergence rates (deformational flow) • C≈0.50 • flux limiter quadratic, conservative linear, non-conservative • Superiority of quadratic reconstruction (absolute error, convergence rates) less pronounced as compared to solid body advection. But still apparent for l∞. • Possible reason: departure region approximation (rhomboidal) does not account for flow deformation.

22. Summary and outlook • Implemented a 2D FFSL transportscheme in ICON (on triangulargrid) • Based on approximationsoriginallyproposedbyMiura (2007) for hexagonal grids • Pursuedhigher order extensionofthe 2nd order ‚Miura‘ schemebyusing a higher order (i.e. quadratic) polynomialreconstruction. • Quadraticreconstructionledto • improvedshapepreservationandreducedmaximumerrors • improvedconvergencerates (in particular L∞) • reduceddependencyoftheerror on CFL-number • Formorechallengingdeformationalflows, thesuperiority was lessmarked • conservativereconstructionis essential whenusinghigher order polynomials • Outlook: • Whichreconstruction (polynomial order) isthemostefficientone? • Possiblegainsfromcheapimprovementsofthedepartureregions? • Implementation on hexagonal gridcomparison

23. Thank you for your attention !!

24. Model error as a function of CFL (solid body) • Fixed horizontal resolution: R2B5 (≈ 70 km) • variable timestep  variable CFL number • flux limiter Flow orientation angle: α=0° Flow orientation angle: α=45° • Linear rec.: increasing error with increasing timestep and strong dependency on flow orientation angle. • Quadratic rec.: errors less dependent on timestep and flow angle.