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.

Sign up for GitHub

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 restart option when the SLM is coupled to MALI #99

Merged
merged 11 commits into from
Mar 1, 2024
Merged

Add restart option when the SLM is coupled to MALI #99

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
35bed8a
Add restart option when the SLM is coupled to MALI
hollyhan Jan 11, 2024
087aace
Improve synchronization of timesteps between MALI and SLM
matthewhoffman Feb 21, 2024
55e579b
Improve error handling, correct other usage of config_uplift_method
matthewhoffman Feb 21, 2024
32ff9f8
Update checks using interval division
matthewhoffman Feb 21, 2024
869f944
Adjust check if adaptive dt is on or not
matthewhoffman Feb 21, 2024
a4d0d68
Update restart check to also use time interval division
matthewhoffman Feb 21, 2024
a3c4337
Add missing arguments to log write statement
matthewhoffman Feb 21, 2024
abaa6ca
Add missing error flag so model actually dies when error occurs
matthewhoffman Feb 24, 2024
6bb3000
Allow restarts at any interval when using SLM
matthewhoffman Feb 24, 2024
be5f341
Add addl info on restart about the calculated time since last SLM call
matthewhoffman Feb 24, 2024
052b15d
Don't call SLM on init of a restart
matthewhoffman Feb 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions components/mpas-albany-landice/src/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -149,6 +149,10 @@
description="Time interval at which the sea-level model is called by MALI. The interval has to be an even multiple of the option 'config_adaptive_timestep_force_interval"
possible_values="Any time interval of the format 'YYYY-MM-DD_HH:MM:SS'"
/>
<nml_option name="config_slm_timestep_restart" type="integer" default_value="0" units="unitless"
description="number of timesteps the sea-level model moved forward since the start of simulation."
possible_values="Any integer value calculated from dividing the difference between simulation start time and restart time by MALI-SLM coupling interval. This poses the restriction on the restart time to be a integer multiple of MALI-SLM coupling interval."
/>
<nml_option name="config_MALI_to_SLM_weights_file" type="character" default_value="mpas_to_grid.nc" units="unitless"
description="File containing the interpolation weights for regridding from MPAS mesh to the Gaussian grid used by the Sea Level Model."
possible_values="Any file name string"
Expand Down
193 changes: 119 additions & 74 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_bedtopo.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,7 @@ subroutine li_bedtopo_init(domain, err)

call mpas_pool_get_config(liConfigs, 'config_uplift_method', config_uplift_method)
if (trim(config_uplift_method)=='sealevelmodel') then
! initialize the 1D sea-level model
! initialize the 1D sea-level model if fresh start
Copy link

Choose a reason for hiding this comment

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

I don't think this change makes sense, because this call is happening whether this is a restart or cold start.

Suggested change
! initialize the 1D sea-level model if fresh start
! initialize the 1D sea-level model

call slmodel_init(domain, err)
endif

Expand Down Expand Up @@ -362,6 +362,7 @@ end subroutine li_bedtopo_finalize
subroutine slmodel_init(domain, err)

#ifdef USE_SEALEVELMODEL
use mpas_timekeeping
use sl_model_mod !< this is part of the SLM code
use sl_io_mod !< this is part of the SLM code
use user_specs_mod, only: nglv !< this is part of the SLM code
Expand All @@ -386,16 +387,21 @@ subroutine slmodel_init(domain, err)
!-----------------------------------------------------------------

#ifdef USE_SEALEVELMODEL
logical, pointer :: config_do_restart
integer, pointer :: config_slm_timestep_restart
character (len=StrKIND), pointer :: config_slm_coupling_interval
character (len=StrKIND), pointer :: xtime, simulationStartTime
type (mpas_pool_type), pointer :: meshPool !< mesh information
type (mpas_pool_type), pointer :: geometryPool
type (MPAS_Time_Type) :: currTime !< current time as time type
real (kind=RKIND), dimension(:), pointer :: thickness, bedTopography
real (kind=RKIND), dimension(:), allocatable :: meshMask
real (kind=RKIND), dimension(nglv,2*nglv) :: ismIceload, ismBedtopo, ismMask
real (kind=RKIND), dimension(nglv*2*nglv) :: thicknessSLgrid1D
real (kind=RKIND), dimension(nglv*2*nglv) :: bedtopoSLgrid1D
real (kind=RKIND), dimension(nglv*2*nglv) :: maskSLgrid1D
integer :: slm_coupling_interval
integer :: slm_coupling_interval, slmRestartYear, slmTimeStep_restart
integer :: currYear, simulationStartYear
integer :: err_tmp
integer :: unit_num_slm ! SLM variable
integer :: itersl, dtime ! SLM variable
Expand Down Expand Up @@ -431,11 +437,15 @@ subroutine slmodel_init(domain, err)

! Gather nCellsOwned
call MPI_GATHER( nCellsOwned, 1, MPI_INTEGER, nCellsPerProc, 1, MPI_INTEGER, &
0, domain % dminfo % comm, err_tmp)
0, domain % dminfo % comm, err_tmp)
err = ior(err, err_tmp)

! Set Displacement variable for GATHERV command
if (curProc.eq.0) then
! First, set SLM unit number to the MALI output log file unit
unit_num_slm = domain % logInfo % outputLog % unitNum
call sl_set_unit_num(unit_num_slm)

! Set Displacement variable for GATHERV command
nCellsGlobal = sum(nCellsPerProc)
allocate(indexToCellIDGathered(nCellsGlobal))
nCellsDisplacement(1) = 0
Expand All @@ -451,89 +461,124 @@ subroutine slmodel_init(domain, err)

! Gather indexToCellID
call MPI_GATHERV( indexToCellID, nCellsOwned, MPI_INTEGER, indexToCellIDGathered, &
nCellsPerProc, nCellsDisplacement, MPI_INTEGER, 0, domain % dminfo % comm, err_tmp)
nCellsPerProc, nCellsDisplacement, MPI_INTEGER, 0, domain % dminfo % comm, err_tmp)
err = ior(err, err_tmp)

call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool)
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
! check if the run is a restart or not
call mpas_pool_get_config(liConfigs, 'config_do_restart', config_do_restart)
if (config_do_restart) then
! find the right slmTimeStep
call mpas_pool_get_config(liConfigs,'config_slm_coupling_interval', config_slm_coupling_interval)
read(config_slm_coupling_interval(1:4),*) slm_coupling_interval
call mpas_pool_get_config(liConfigs,'config_slm_timestep_restart', config_slm_timestep_restart)
Copy link

Choose a reason for hiding this comment

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

So is the user required to manually set config_slm_timestep_restart before doing a restart? What purpose does this serve if the code is able to calculate slmTimeStep_restart on its own? I would prefer to not require an additional namelist change to do a restart with the SLM if it's not necessary. I think it's fair to assume that a user hasn't has not messed around with any output files before a restart.

Copy link
Author

Choose a reason for hiding this comment

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

@matthewhoffman
The purpose is to make sure the restart time of MALI (taken either from config_start_time or restart_timestamp is setup correctly) is set correctly such that the restart year is a correct multiple of the year from config_slm_coupling_interval. Note that slmTimeStep_restart is calculated using currYear, which comes from MALI. This wouldn't be needed if I'd used the mpas alarm functionality and make it automatic in figuring out what the proper restart year should be for the SLM, which could be done in the next PR at some point.

Copy link

Choose a reason for hiding this comment

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

I'm still confused. If you are checking to make sure these are identical, and then using the value if they are identical, can't we just use the internal calculation and not require the user to set this?

We could instead add a check that the slm coupling interval matches the restart interval. That is the requirement you are worried about, if I understand this correctly.

slmTimeStep = config_slm_timestep_restart

! check if the timestep is setup properly
! get the simulation start time
call mpas_pool_get_array(meshPool, 'simulationStartTime', simulationStartTime)
read(simulationStartTime(1:4),*) simulationStartYear

! get the current (restart) time
call mpas_pool_get_array(meshPool, 'xtime', xtime)
read(xtime(1:4),*) currYear

! take the difference and divide by the coupling interval
slmTimeStep_restart = ( currYear - simulationStartYear ) / slm_coupling_interval

! if the result is not equal to the defined slmTimeStep, raise an error.
if (slmTimeStep /= slmTimeStep_restart) then
slmRestartYear = ( slmTimeStep * slm_coupling_interval ) + simulationStartYear
call mpas_log_write("Restart year from SLM restart timestep and MALI restart time &
are inconsistent", MPAS_LOG_ERR)
call mpas_log_write("Problem: Restart year based on `config_slm_timestep_restart` &
is $i, but MALI restart year is set to $i.", &
intArgs=(/slmRestartYear, currYear/))
err = ior(err,1)
else
call mpas_log_write("Calling the SLM. SLM timestep $i", intArgs=(/slmTimeStep/))
call slmodel_solve(slmTimeStep, domain)
Copy link

Choose a reason for hiding this comment

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

@hollyhan , do you know if this call to slmodel_solve on init on a restart is necessary? The time level on init is a time level that was already fully run on the previous run, so I don't think we need to run the SLM again. It wouldn't hurt, but I think it can be skipped. Unless you are aware of a reason it is necessary.

To get restarts to work when using restart intervals shorter than the coupling interval, this call would need to be deleted or only made if the restart time happened to coincide with the coupling interval.

endif

if (curProc.eq.0) then
matthewhoffman marked this conversation as resolved.
Show resolved Hide resolved
allocate(globalArrayThickness(nCellsGlobal), gatheredArrayThickness(nCellsGlobal))
allocate(globalArrayBedTopography(nCellsGlobal), gatheredArrayBedTopography(nCellsGlobal))
allocate(meshMask(nCellsGlobal))
ismIceload(:,:) = 0.0
ismBedtopo(:,:) = 0.0
ismMask(:,:) = 0.0
bedtopoSLgrid1D(:) = 0.0
thicknessSLgrid1D(:) = 0.0
maskSLgrid1D(:) = 0.0
else
! Intel requires these be allocated even though they are not meaningful on the non-destination procs
allocate(globalArrayThickness(1), gatheredArrayThickness(1))
allocate(globalArrayBedTopography(1), gatheredArrayBedTopography(1))
allocate(meshMask(1))
endif

! Gather only the nCellsOwned from thickness and bedtopo (does not include Halos)
call MPI_GATHERV((thickness*real(li_mask_is_grounded_ice_int(cellMask),RKIND)), &
nCellsOwned, MPI_DOUBLE, gatheredArrayThickness, nCellsPerProc, &
nCellsDisplacement, MPI_DOUBLE, 0, domain % dminfo % comm, err_tmp)
err = ior(err, err_tmp)
call MPI_GATHERV(bedTopography, nCellsOwned, MPI_DOUBLE, gatheredArrayBedTopography, nCellsPerProc, &
nCellsDisplacement, MPI_DOUBLE, 0, domain % dminfo % comm, err_tmp)
err = ior(err, err_tmp)

if (curProc.eq.0) then
call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool)
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)

! First, check consistency in coupling interval set up in MALI and SLM
err = 0
call mpas_pool_get_config(liConfigs, 'config_slm_coupling_interval', config_slm_coupling_interval)
read(config_slm_coupling_interval(1:4),*) slm_coupling_interval
call sl_drive_readnl(itersl, dtime, starttime) !SLM subroutine
if (slm_coupling_interval .NE. dtime) then
call mpas_log_write("The coupling interval in MALI and SLM settings are inconsistent", &
MPAS_LOG_ERR)
err = ior(err,1)
if (curProc.eq.0) then
allocate(globalArrayThickness(nCellsGlobal), gatheredArrayThickness(nCellsGlobal))
allocate(globalArrayBedTopography(nCellsGlobal), gatheredArrayBedTopography(nCellsGlobal))
allocate(meshMask(nCellsGlobal))
ismIceload(:,:) = 0.0
ismBedtopo(:,:) = 0.0
ismMask(:,:) = 0.0
bedtopoSLgrid1D(:) = 0.0
thicknessSLgrid1D(:) = 0.0
maskSLgrid1D(:) = 0.0
else
! Intel requires these be allocated even though they are not meaningful on the non-destination procs
allocate(globalArrayThickness(1), gatheredArrayThickness(1))
allocate(globalArrayBedTopography(1), gatheredArrayBedTopography(1))
allocate(meshMask(1))
endif
! Rearrange data into CellID order
do iCell = 1,nCellsGlobal
globalArrayThickness(indexToCellIDGathered(iCell)) = gatheredArrayThickness(iCell)
globalArrayBedTopography(indexToCellIDGathered(iCell)) = gatheredArrayBedTopography(iCell)
meshMask(indexToCellIDGathered(iCell)) = 1
enddo

! interpolate thickness, bedTopograpy, mesh mask to the Gaussian grid
call interpolate(toColValues, toRowValues, toSvalues, globalArrayThickness, thicknessSLgrid1D)
call interpolate(toColValues, toRowValues, toSvalues, globalArrayBedTopography, bedtopoSLgrid1D)
call interpolate(toColValues, toRowValues, toSvalues, meshMask, maskSLgrid1D)
! Gather only the nCellsOwned from thickness and bedtopo (does not include Halos)
call MPI_GATHERV((thickness*real(li_mask_is_grounded_ice_int(cellMask),RKIND)), &
nCellsOwned, MPI_DOUBLE, gatheredArrayThickness, nCellsPerProc, &
nCellsDisplacement, MPI_DOUBLE, 0, domain % dminfo % comm, err_tmp)
err = ior(err, err_tmp)
call MPI_GATHERV(bedTopography, nCellsOwned, MPI_DOUBLE, gatheredArrayBedTopography, nCellsPerProc, &
nCellsDisplacement, MPI_DOUBLE, 0, domain % dminfo % comm, err_tmp)
err = ior(err, err_tmp)

! reformat the interpolated data
ismIceload = reshape(thicknessSLgrid1D, [nglv,2*nglv])
ismBedtopo = reshape(bedtopoSLgrid1D, [nglv,2*nglv])
ismMask = reshape(maskSLgrid1D, [nglv,2*nglv])
if (curProc.eq.0) then

! initialize coupling time step number. initial time is 0
slmTimeStep = 0
! First, check consistency in coupling interval set up in MALI and SLM
err = 0
call mpas_pool_get_config(liConfigs, 'config_slm_coupling_interval', config_slm_coupling_interval)
read(config_slm_coupling_interval(1:4),*) slm_coupling_interval
call sl_drive_readnl(itersl, dtime, starttime) !SLM subroutine
if (slm_coupling_interval .NE. dtime) then
call mpas_log_write("The coupling interval in MALI and SLM settings are inconsistent", &
MPAS_LOG_ERR)
err = ior(err,1)
endif
! Rearrange data into CellID order
do iCell = 1,nCellsGlobal
globalArrayThickness(indexToCellIDGathered(iCell)) = gatheredArrayThickness(iCell)
globalArrayBedTopography(indexToCellIDGathered(iCell)) = gatheredArrayBedTopography(iCell)
meshMask(indexToCellIDGathered(iCell)) = 1
enddo

! set SLM unit number to the MALI output log file unit
unit_num_slm = domain % logInfo % outputLog % unitNum
! interpolate thickness, bedTopograpy, mesh mask to the Gaussian grid
call interpolate(toColValues, toRowValues, toSvalues, globalArrayThickness, thicknessSLgrid1D)
call interpolate(toColValues, toRowValues, toSvalues, globalArrayBedTopography, bedtopoSLgrid1D)
call interpolate(toColValues, toRowValues, toSvalues, meshMask, maskSLgrid1D)

! series of calling SLM routines
call sl_set_unit_num(unit_num_slm)
call sl_call_readnl
call sl_solver_checkpoint(itersl, dtime)
call sl_timewindow(slmTimeStep)
call sl_solver_init(itersl, starttime, ismIceload, ismBedtopo, ismMask)
call sl_deallocate_array
! reformat the interpolated data
ismIceload = reshape(thicknessSLgrid1D, [nglv,2*nglv])
ismBedtopo = reshape(bedtopoSLgrid1D, [nglv,2*nglv])
ismMask = reshape(maskSLgrid1D, [nglv,2*nglv])

endif
deallocate(globalArrayThickness)
deallocate(gatheredArrayThickness)
deallocate(globalArrayBedTopography)
deallocate(gatheredArrayBedTopography)
deallocate(meshMask)
! initialize coupling time step number. initial time is 0
slmTimeStep = 0

! series of calling SLM routines
call sl_call_readnl
call sl_solver_checkpoint(itersl, dtime)
call sl_timewindow(slmTimeStep)
call sl_solver_init(itersl, starttime, ismIceload, ismBedtopo, ismMask)
call sl_deallocate_array

endif
deallocate(globalArrayThickness)
deallocate(gatheredArrayThickness)
deallocate(globalArrayBedTopography)
deallocate(gatheredArrayBedTopography)
deallocate(meshMask)

endif ! endif restart or not

# else
call mpas_log_write("The sea-level model needs to be included in the compilation with 'SLM=true'", &
Expand Down