Skip to content

Add FastIsostasy#183

Open
matthewhoffman wants to merge 8 commits into
developfrom
matthewhoffman/mali/fastisostasy
Open

Add FastIsostasy#183
matthewhoffman wants to merge 8 commits into
developfrom
matthewhoffman/mali/fastisostasy

Conversation

@matthewhoffman

@matthewhoffman matthewhoffman commented Jun 13, 2026

Copy link
Copy Markdown

This PR adds FastIsostasy as an alternate GIA model to the existing global 1d Sea Level Model. The implementation makes use of a lot of the existing coupling code and builds off of FastIsostasy documentation on coupling here:
https://github.com/palma-ice/FastIsostasy/blob/main/README.md#coupling-fastisostasy-to-your-favorite-ice-sheet-model

Note: Currently based on #180 and needs to be rebased after that is merged.

@matthewhoffman

matthewhoffman commented Jun 13, 2026

Copy link
Copy Markdown
Author

Status: Code compiles with FASTISOSTASY=true but has not been run or tested in any way. Code was written by Copilot w/ GPT 5.3-Codex using these instructions:

Plan instructions

Plan: Couple MALI with FastIsostasy (Fortran)

TL;DR: Use the Fortran version of FastIsostasy (palma-ice/FastIsostasy) — not the Julia version — and follow the existing SLM coupling pattern already in mpas_li_bedtopo.F. The Fortran library has a clean 4-call API and the SLM code provides a near-complete template for the MPI gather/interpolate/solve/scatter workflow.


Key Discovery: A Fortran version already exists

FastIsostasy includes a pure Fortran implementation specifically designed for coupling with ice-sheet models. Its API is:

call isos_init(isos, path_par, "isos", nx, ny, dx, dy)
call isos_init_ref(isos, z_bed_ref, H_ice_ref)
call isos_init_state(isos, z_bed, H_ice, time, bsl)
! In time loop:
call isos_update(isos, H_ice, time, bsl)
z_bed = isos%out%z_bed

The docs for it include specific guidance for coupling to an ice-sheet model: https://github.com/palma-ice/FastIsostasy#coupling-fastisostasy-to-your-favorite-ice-sheet-model


Why this is tractable

MALI already has a complete coupling template in mpas_li_bedtopo.F via the SLM (sea-level model) coupling. That code handles:

  • Alarm-based coupling intervals
  • MPI gather of ice thickness from all ranks to rank 0
  • Interpolation from unstructured MPAS mesh → regular grid (sparse matrix multiply using weight files)
  • External solver call on rank 0
  • Interpolation back to MPAS mesh
  • MPI scatter of bed topography changes
  • Halo exchange + geometry/mask update

The FastIsostasy coupling would replicate this pattern almost exactly, substituting FastIsostasy's isos_update() for the SLM solver calls.


Phases

Phase 1: Build System

  • Add palma-ice/FastIsostasy as a submodule or vendored source
  • Add USE_FASTISOSTASY conditional compilation flag (mirrors USE_SEALEVELMODEL)
  • Link against libisostasy.a + FFTW dependency

Phase 2: Registry & Config

  • New config_uplift_method = 'fastisostasy' option
  • Config params: coupling interval, namelist path, weight file paths

Phase 3: Core Couplingmost complex, but heavily templated

  • fastisostasy_init(): setup alarm, read weights, MPI layout, call isos_init/isos_init_state
  • fastisostasy_solve(): gather → interpolate → isos_update() → interpolate back → scatter → apply
  • New branch in li_bedtopo_solve() dispatch

Phase 4: Grid/Weight Files

  • Generate ESMF/TempestRemap interpolation weights between MPAS mesh and FastIsostasy's regular Cartesian grid
  • FastIsostasy needs a square domain with ≥512 km padding beyond ice extent

Phase 5: Restart

  • FastIsostasy has built-in isos_restart_write/read; wire into MALI checkpoint

Phase 6: Testing

  • Idealized dome test with known analytic GIA solution
  • MPI correctness (1/4/16 procs)
  • Bit-for-bit restart verification

Relevant Files

  • mpas_li_bedtopo.F — primary file to modify (SLM coupling template is here)
  • Registry.xml — add config options
  • Makefile — build system (has USE_SEALEVELMODEL pattern to copy)
  • mpas_li_time_integration_fe_rk.F — calls li_bedtopo_solve() at line 634 (no changes needed here)

Main Challenge: Grid Interpolation

MALI uses an unstructured MPAS mesh; FastIsostasy uses a regular Cartesian (stereographic) grid. The SLM coupling already solves this exact problem using pre-computed sparse interpolation weight files (generated by ESMF_RegridWeightGen or TempestRemap). The same tooling applies — you just need to generate weights for the FastIsostasy target grid rather than the SLM's Gauss-Legendre grid. This will be done outside of MALI in a separate preprocessing step and does not need to be implemented here. However, the MALI implementation should assume a similar format for mapping files as for the existing SLM.

Build and test loop

  • I want the work divided into logical commits with descriptive commit messages.

  • To get a usable environment for compiling MALI, do:
    source /global/cfs/cdirs/fanssie/users/hoffman2/compass/ismip7-devs//load_compass_pm-cpu_gnu_mpich.sh

  • Before each commit, make sure the code builds:
    cd /global/cfs/cdirs/fanssie/users/hoffman2/compass/ismip7-devs/MALI-Dev/components/mpas-albany-landice; make clean; make -j 4 gnu-cray DEBUG=true

  • If code compiles, you may commit and move on to next step.

  • If you run into any environment problems or have any questions, stop and ask for more instructions. Otherwise continue working through the implementation steps until complete.

I made a cursory review of the implementation. The coupling code looks correct in that it mirrors the coupling code for the SLM heavily, as intended, but I have not done a line by line review. The Makefile/build changes look reasonable, but some cleanup is needed to generalize it to any machine. Note that the FastIsostasy build is really slow for two reasons: 1. it uses config which is generally slow and 2. the FFTW library it depends on is very slow to compile and we may want to support optionally pointing to a pre-build lib of it.

Next steps are to build mapping files, set up a FI namelist, and attempt a run.

@matthewhoffman matthewhoffman added the in progress Still being worked on, don't merge yet! label Jun 13, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

in progress Still being worked on, don't merge yet! testing needed

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant