Arbitrary height-dependent blending between advection schemes

[ewquon]

For LES studies, this capability is needed to maximize the effective resolution while preventing numerical instability and spurious numerical artifacts. Spectral resolution for WRF is ~7dx whereas other codes can achieve much better effective resolution of ~4dx. Higher effective resolution can be achieved with higher-order schemes with judicious use of upwinding in selected parts of the flow (e.g., in the free atmosphere above the ABL).

An example numerical setup (from @mchurchf):

The interpolated quantity should be f(z)*adv_scheme1 + (1-f(z))*adv_scheme2. Make f(z) = nearly 1 up to the capping inversion height, then have it ramp to something like 0.75 over the next 100 m and hold constant at 0.75 above that. It would be nice if it can be arbitrary enough that you can blend any 2 of the available advection schemes.

This should be pretty straightforward to implement but I'm not sure what the best way to do it while maximizing computational efficiency would be.

[mchurchf]

The reason this is suggested is because this is what we've done with good success over the past decade or more in SOWFA. We use the local advection scheme blending capability in OpenFOAM. The code that provides this local blending capability is here: https://cpp.openfoam.org/v6/localBlended_8H_source.html. You just provide the localBlended scheme with a "blendingField." Leaving that "blendingField" very arbitrary is useful. At some point in the future, you may want to make the "blendingField" a function of some other arbitrary thing like x,y,z,t, grid cell size, or even some solution variable, so keeping it general is good. Below is a picture from a recent SOWFA/OpenFOAM simulation showing the blended advection scheme. The contours on the left are velocity into the screen. They are nice and crisp in the ABL and oscillation free above the ABL. The contours on the right are the blending factor field which in this case is only z-varying, but we can make it completely arbitrary and user-defined and even something that updates in time.

[AMLattanzi]

Thanks @ewquon and @mchurchf for the suggestion and examples. This sounds like a very promising pathway. From an ERF perspective, it sounds like we want to grab 2 horizontal advection operators and 2 vertical advection operators to do the blending. Since we have modularized the advection routines, select cases are currently used to grab the 1 horizontal advection operator here and 1 vertical advection operator here.

It seems to me that we could add additional functions to the AdvectionSrcForMom_N/T.H files that choose the second operators in the horizontal and vertical direction. This should be mainly copy and paste of code and the addition of template parameters. Furthermore, this should be compatible with all the advection operators ERF currently supports. With regard to the f(z) blending function, do you anticipate this being updated once per whole time step or once per RK stage? Additionally, would a device vector suffice or do you want this to vary in the horizontal direction? The only computational drawback that is jumping out to me is the memory access for adv_scheme2 that could occur in regions where (1-f(z)) ~= 0; for high-order methods you will touch a lot of data but the overall impact to the interpolated quantity might be low. Perhaps we could determine a stability/accuracy threshold for f(z) and choose when to call adv_scheme2 versus using the value already obtained from adv_scheme1.