Update MALI version to include higher-order advection and time integration #6126
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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.
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!
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.
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.
Previous commit neglected to actually update tracers after calculating tendencies. This corrects that.
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.
Use cellMask_dynamicMargin to calculate advMaskHighOrder. Previously, the mask defining where higher order advection was allowed was incorrectly defined, which led to higher order advection being used at the boundaries where only first-order advection should be used.
Use dynamic ice mask instead of dynamic margin to define advMaskHighOrder. This does not seem to fix issues found with fct thickness advection, but it avoids the stencil potentially using ice-free cells in the mask. Also remove definition of cellMask, edgeMask, and vertexMask in li_mesh, which were causing problems when using cellMask in li_advection_fct_shared_init.
A boundaryCell is one that has at least one empty or non-dynamic neighbor.
Update nTracers as tracers are added in tracer_setup instead of assuming all tracers are present. This makes it easy to index when getting layerThickness back from fct
Allocate arrays after nTracers has been calculated. Previous commit was using nTracers to define dimensions for allocatable arrays before it was calculated.
Pass array of 1s as dummy tracer for fct thickness advection. Passing layerThickness itself as a tracer gave very wrong results.
Fix bug in marking boundaryCell that arose from using li_mask_is_dynamic_ice(jCell) instead of li_mask_is_dynamic_ice(cellMask(jCell)). Also write out number of boundary cells marked on each proc to log files.
Add passiveTracer to help verify advection scheme.
Pass layerThickness as tracer instead of array of 1s. The array of 1s seems to only give first order advection. Using layerThickness as a tracer is currently problematic at the ice edge (causes too much advance, but not clear why), but is successful in the ice interior.
This reverts commit 813ee05. After testing and discussion, it is apparent that including calving within the Runge-Kutta loop is not desirable because it changes the domain between stages of the Runge-Kutta integration.
Calculate fluxAcrossGroundingLineOnCells in advection module for each RK stage and then take weighted average after RK integration instead of calcuating fluxAcrossGroundingLineOnCells on final state. This is more consistent with how fluxAcrossGroundingLine (on edges) is treated.
…alving Add logic to control whether velocity is solved before and/or after calving. This allows both for backward compatibility and for the abilility to use self-consistent velocty, stress, and strain-rate fields when calculating calving.
Make config_update_velocity_before_calving = .false. the default. This is not necessarily the desired behavior for production runs, but it is necessary for baseline comparison when using forward Euler time integration.
…ving. Always update velocity after calving if it was not updated before calving.
… positions Move subglacial hydro and bed topography updates back to their original positions. These were initially moved due to the complexity of interaction with RK integration, but I think subsequent changes have made that move unnecessary.
…velop This merge adds Runge-Kutta time integration to MALI, which is likely necessary when using higher-order advection. Currently, Runge-Kutta integration is used for damage evolution and thickness and tracer advection. Surface and basal mass balances are applied within the Runge-Kutta loop and budgets are calculated at the end of the loop. Front ablation (i.e., calving and face-melting) and bed topography updates are applied via forward Euler outside the Runge-Kutta loop. This could conceivably be updated for full consistency. The Runge-Kutta formulations used here are the two-stage second-order and and three-stage third-order strong stability preserving schemes of Shu and Osher (1988), as well as the four-stage third-order strong stability preserving scheme of Spiteri and Ruuth (2002). See equations 2.47–2.49 in Durran (2010) for an overview. There is one velocity solve per Runge-Kutta stage and an optional solve before calving (more on that below), meaning that the four-stage RK3 scheme is ~25–33% more expensive than the three-stage RK3 scheme for a given time step length. However, the four-stage scheme theoretically allows for an effective CFL fraction of 2.0, while the three-stage scheme is limited to CFL fraction of 1.0. Testing indicates stable results using the four-stage scheme for an effective CFL fraction of 1.8 for a simple grounded slab advecting at uniform speed. Note that when using the four-stage scheme in MALI, the maximum allowable time step is updated, so config_adaptive_timestep_CFL_fraction should still be between 0.0 and 1.0. There is now an optional extra velocity solve between advection and calving in addition to the final velocity solve at the end of the time step, for any choice of time integration scheme. The extra velocity solution prior to calving ensures a self-consistent set of inputs to calving routines; i.e., the velocity, stress, and strain rate fields will be consistent with the geometry used for calving, which is not the case when velocity is solved only after calving. This makes calving more accurate, but is obviously significantly more expensive. The extra velocity solution is enabled using config_update_velocity_before_calving = .true.. It is disabled by default to allow regression tests to pass baseline comparison. Note that the extra velocity solution is not necessary for some calving laws, such as damage threshold calving, so the user should carefully consider this option when setting up simulations. There is currently no check in the code to ensure the combination of calving and extra velocity solve settings are sensible. Also of note is that li_restore_calving_front (if enabled) is called before the extra velocity solution to prevent solving velocities on an extended geometry, while the remaining calving routines are called after the velocity solve. If the volume of dynamic ice is not changed by calving, then the final velocity solution is automatically skipped. References: Durran, Dale R. Numerical Methods for Fluid Dynamics. Vol. 32. Texts in Applied Mathematics. New York, NY: Springer New York, 2010. https://doi.org/10.1007/978-1-4419-6412-0. Shu CW, Osher S (1988) Efficient implementation of essentially nonoscillatory shock-capturing schemes. J Comp Phys 77:439–471 Spiteri RJ, Ruuth SJ (2002) A new class of optimal high-order strong-stability-preserving time discretization methods. SIAM J Numer Anal 40:469–491 * trhille/implement_rk_time_integration: (36 commits) Move subglacial hydro and bed topography updates back to their former positions Always update velocity after calving if it was not updated before calving. Make config_update_velocity_before_calving = .false. the default Add logic to control whether velocity is solved before and/or after calving Average fluxAcrossGroundingLineOnCells in time integration module Revert "Include calving in Runge-Kutta loop" Revert "Calculate layerThickness after calving in RK loop" Calculate layerThickness after calving in RK loop Include calving in Runge-Kutta loop Remove support for config_rk_order = 1 Restore calving front within RK loop Initialize RK weights to large negative values Simple clean-up after code review Update temperature and waterFrac in RK loop, but not enthalpy Update thickness and tracer halos before velocity solve Update stresses and strain rates after calving Reset layerThickness when using config_restore_thickness_after_advection Calculate groundedToFloatingThickness and grounding line flux after RK Fix bug in groundedSfcMassBalApplied after Runge-Kutta loop Some more cleanup ...
| @@ -3,6 +3,8 @@ | |||
|
|
|||
| OBJS = mpas_li_constants.o \ | |||
| mpas_li_mask.o \ | |||
| mpas_li_mesh.o \ | |||
There was a problem hiding this comment.
Noting another place where E3SM build system may need updating
|
@jonbob , this PR will also need updates to the namelist generation files in E3SM because we've added some new namelist options. Let me know if you are able to do that or prefer that I do it. |
|
Thanks for pointing that out, @matthewhoffman -- I'm happy to do that |
|
@matthewhoffman -- I did update the bld files to match the Registry changes, and I think they're mostly straightforward. The only part I'm not sure of is in We're adding some new advection configs and it's not clear to me if they need to be wrapped in logic for MALI running in its dynamic mode like we do for config_thickness_advection and config_tracer_advection: |
@jonbob, the new advection settings should not have any effect on the default MALI dynamic mode. |
When the velocity solver fails, the internal dissipation field it returns can have large artifacts in it that lead to unphysical melt rates being applied. This commit skips updating the heatDissipation field when the solver fails so that the previous timestep's value is retained.
…into MALI-Dev/develop This merge fixex dissipation when solver fails. When the velocity solver fails, the internal dissipation field it returns can have large artifacts in it that lead to unphysical melt rates being applied. This merge skips updating the heatDissipation field when the solver fails so that the previous timestep's value is retained.
matthewhoffman
force-pushed
the
mali/develop-231220
branch
from
fe8d760 to
2c646ac
Compare
matthewhoffman
force-pushed
the
mali/develop-231220
branch
from
2c646ac to
9af3b1e
Compare
matthewhoffman
approved these changes
Jan 16, 2024
There was a problem hiding this comment.
@trhille , thanks for putting this together. This is good timing with our push starting up to get MALI runs more operational in E3SM.
@jonbob , I rebased this to bring in a bug fix we just added to MALI that will be useful to have right away in E3SM. I also added a new commit that adds new modules to the landice.cmake file used for building MALI in E3SM. This looks good to me, pending final testing in E3SM.
|
Passes BFB with expected NML DIFFs for:
BFB and no NML differences for:
ready to be merged |
jonbob
added a commit
that referenced
this pull request
Jan 24, 2024
Update MALI version to include higher-order advection and time integration This merge includes numerous updates to MALI The major updates include: * Higher-order, flux-corrected thickness and tracer transport. * Three strong stability-preserving Runge-Kutta time integrations schemes (a two-stage, second-order scheme; three- and four-stage third-order schemes) Minor updates include: * Earlier update of clock within time step * Prevent unrealistic calving into holes within an ice shelf * Change default start and stop times to avoid unpredictable I/O behavior * Higher-order grounding-line flux calculation * Several small bug fixes [NML] - only for cases with MALI [BFB]
|
merged to next |
|
merged to master and expected NML DIFFs blessed |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This merge includes numerous updates to MALI
The major updates include:
• Higher-order, flux-corrected thickness and tracer transport.
• Three strong stability-preserving Runge-Kutta time integrations schemes (a two-stage, second-order scheme; three- and four-stage third-order schemes)
Minor updates include:
• Earlier update of clock within time step
• Prevent unrealistic calving into holes within an ice shelf
• Change default start and stop times to avoid unpredictable I/O behavior
• Higher-order grounding-line flux calculation
• Several small bug fixes
[NML] - only for cases with MALI
[BFB]