Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add flux-corrected transport scheme to MALI #70

Merged
merged 46 commits into from
Sep 6, 2023

Conversation

trhille
Copy link

@trhille trhille commented Mar 7, 2023

This merge adds a flux-correct transport (FCT) scheme to MALI for thickness and tracer advection, ported over with MALI-relevant modification from MPAS-Ocean's routines, which are based on Skamarock and Gassmann (2011) (https://doi.org/10.1175/MWR-D-10-05056.1). This uses a blend of 3rd and 4th order fluxes to achieve monotonicity. The FCT routine is only used for tracers in MPAS-Ocean, whereas here it is modified for use with both thickness and tracers. The user can specify 2nd, 3rd, or 4th order advection with the config_horiz_tracer_adv_order option, but only config_horiz_tracer_adv_order = 3 with 0 < config_advection_coef_3rd_order < 1 is truly FCT. The config_advection_coef_3rd_order option specifies the blend between 3rd and 4th order fluxes used in the flux correction. config_advection_coef_3rd_order = 1.0 is purely 3rd order, while config_advection_coef_3rd_order = 0.0 is purely 4th order. The default value of 0.25 is taken from the MPAS-Ocean default and may not be appropriate for all situations. Note that all higher-order advection must reduce to 1st order at the boundaries.

This also adds a new variable passiveTracer2d, that can be used to verify advection schemes.

Currently supported combinations of thickness and tracer advection with fct include:

  1. config_thickness_advection = 'fo'; config_tracer_advection = 'fct'
  2. config_thickness_advection = 'fct'; config_tracer_advection = 'fct'
  3. config_thickness_advection = 'fct'; config_tracer_advection = 'none'
    FCT tracer advection with no thickness advection and FCT thickness advection with FO tracer advection could be added, but we do not currently anticipate using them. Therefore, we have left them out of this PR as they would add unnecessary complexity to the code.

When used with forward Euler time integration, FCT requires a severely reduced time step relative to first order advection. Testing shows that CFL fraction on the order of 0.1 is likely sufficiently small, but this is probably case-dependent. Once Runge-Kutta time integration is operational, it is recommended to use that instead of forward Euler.

At the time of merging, there is a very slight conservation error (<0.1% in the grounded slab test shown below) for tracers when using FCT for both thickness and tracer advection, while mass is fully conserved. Tracers are conserved when using FO thickness advection with FCT tracer advection. The convergence order has not yet been verified, but testing is under way using a new mesh convergence COMPASS test: MPAS-Dev/compass#690.

Add an option to use the flux-corrected transport advection scheme
provided in MPAS framework, and make calls to those modules from
mpas_li_advection. Compiles, but not yet tested. No expectation that
it will run yet.
We have decided to use modified versions of the routines for
higher-order advection from MPAS-Ocean. This commit copies over the
modules from the ocean model, but does not modify or call them yet.
@trhille trhille added enhancement New feature or request in progress Still being worked on, don't merge yet! labels Mar 7, 2023
Update Makefile to include fct modules. Currently does not compile
because of undefined variables.
Initialize fct from tracer_setup. Also rename public routines to
start with "li_"
Add li_mesh module that is called during init to make all mesh
parameters public. This is copied from mpas-ocean, with modifications
appropriate for MALI.
Comment out irrelevant parameters and routine in li_mesh that are
specific to MPAS-Ocean and caused build errors.
Add li_config module so other modules can get all configs with
'use li_config'
Also general cleanup to get rid of build errors stemming from
these modules.
Some fields need to be revisited for accuracy, but these changes
fix build errors.
Remove call to tracer_advection_vert_flx from li_tracer_advection_fct_tend.
We will probably want to add this back, but for now just trying
to get this to compile without vertical fct. And now it compiles!
@trhille
Copy link
Author

trhille commented Mar 9, 2023

Code finally compiles as of 414eac6!

Differences between ocean and land ice mesh variables were causing
segmentation faults at run time. This clean-up commit allows full_integration
to run and pass, which does not test the fct code at all, but shows
that li_mesh is now operational.
@trhille
Copy link
Author

trhille commented Mar 10, 2023

full_integration passes as of 83f5168, which does not test the fct at all, but indicates that the new li_mesh and li_config modules are operational.

Fix issues with public parameter allocations. This now allows the code
to run using config_tracer_advection='fct' without invalid memory
reference errors, although results are not yet valid.
Clean up li_advection_fct_tend to get rid of runtime errors. All
vertical advection code is currently commented out and not operational.
This runs with config_thickness_advection='fo' and config_tracer_advection='fct',
but it currently seems to give the same results for the humboldt 3km test
case in full_integration as using 'fo' for both thickness and tracer advection.
Need to look further into that.
Currently fct advection for thickness is not supported, but fo
thickness advection wtih fct tracer advection is operational, although
not validated.
@trhille
Copy link
Author

trhille commented Mar 14, 2023

Testing on a364f7a shows that the code will run with config_thickness_advection = 'fo' and config_tracer_advection = 'fct', but from examining damage output, tracers don't seem to actually be advecting.

Previous commit neglected to actually update tracers after calculating
tendencies. This corrects that.
@trhille
Copy link
Author

trhille commented Mar 16, 2023

As of e85d3fd, fct tracer advection is at least doing something, although is not yet verified. Here are three MISMIP runs after 100 years with velo solver turned off, no sub-shelf melt, calculate-damage = true. All use 'fo' thickness advection. Left: no tracer advection; Middle: fo tracer advection; right: fct tracer advection:
image
White contour is damage = 0.4 from the 'fo' case.

Pass layerThicknessOld instead of layerThickness to fct
tracer advection routine because the fct routine calculates a
provisional advected layerThickness based on the normalThicknessFlux,
which is layerThicknessEdge (before advection) * layerNormalVelocity.
Add support for fct thickness advection by passing layerThickness
to li_tracer_advection_fct_tend as a tracer. Note that this seems
to require small time steps to avoid terracing effects in the
thickness field.
@trhille
Copy link
Author

trhille commented Mar 16, 2023

As of 8734303, fct thickness advection is operational, and seems to be working correctly. One issue is that it seems to require a very small CFL fraction to avoid terracing effects.
From left to right: fo thickness advection with CFL fraction = 0.8; fct thickness advection w/ CFL fraction 0.8; fct thickness advection w/ CFL fraction 0.1:
image

CFL fraction of 0.1 (~40 day time steps) seems to be sufficiently low here. I tested another run w/ a CFL fraction of 0.01 (~4 days) and the differences at year 100 are generally <1 m:
image

What we are seeings seems to be a known issue w/ FCT: "Another infamous byproduct of FCT manifests itself in distortions of a smooth profile. This phenomenon is known as terracing and represents ‘an integrated, nonlinear effect of residual phase errors’ [64] or, loosely speaking, ‘the ghosts of departed ripples’ [7]. A particularly severe form of terracing is caused by the linear instability of the high-order scheme. For this reason, we do not recommend the use of the forward Euler time-stepping (θ = 0) even though the flux-corrected scheme proves positivity-preserving under the CFL-like time step restriction (62)."
(found on page 178 of https://link.springer.com/content/pdf/10.1007/978-94-007-4038-9.pdf)

@trhille
Copy link
Author

trhille commented Mar 16, 2023

Adding to comment above: CFL fraction of 0.4 also drastically reduces the terracing and gives a thickness field visually indistinguishable from CFL fractions of 0.1 and 0.01, so that's great news! Differences between runs w/ CFL fractions of 0.4 and 0.1:
image

@trhille
Copy link
Author

trhille commented Mar 17, 2023

Just for comparison, here are the thickness differences at year 100 when using FO thickness advection with CFL of 0.8 vs 0.1:
image
So the FCT solution is definitely more sensitive to the time step than the FO solution, even when the time step is quite small.

@trhille
Copy link
Author

trhille commented Mar 17, 2023

Testing on Thwaites mesh indicates issue with FCT thickness advection during advance, which is not surprising. It turns out that boundaryCell is not being calculated, and thus neither is advMaskHighOrder. This leads to very high velocities and strange geometries in areas where ice advances beyond the initial extent.

First order flux needs to be masked when subtracted from higher
order flux, because they don't have the same valid spatial extent.
This seems to largely fix the issue with thickness ramps growing
at the ice margin in tests on the Thwaites domain. It also leads to
a bit more diffusion of the ice mask in the slotted cylinder test,
but it also gets rid of the noisy artifacts in the thickness field,
so that's great!
@trhille
Copy link
Author

trhille commented Apr 24, 2023

I have confirmed that results are BFB for restarts and decomposition tests (64 and 128 cores) using the slotted cylinder test case. I should also put together a version of full_integration that uses fct.
Here are the restart test results, with the full run on top, a restart run branching from year 10 in the middle and the difference at year 25 on bottom. This uses 3rd order, beta = 0.5, and CFL fraction = 0.25 (hence the ugliness)

image

And here is the same for the decomp test. Top is 128 procs, middle is 64 procs, bottom is difference:
image

@trhille trhille removed the in progress Still being worked on, don't merge yet! label Apr 24, 2023
@trhille
Copy link
Author

trhille commented Apr 24, 2023

Here are results at year 25 with CFL fractions of 0.25 (top), 0.1 (middle) and 0.05 (bottom):
image

@trhille
Copy link
Author

trhille commented Aug 1, 2023

Uh oh, it appears that passiveTracer is not being conserved in test runs with fct. I ran the notched slab test that I was using above for 100 years and summed up passiveTracer at each time step: ncap2 -s 'tracer_sum = (passiveTracer * thickness * areaCell).sum($nCells)' output_all_timesteps.nc tmp.nc.
fct:
image
fo conserves just fine:
image

Both runs do conserve mass, however. Blue is fo, orange is fct:
image

Update:
I rewound the branch to 4c87809 and the issue is still there, but it is less egregious and drifts the opposite way:
image

Fix a bug in fct tracer advection that caused tracers to not be conserved
due to an inconsistent treatment of the thickness flux. The previous
implementation inadvertently used first-order upwind thickness flux
to calculate higher-order tracer flux, which led to loss of conservation
when both thickness and tracers used fct. Note that this still does
not conserve tracers to machine precision, even though mass is conserved
to machine precision. But testing shows the tracer volume conserved
to within about 0.02%.
@trhille
Copy link
Author

trhille commented Aug 3, 2023

8d2f023 fixed the tracer conservation issue detailed above! It still doesn't conserve to machine precision, but a 100 year test run on the grounded slab test case showed about 0.06% drift by the end of the simulation. That seems pretty close!

Compare this with the drift above
image

@matthewhoffman
Copy link

@dengwirda, see note about slight non-conservation of tracers in Trevor's current implementation (#70 (comment)). Should we be expecting machine precision conservation? If so, do you think we should worry that we're not achieving that given how small the non-conservation is?

I suppose we would be introducing energy conservation errors, given that we will be advecting temperature with this scheme.

Copy link

@matthewhoffman matthewhoffman left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@trhille , I've made it through these changes and have a number of requests and questions to address, but this seems very close to being ready to merge. Thanks for you hard work and what turned out to be a much more involved project than initially anticipated!

components/mpas-albany-landice/src/Registry.xml Outdated Show resolved Hide resolved
! Determine bounds on tracer (tracerMin and tracerMax) from
! surrounding cells for later limiting.

#ifdef MPAS_OPENACC

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove accelerator stuff?

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We decided to save removing accelerator stuff for a separate PR.

! This does conserve mass:
layerThickness(:,:) = layerThickness(:,:) + tend(nTracers,:,:) * dt
! Reset tracers after thickness advection for tracer advection
! TODO: determine whether we need to treat other tracer tendencies

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@trhille , have you looked into this question?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, it doesn't work for tracers, for some reason.

call li_tracer_advection_fct_tend(&
tend, advectedTracers, layerThicknessOld, &
layerThicknessEdge * layerNormalVelocity, 0 * normalVelocity, dt, &
nTracers, computeBudgets=.false.)!{{{
nTracers, activeTracerHorizontalAdvectionEdgeFlux, computeBudgets=.false.)!{{{

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We may want to move activeTracerHorizontalAdvectionEdgeFlux to Registry so it could be output given its importance in the algorithm. But we'd want to refactor things to avoid it being nTracers big first. So let's wait until we encounter a need to visualize it before worrying about that.

Add explanatory comments, remove commented code, clarify desciriptions
in Registry.xml
Throw error forr fct thickness with fo tracer, as this combination
is currently not supported.
Make activeTracerHorizontalAdvectionEdgeFlux an optional argument
when calling fct code.
This fixes an error introduced in the previous commit that prevented
code from compiling.
Fix a bug that disabled thickness advection when tracer advection
is turned off.
Only pass layer thickness tracer instead of all tracers for thickness
advection step, which is the first call to fct when using fct for
both thickness and tracers.
Change passiveTracer to passiveTracer2d
Copy link

@matthewhoffman matthewhoffman left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@trhille , thanks for all the clean up on this PR. I would like to proceed with merging the version you have now and resolve additional issues related to conservations and convergence in follow up PRs as we learn more.

However, it would be nice to first add the logic to support FCT thickness advection with no tracer advection. That would be useful for the horizontal_advection convergence test using the thickness bump and for the halfar convergence test. I think disabling tracer advection would speed those tests up considerably.

One solution would be to fold it into this option:

         elseif ((trim(config_thickness_advection) .eq. 'fct') .and. &
                  trim(config_tracer_advection) .eq. 'fct') then

and then only do the tracer bit if tracer advection is set to fct.

Enable fct thickness advection without tracer advection and throw error
if using unsupported  combination of thickness and tracer advection. Also
initialize some previously uninitialized allocatable arrays.
@matthewhoffman
Copy link

matthewhoffman commented Sep 6, 2023

Testing

@trhille , thanks for making the final (more complicated than expected) adjustment. I tested one more time and all tests pass:

Test Runtimes:
00:19 PASS landice_dome_2000m_sia_restart_test
00:19 PASS landice_dome_2000m_sia_decomposition_test
00:16 PASS landice_dome_variable_resolution_sia_restart_test
00:06 PASS landice_dome_variable_resolution_sia_decomposition_test
00:26 PASS landice_enthalpy_benchmark_A
00:24 PASS landice_eismint2_decomposition_test
00:23 PASS landice_eismint2_enthalpy_decomposition_test
00:29 PASS landice_eismint2_restart_test
00:28 PASS landice_eismint2_enthalpy_restart_test
00:22 PASS landice_greenland_sia_restart_test
00:12 PASS landice_greenland_sia_decomposition_test
00:20 PASS landice_hydro_radial_restart_test
00:12 PASS landice_hydro_radial_decomposition_test
00:44 PASS landice_humboldt_mesh-3km_decomposition_test_velo-none_calving-none_subglacialhydro
00:30 PASS landice_humboldt_mesh-3km_restart_test_velo-none_calving-none_subglacialhydro
00:26 PASS landice_dome_2000m_fo_decomposition_test
00:30 PASS landice_dome_2000m_fo_restart_test
00:11 PASS landice_dome_variable_resolution_fo_decomposition_test
00:16 PASS landice_dome_variable_resolution_fo_restart_test
00:18 PASS landice_circular_shelf_decomposition_test
00:44 PASS landice_greenland_fo_decomposition_test
00:44 PASS landice_greenland_fo_restart_test
00:41 PASS landice_thwaites_fo_decomposition_test
00:32 PASS landice_thwaites_fo_restart_test
00:15 PASS landice_thwaites_fo-depthInt_decomposition_test
00:22 PASS landice_thwaites_fo-depthInt_restart_test
00:58 PASS landice_humboldt_mesh-3km_restart_test_velo-fo_calving-von_mises_stress_damage-threshold_faceMelting
00:32 PASS landice_humboldt_mesh-3km_restart_test_velo-fo-depthInt_calving-von_mises_stress_damage-threshold_faceMelting
Total runtime 12:04
PASS: All passed successfully!

I will also test the FCT tests you have in MPAS-Dev/compass#650

@matthewhoffman
Copy link

Testing, Part 2

The FCT suite all passes. I did not look at any of the output, but they all run and pass validation. For the sake of merging this PR as a checkpoint, I think that is sufficient.

Test Runtimes:
00:19 PASS landice_dome_2000m_fo_fct_decomposition_test
00:35 PASS landice_dome_2000m_fo_fct_restart_test
00:16 PASS landice_dome_variable_resolution_fo_fct_decomposition_test
00:14 PASS landice_dome_variable_resolution_fo_fct_restart_test
00:42 PASS landice_greenland_fo_fct_decomposition_test
00:47 PASS landice_greenland_fo_fct_restart_test
00:38 PASS landice_thwaites_fct_decomposition_test
00:34 PASS landice_thwaites_fct_restart_test
00:35 PASS landice_humboldt_mesh-3km_restart_test_velo-fo_advec-fct_calving-von_mises_stress_damage-threshold_faceMelting
Total runtime 04:46
PASS: All passed successfully!

Copy link

@matthewhoffman matthewhoffman left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The long-awaited for day is upon us! Awesome work @trhille !

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants