Vertical Lagrangian-remap design doc

[cbegeman]

This design doc is for the implementation of the vertical Lagrangian-remap method. This development is motivated by the occurrence of instabilities in cases of evolving ice-shelf geometry, but should have broad advantages including higher-order accuracy and lifting the vertical CFL constraint.

[xylar]

Looks good, thanks for the new PR!

[dengwirda]

Thanks for putting this together @cbegeman!

[xylar]

@cbegeman, thanks for getting this going. We will just keep this PR open and the conversation going until the design doc is completely finished, so it won't get merged anytime soon. For now, let's focus on the requirements and maybe start sketching out the Algorithmic Formulation / Design Solution for each of them. Implementation and testing usually get left until later but it doesn't hurt to start playing around with prototype implementation even before we get all the requirements and algorithms nailed down.

[xylar]

Ooh, good ones! Conservation and monotonicity should both be added. I would think we would require conservation in an integral sense (e.g. column volume is the same to machine precision before and after remapping). Maybe there would be other requirements (e.g. conservation in a layer-wise sense, or at least the option for it). @dengwirda, what is possible with PPR?

Regarding monotonicity, @dengwirda would be in a better position than I am to suggest the requirements. But we would certainly want at least user control over how monotonic the reconstruction is.

[cbegeman]

My understanding of PPR is that there is conservation over the layer integral.

[cbegeman]

I'm not sure how detailed it makes sense to get about the details of the remapping operation in the requirements section. Maybe just include layer-wise conservation?

[dengwirda]

@cbegeman, @xylar, yes, the remapping methods should conserve the vertical integrals of whatever is being remapped (sum(h_k * psi_k) for example, where h_k is layer thickness and psi_k is the layer mean value of a given variable).

There are a few different options for order-of-accuracy + monotonicity / slope-limiting in PPR: PLM (2nd-order), PPM (3rd-order), PQM (5th-order), and either monotone or WENO flux / slope limiting.

I agree that there's no need to get into all of these details in the doc / requirements, but just for reference:

• PLM: Piecewise Linear Method
• PPM: Piecewise Parabolic Method
• PQM: Piecewise Quartic Method
• Monotone limiting: No under/overshoots, but perhaps more numerical diffusion than is desired, can limit order-of-accuracy to 2nd-order.
• WENO: Weighted Essentially Non-oscillatory limiting: "eps"-sized under/overshoots, less numerical diffusion, should generally achieve full order-of-accuracy.

[mark-petersen]

@cbegeman thanks for starting this design document. I was trying to read a compiled version. @xylar what's the easiest way to do that for a single page like this? I tried the compass instructions, but the Makefile is missing here (might be in another PR). Looks like rst2html works for a quick html version, but not for equations.

[cbegeman]

@cbegeman thanks for starting this design document. I was trying to read a compiled version. @xylar what's the easiest way to do that for a single page like this? I tried the compass instructions, but the Makefile is missing here (might be in another PR). Looks like rst2html works for a quick html version, but not for equations.

Thanks @mark-petersen. I rebased onto the @xylar's design doc template branch so the Makefile should be there now.

[xylar]

@cbegeman, out of curiosity, what does this do?

[cbegeman]

That's leftover from my personal markup. I'll be removing it. If you add these lines at the top or in a stylesheet (building with rst2html) then you can use :red: to get red text:

.. raw:: html

<style> .red {color:red} </style>

[xylar]

Ah, good to know.

[xylar]

I'm a bit confused by this notation. Since k is the vertical index, shouldn't this be z_1, z_2, etc. instead of z_k^1, z_k^2, etc?

Separate from that confusion, I'm not clear on why the top layer is handled differently than subsequent layers. I would expect one would have a list of layer-interface locations z_k^\textrm{target} that could z-level, z-star or z-tilde for now (and more general in the future). The top layer would not be special.

[xylar]

Ideally, there is a single function that is used for both ALE implementations, with and without VLR.

I think this should be a requirement (maybe to add to the requirement section even). Redundant code would be quite risky in my opinion.

[dengwirda]

@xylar, I believe the surface layer is currently treated differently for coastal wetting + drying (as I understand it, this is the only layer that can flood the coastal zone), but agree that we want the general formulation for VLR --- setting a target grid for all layers.

[xylar]

I see, the coastal wetting-and-drying case is certainly one we should support but maybe not one that needs special attention in this design document, except to say that it would be another target grid in addition to z-star and z-tilde that we would want to support.

[cbegeman]

Since k is the vertical index, shouldn't this be z_1, z_2, etc. instead of z_k^1, z_k^2, etc?

This isn't the best notation choice. In general the superscript is the timestep marker. These numbers are just referring to temporary, intermediate grid choices, but they're maybe better expressed as updates to z_k^{target} prior to h_k actually being reassigned to be consistent with z_k^{target}.

[xylar]

I guess I'm still not following but I agree that the numbers are not a good choice for provisional variables. The POP convention was to use a * or some similar symbol to indicate a provisional variable. Perhaps something like that would be more typical.

[xylar]

At least for purposes of this design, we should probably settle on a single term ("regridding" or "remapping"). My sense is that, in general, "remapping" is more general and includes the potential for unstructured meshes. Either term is correct in this context but I might lean toward the more restrictive term "regridding".

[xylar]

Reading further, I gather that "regridding" here refers to determining the new target grid, while "remapping" is putting the variables onto that target grid. These uses of the terms would not be consistent with my understanding but may have been adopted in the ALE literature.

[dengwirda]

@xylar, yes, I think the "regrid" + "remap" terms are consistent with e.g. HYCOM, MOM6 terminology, etc. See Griffies et al, for example:

• "regrid": select new target grid levels.
• "remap": do the integral interpolation.

I'd personally be happy to use different terminology (others have also found the regrid / remap combination to be a little inscrutable in the past!) Perhaps: "grid selection" and "remap"?

[xylar]

Yes, I'd prefer "grid selection" and "remap".

[cbegeman]

I'm fine with changing the terminology as suggested. However, I think this does raise the question of where we want the layerThickness update to take place because in Griffies et al. this is at least conceptually associated with the regridding step. I think it does make sense to have the remapping module update both prognostic variables and layer thickness and just do grid selection in the regridding module. Sound good to you @xylar and @dengwirda ?

[xylar]

Yes, I think grid selection would merely determine what the target vertical grid is. This would probably be specified in terms of the target layer thickness on cells. But this would not necessarily involve actually updating layerThickness or layerThicknessEdge.

I agree, then, that it would make sense to update layerThickness (and make the call to update layerThicknessEdge) within the remapping step.

So there would not be a "regridding module" but rather something like a "vertical grid" module, I would think. We have something like that already in mode_init but it's nothing to use as a starting point.

[dengwirda]

Yes, I think this makes sense: "grid selection" would determine the position of the new grid levels, and "remap" would update the prognostic variables, thicknesses, etc for the next time step.

[philipwjones]

Keeping thickness and tracer remapping together also helps if you require tracer consistency, so another reason to do it this way.

[xylar]

If this is the only usage, it would be worth discussing with @vanroekel, @milenaveneziani and @mark-petersen whether the MOC itself should be computing the this variable as needed when we're in VLR mode.

One consideration is that the MOC also really needs to be computed on an Eulerian (fixed-in-time) or isopicnal grid. We have treated z-star as sufficiently Eulerian up to now but more complex ALE coordinates would not be. The vertical velocity through ALE surfaces is not what is needed to compute an accurate MOC. Instead, the MOC would want to compute the Eulerian vertical velocity (w = - \nabla \cdot (h_k \mathbf{u}_k) if I'm not mistaken). This work is outside the scope of the current design doc but might be worth mentioning as a .. note::.

[xylar]

I would think all of these should be noted as avenues for exploration in a future project but would not be in the scope of the current design. Instead, I would think our fallback for now if performance isn't great would be that VLR would only be used in contexts where it is required to support a new vertical coordinate, and where the performance hit would be deemed acceptable.

[xylar]

One way of detecting these would be to run a test where you set all the prognostic variables at the previous time level to garbage (-1e34 is a popular choice in MPAS-Ocean) and make sure nothing breaks.

[cbegeman]

Yes, I agree. I wrote something along these lines in the testing section.

[xylar]

Ah, sorry, I didn't get that far.

[xylar]

👍🏻👍🏻👍🏻

And this separate routine should be shared with the non-VLR implementation as well.

[xylar]

@cbegeman, things are looking really good. As far as I'm concerned, it's time to jump into a draft implementation (as I think you already are).

As I noted, one requirement that might need to be made explicit (in an existing requirements section or a new one) is that the function where the target grid is computed should be shared with the non-VLR implementation of the same vertical coordinate (if there is one). So we should only have one place where we define z-star, for example, and it should be shared. Otherwise, there is the potential for bug fixes to be applied to either the VLR version of the non-VLR but not both.

[xylar]

Ah, yep, that's exactly what I was thinking, perfect!

[mark-petersen]

@cbegeman and others, I added a commit with some text, I thought that was easier than marking up with comments. Feel free to revise.

[mark-petersen]

I see this as a choice, so I would prefer being more explicit about false:
config_tracer_vert_advection_method = 'transport' or 'remap'
Then in the code, if (config_tracer_vert_advection_method.eq.'transport') is more clear.

The problem is that this list is now confusing:

&advection
    config_vert_tracer_adv = 'stencil'
    config_vert_tracer_adv_order = 3
    config_horiz_tracer_adv_order = 3
    config_coef_3rd_order = 0.25
    config_monotonic = .true.
/

because the first two flags look like they apply to remapping. I see the first is not even used. So I would do this:

&advection
    config_vert_tracer_advection_method = 'transport' or 'remap' (new)
    config_vert_tracer_adv = 'stencil' --- delete
    config_vert_tracer_adv_order = 3 --> config_vert_tracer_transport_order = 3
    config_vert_tracer_remap_order = 3 (new)
    config_horiz_tracer_adv_order = 3
    config_coef_3rd_order = 0.25
    config_monotonic = .true.
/

I think the new flags belong in the advection list, not the ALE list.

[cbegeman]

@mark-petersen I see. That seems appropriate to me.

[dengwirda]

I think we may actually need the single option to toggle lagrangian-remap on / off, since this'll be changing the vertical advection scheme / terms for tracers, thickness and momentum (rather than just tracers) and it'll be necessary for everything to be done either as remap or transport consistently.

No problem to opt for more explicitness in the options if that's desired, though config_vert_tracer_adv_order = {1,2,3,5} would definitely make direct sense (to me) in the remap case too.

[mark-petersen]

@dengwirda I see what you mean about thickness and momentum. We could instead use

config_vert_advection_method = 'transport' or 'remap'

for tracers, thickness and momentum.

My point was just that if you use a Boolean, it does not communicate what happens for false. A Boolean is appropriate if it is really on/off, but here false implies that you turn on transport instead.

[cbegeman]

I think we may actually need the single option to toggle lagrangian-remap on / off,since this'll be changing the vertical advection scheme / terms for tracers, thickness and momentum (rather than just tracers)

Is this an appropriate fix, make some of the remap-related names more general:

&advection
    config_vert_tracer_advection_method = 'transport' or 'remap' --> config_vert_advection_method = 'transport' or 'remap'
    config_vert_tracer_transport_order = 3
    config_vert_tracer_remap_order = 3 --> config_vert_remap_order = 3
    config_horiz_tracer_adv_order = 3
    config_coef_3rd_order = 0.25
    config_monotonic = .true.

and should we add all of the remap options here:

   config_vert_remap_method = 'pqm'
   config_vert_remap_edge_interp_method = `p5e`
   etc.

[dengwirda]

@cbegeman, yes I think config_vert_advection_method = transport or remap makes sense, although what does everyone think about using flux-form in place of transport?

In terms of the remap options, we may be able to get away with simplifying a little, perhaps something like:

(method = pcm) or (order = 1) --> pcm_method
(method = plm) or (order = 2) --> plm_method
(method = ppm) or (order = 3) --> ppm_method + p3e_method
(method = pqm) or (order = 5) --> pqm_method + p5e_method

It is possible to mix and match some of these options in PPR, but I think I'd expect the combinations above to be the ones typically used in practice.

[cbegeman]

@dengwirda flux-form sounds clear to me.

I'm in favor of simplifying. Are you amenable to cutting config_vert_remap_method and only using config_vert_remap_order since we don't expect folks to mix and match?

If not, do you expect that logic so be applied sequentially so that if order is higher than method, the method is changed; and if method is higher than order, the method stays the same? Given this logic it would be pretty hard to mix and match anyhow, right?

[dengwirda]

Oh, I just meant that I think the edge interp scheme can be inferred, so that ppm+p3e and pqm+p5e always go together, and there's not a separate edge option necessary.

I don't have strong feelings either way in terms of the order vs method question. As you say, I think config_vert_remap_order would work well, and is probably more consistent with the way other MPAS options are set up.

[mark-petersen]

Add to ALE_vertical_grid, config_ALE_target_grid = 'z_initial'

[philipwjones]

While this is a design/req doc, just a reminder that for the implementation, while string choices are fine/encouraged in the input file or namelist, once we're outside the initialization, we should avoid the use of string comparisons and replace with enums (integer params in Fortran) as the choice