Documenting advection schemes

[glwagner]

We should document the advection schemes that are available.

[navidcy]

I can help with that if you point me to a source/wiki/somewhere explaining exactly what each scheme does.

[glwagner]

I don't know a good reference for everything...

To start with, our object is to calculate advective fluxes across cell interfaces. For this purpose we need to know the velocity field and the advected quantity (usually a tracer or velocity component) at the cell interfaces.

The tracer schemes are simple, because by construction the velocity field is located at cell interface and no interpolation is required for the velocity field. Only the tracer field is interpolated. In this case the advection scheme corresponds to the interpolation method. The interpolation method "reconstructs" the value of the tracer field at cell interfaces, given knowledge of the cell-averaged tracer values on a regularly spaced grid. It's important that the stencils are specific to reconstruction using cell-averaged values (ie, a finite volume distribution) rather than tracer values at nodes. CenteredSecondOrder just averages the cell-averaged tracer values on either side of the interface. UpwindBiasedThirdOrder incorporates an addition cell in the upwind direction. CenterdFourthOrder uses a fourth order stencil for a cell-averaged finite volume representation. UpwindBiasedFifthOrder uses 5 cells around the cell interface, biased in the upwind direction. WENO5 combines three third-order stencils using a smoothness indicator that biases the stencil away from discontinuities and rapid changes in the underlying tracer distribution.

The advection scheme for momentum is more complicated because both the advected quantity and advecting velocity field have to be interpolated. For this we use symmetric interpolation (even-ordered interpolation) of the advecting velocity, and the same interpolation scheme used for tracers for the advected quantity.

For a centered advection scheme (CenteredSecondOrder, CenteredFourthOrder) the symmetric interpolation for advecting velocities is identical to interpolation used for advected quantities.

For an upwind-biased advection scheme (UpwindBiasedThirdOrder, UpwindBiasedFifthOrder, WENO5), the symmetric interpolation for advecting velocities is a symmetric scheme of order n-1. So that's CenteredSecondOrder, CenteredFourthOrder, and CenteredFourthOrder for the three schemes, respectively.

The symmetric interpolation of advecting velocities for upwind biased schemes is described in Ghosh and Baeder 2010:

[navidcy]

This is helpful. I didn't really mean one reference for all... What you've put here is a good start. I'll start a PR and we can continue chatting there.

[francispoulin]

I am glad that we are talking about this and I'm also happy to help if I can.

One thing that I noticed on the docs is that we say we are using the finite volume method but then we never actually integrate the PDEs to obtain the equations in terms of the cell averaged quantities. I think this would be helpful to the user as it would point out the differences between the value at a point and the cell average of that quantity. I don't think there needs to be a lot but integrating the tracer equation first, since that's easier, and the integrating the momentum equation would add someting which at that moment I don't think is present.

What do you think @navidcy ?

[glwagner]

The classic WENO reference is https://link.springer.com/chapter/10.1007/BFb0096355

I think there are some references in the code.

A derivation of discrete equations for cell-averaged quantities is a good idea too.

[francispoulin]

This is a rough first attempt to show how we derive the finite volume formulation for the tracer equation.

Two questions:

• Is this at all what we might want?
• It seems to me that the velocities are evaluated at the edges but in the case we use the volume averaged velocities. How do we get this in terms of volume averages of velocity?

[navidcy]

This is a rough first attempt to show how we derive the finite volume formulation for the tracer equation.

Two questions:

• Is this at all what we might want?
• It seems to me that the velocities are evaluated at the edges but in the case we use the volume averaged velocities. How do we get this in terms of volume averages of velocity?

Finite_Volume_Tracer_Equation.pdf

Why do you say that you "assuming regular spacings"? I see symbols like Δx_{i,j,k} which does not imply that the spacing is the same for every cell...

[francispoulin]

In the last equation I am basically cancelling out \Delta A/ V to get a length scale but if the area is different at either ends then this isn't quite true. For a rectilinear case, which is the one to explain hear, there is no problem. I was maybe worried about the curvlinear case but that is something to think about for the future.

I will remove that bullet point.

[navidcy]

oh I see...

This is a good, but it should come with a picture defining what each index means. :) Don't worry though for now. Or perhaps a slightly different notation instead of i, j, k?

And then we can explain the various ways advection schemes work just by focusing on one of the terms, e.g., ∂ₓ(uc).

[navidcy]

I'm tempted to have the derivation in the Docs only for 1D case...

3D/ND generalizations follow trivially from that... They are just laborious. Am I missing something?

[francispoulin]

You are right that the notation needs to be cleaned up and it would be nice of it's is consistent with the code.

I also agree that 3D is more complicated then we need. How about a compromise and we do 2D? It is true that there is nothing really new going from 1D to higher D, but its easy enough to do a picture in 2D and I think that it's important to show the different area elements. But I'm happy for us to do anything so if people prefer 1D, then lets do 1D.

[navidcy]

Perhaps 2D might be better indeed...

[glwagner]

I think even 1D would be ok.

On the question of interpolation of velocities: to compute fluxes we need to reconstruct the fields at the cell interfaces. The velocity field is staggered with respect to the tracer field, such that we can trivially reconstruct the velocity field by using the cell-averaged velocity component that's centered on the location it's needed.

[glwagner]

I think your derivation is valid for locally orthogonal grids, not just rectilinear grids (for example, in a curvilinear orthogonal domains you can also use the three velocity components to compute the advective tracer flux).

[francispoulin]

Having to reconstruct the fields from the cell averages makes sense. But when I look in tracer_advection_operators.jl, it seems like we are using the velocitys that are stored, which I thought were the cell-averaged quantities. I believe I'm missing something?

    1/Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, advective_tracer_flux_x, advection, U.u, c) +
                                   δyᵃᶜᵃ(i, j, k, grid, advective_tracer_flux_y, advection, U.v, c) +
                                   δzᵃᵃᶜ(i, j, k, grid, advective_tracer_flux_z, advection, U.w, c))

[navidcy]

You are saying that you’d expect to have estimates of the velocity values at the interfaces rather than their cell-averaged values (centered around the interfaces)?

[francispoulin]

Yes. Given the derivation I posted, and what @glwagner said, it seems like we want to be using the velocity at an edge not the cell-averaged values, which is what we are storing as fields. But I could be missing something here.

[francispoulin]

I am including an update of my notes. Differences of note:

• Doing everything in 2D, but can downgrade to 1D if people feel this is confusing. I thought this would work well with the 2D picture we have already
• Using overline to denote a cell-averaged quantity. I know the code doesn't do this but this is the convention in the literature and seems good to have something to distintuish between cell-average and simply value
• The final answer is indeed in terms of quantities are the edges but the final equation is exact. The approximations come in when we approximate those fluxes using the cell-averaged quantities. We have a bunch of ways of doing this and I for one would like to discuss this more, but thought we should discuss what I have so far before going any further.

[glwagner]

Might make sense to include a diffusive term too? The second-order fluxes are crucial and usually mediate how boundary conditions are prescribed.

[francispoulin]

Good idea @glwagner. Below is a revised copy that includes simple diffusion.

If we think this is a good start I can also create a PR with this in the docs.

[glwagner]

Nice! We usually use \nu for viscosity and \kappa for diffusivity --- so should we use \kappa here? You can allow it to depend on x, y, z by writing \nabla \cdot (\kappa \nabla c) without loss of generality (it's also nice to see the divergence operator show up twice I think). Thanks for putting this together...

[glwagner]

@francispoulin if this doc is on overleaf I'm happy to help too.

[francispoulin]

My bad, yes, \kappa is much better. Fixed.

Also, I rewrote the Laplacian as the divergence of the product of \kappa and the gradient of c. That is also better I agree.

Should I put together a PR for this part and then work on the momentum? I'm also happy to wait for more feedback. @navidcy , what do you think?

@glwagner : I will send you a link to the overleaf document.

[glwagner]

Yes. Given the derivation I posted, and what @glwagner said, it seems like we want to be using the velocity at an edge not the cell-averaged values, which is what we are storing as fields. But I could be missing something here.

I think we are making a particular choice: because we are using a staggered grid, we can reconstruct the velocity field at tracer cell interfaces simply be evaluating the cell-averaged velocity field there. Perhaps this reconstruction has a particular order of accuracy (eg second-order?) I'm not sure.

[glwagner]

Should I put together a PR for this part and then work on the momentum? I'm also happy to wait for more feedback. @navidcy , what do you think?

My two cents is that it might be better to add as much content as possible in one PR. I thinks docs PRs can be a bit painful. It's not as crucial either for docs updates that the PRs are tight and focused.

I think it's simpler to keep things general for curvilinear grids than to assume constant cell spacings by the way.

We should discuss 2D versus 3D a bit more. While 2D is simpler, we will almost certainly need 3D as well at some point (when we resolve #1679). So we can have a "tutorial" 2D section, and then a "fully descriptive" 3D development --- or we can just do 3D from the start.

[francispoulin]

I'm still trying to figure out exactly how to go from cell value to cell-averaged value and back again. The forward approach is the quadrature scheme that we use to integrate on our cells. The backwards approach is the reconstructor.

The standard choices that I've seen are to assume that the function is constant, linear and parabolic on each cell. If it's constant then the two coincide and there is nothing to worry about, but I presume this yields a lower order of accuracy.
On top of what , we then have another choice on how to choose the flux based on the advection scheme.

Can you help me figure out what we are using as a reconstruction?

Everything else you suggest sounds good to me.

[glwagner]

Can you help me figure out what we are using as a reconstruction?

For the velocity field or tracer?

For the tracer it depends on the advection scheme.

For the velocity field, I think our method is consistent with either a constant or linear reconstruction. I'm not sure about higher.

Not that if we use a different reconstruction for computing tracer fluxes, we'd have to use the same reconstruction for calculating the divergence of the velocity field so that they are consistent.

[glwagner]

I think section 2 of this note has a useful description of the reconstruction process:

http://ftp.cc.ac.cn/lcfd/DEWENO/paper/icase-1997-65.pdf

I haven't found any references regarding velocity field reconstruction for tracer advection / mass conservation in the context of a staggered grid scheme.

[navidcy]

I missed all this in my sleep :( Sleeping makes you miss the fun..

[navidcy]

My bad, yes, \kappa is much better. Fixed.

omg, definitely κ :)

Should I put together a PR for this part and then work on the momentum? I'm also happy to wait for more feedback. @navidcy , what do you think?

Yeap, why not. Start a PR.

@glwagner : I will send you a link to the overleaf document.

Add me to that. Or better just start the PR and we can all edit there.

[navidcy]

~I think this is a bit better reference.~

~Shu, Chi-Wang. "High order weighted essentially nonoscillatory schemes for convection dominated problems." SIAM review 51 (1) (2009): 82-126.~

On second thought, this ref discusses also finite differences schemes as well as finite volume schemes. So let's stick with FVM-related refs.

[navidcy]

[cc @dhruvbhagtani2105 in case he's missing the fun also.. :)]

[navidcy]

OK. Sorry for not chiming in earlier @francispoulin, but indeed you are correct. The derivation you have for the advective terms of the tracer equation requires, e.g., u at the cell's interfaces. But we don't have those values, as you point out. In a finite volume formulation we only have the cell-averages of u around cells with centers the cell interfaces. But there is a way to reconstruct the value of u @ the interface from knowledge of the cell-averages of ū. This is usually referred as "Reconstruction" (section 2 of Shu 2009). I guess the best way to reconstruct u @ interface from all ū's is to take u as the cell average of ū at the particular cell? Have to think a bit on that (and wait for the coffee to kick in), but I guess this holds.

[navidcy]

Let's do this exercise: if we have the cell averages of a field ūⱼat points xⱼ then what is the n-th order estimation of u(xⱼ)? I guess this will answer the question?

[francispoulin]

@navidcy

I just sent you an invite to the overleaf document. I thought maybe we want to have both the tracer and momentum equations before we create a PR, but I'm flexible.

Thanks also for the reference. I remember seeing this paper before and liking it so I will look at it again.

It seems that our advection schemes use the cell-averaged velocity to approximate the velocity at the edge. This is a perfectly valid choice and corresponds to a certain order of accuracy. I don't know what that is but I imagine using more cell averaged velocities would give a higher order accuracy. I'm not suggesting we do this, but I'm just trying to get a better understanding of the order of accuracy.

[francispoulin]

Both the references that @glwagner and @navidcy shared are the classical ones so you can't go wrong. However, their formulas don't apply directly as they do not use a staggered grid.

An approach that uses WENO on a staggered C-grid can be found here. In this work they also take the cell-averaged velocity to compute the flux and only use multiple values of the tracer to get high order accuracy in WENO. They cite the same papers so maybe I should read them more carefully to figure out how to deal with a staggered grid.

[glwagner]

Mishra et al 2020 (referenced in #1705) discusses how the velocity reconstruction scheme we use for tracers is second-order accurate (in Appendix A they give additional references for such a discussion). In section 5 they define a higher-order scheme for this purpose.

As explained there, almost all implementations with WENO are effectively 2nd order (for both momentum and tracers). It's unclear how much is gained from higher-order accuracy, either for momentum or for tracers. The asymptotic convergence rate is just one factor in determining the accuracy of a numerical solution, I guess. Resolving #1705 would be one step toward figuring that out.

[glwagner]

Skamarock and Klemp (2008) state:

The order of accuracy is, however, only formally third order for constant U; for non-constant velocities the scheme is only second-order accurate.

Text here:

[navidcy]

Even just something within the docstrings of the advection scheme composite type definitions/constructors, like we did for WENO would be very beneficial.

cc @simone-silvestri