Skip to content

Update glacierdynamics.MassRedistributionCurveModel to work with any timestep #147

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

Merged
btobers merged 10 commits into from
Oct 22, 2025
Merged

Update glacierdynamics.MassRedistributionCurveModel to work with any timestep #147

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
26d5a94
Adapt to work with any timestep
btobers Oct 20, 2025
b4c0df6
Add massredistributioncurves and no dynamics simulation tests
btobers Oct 20, 2025
d6bd105
Update test_05
btobers Oct 20, 2025
5b874bf
Removed unused daysinmonth() fxn
btobers Oct 21, 2025
b84d5de
Use mbmod.dates_table directly
btobers Oct 21, 2025
e5e8fa8
Introduce get_step_inds fxn to return start and stop time indices of …
btobers Oct 21, 2025
7e3b245
Add diagnostic plot legend
btobers Oct 21, 2025
cbba541
Account for nans in comparison
btobers Oct 21, 2025
e38f5aa
Indexing bug fix
btobers Oct 21, 2025
9c422ca
Bug fix
btobers Oct 21, 2025
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
3 changes: 2 additions & 1 deletion .github/workflows/test_suite.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,4 +73,5 @@ jobs:
python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_01_basics.py
python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_02_config.py
python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_03_notebooks.py
python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_04_postproc.py
python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_04_auxiliary.py
python3 -m pytest --cov=pygem -v --durations=20 pygem/tests/test_05_postproc.py
2 changes: 0 additions & 2 deletions pygem/bin/postproc/postproc_binned_subannual_thick.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -333,12 +333,10 @@ def main():
if args.simpath:
# filter out non-file paths
simpath = [p for p in args.simpath if os.path.isfile(p)]

elif args.simdir:
# get list of sims
simpath = sorted(glob.glob(args.simdir + '/*.nc'))
if simpath:
print(simpath)
# number of cores for parallel processing
if args.ncores > 1:
ncores = int(np.min([len(simpath), args.ncores]))
Expand Down
42 changes: 19 additions & 23 deletions pygem/bin/run/run_simulation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -986,7 +986,9 @@ def run(list_packed_vars):

if debug:
fig, ax = plt.subplots(1)
graphics.plot_modeloutput_section(ev_model, ax=ax)
graphics.plot_modeloutput_section(
ev_model, ax=ax, lnlabel=f'Glacier year {args.sim_startyear}'
)

diag = ev_model.run_until_and_store(args.sim_endyear + 1)
ev_model.mb_model.glac_wide_volume_annual[-1] = diag.volume_m3[-1]
Expand Down Expand Up @@ -1056,9 +1058,10 @@ def run(list_packed_vars):
)

if debug:
print('New glacier vol', ev_model.volume_m3)
graphics.plot_modeloutput_section(ev_model)
plt.show()
fig, ax = plt.subplots(1)
graphics.plot_modeloutput_section(
ev_model, ax=ax, lnlabel=f'Glacier year {args.sim_startyear}'
)

_, diag = ev_model.run_until_and_store(args.sim_endyear + 1)
# print('shape of volume:', ev_model.mb_model.glac_wide_volume_annual.shape, diag.volume_m3.shape)
Expand Down Expand Up @@ -1113,25 +1116,20 @@ def run(list_packed_vars):
years = np.arange(args.sim_startyear, args.sim_endyear + 1)
mb_all = []
for year in years:
mb_annual = mbmod.get_annual_mb(
# Calculate annual mass balance
mbmod.get_annual_mb(
nfls[0].surface_h,
fls=nfls,
fl_id=0,
year=year,
debug=True,
debug=debug,
)
mb_mwea = (
mb_annual
* 365
* 24
* 3600
* pygem_prms['constants']['density_ice']
/ pygem_prms['constants']['density_water']
# Record glacierwide annual mass balance
t_start, t_stop = mbmod.get_step_inds(year)
mb_all.append(
mbmod.glac_wide_massbaltotal[t_start : t_stop + 1].sum()
/ mbmod.glacier_area_initial.sum()
)
glac_wide_mb_mwea = (
mb_mwea * mbmod.glacier_area_initial
).sum() / mbmod.glacier_area_initial.sum()
mb_all.append(glac_wide_mb_mwea)
mbmod.glac_wide_area_annual[-1] = mbmod.glac_wide_area_annual[0]
mbmod.glac_wide_volume_annual[-1] = mbmod.glac_wide_volume_annual[0]
diag['area_m2'] = mbmod.glac_wide_area_annual
Expand All @@ -1149,15 +1147,13 @@ def run(list_packed_vars):
np.round(np.median(mb_all), 3),
)

# mb_em_mwea = run_emulator_mb(modelprms)
# print(' emulator mb:', np.round(mb_em_mwea,3))
# mb_em_sims.append(mb_em_mwea)

# Record output for successful runs
if successful_run:
if args.option_dynamics is not None:
if debug:
graphics.plot_modeloutput_section(ev_model, ax=ax, srfls='--')
graphics.plot_modeloutput_section(
ev_model, ax=ax, srfls='--', lnlabel=f'Glacier year {args.sim_endyear + 1}'
)
plt.figure()
diag.volume_m3.plot()
plt.show()
Expand Down Expand Up @@ -1192,7 +1188,7 @@ def run(list_packed_vars):
print(' mb_mbmod [mwea]:', np.round(mb_mwea_mbmod, 2))

if np.abs(mb_mwea_diag - mb_mwea_mbmod) > 1e-6:
ev_model.mb_model.ensure_mass_conservation(diag, dates_table)
ev_model.mb_model.ensure_mass_conservation(diag)

if debug:
print(
Expand Down
36 changes: 17 additions & 19 deletions pygem/glacierdynamics.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,6 @@ def __init__(
inplace=False,
debug=True,
option_areaconstant=False,
spinupyears=0,
constantarea_years=0,
**kwargs,
):
Expand Down Expand Up @@ -87,12 +86,10 @@ def __init__(
)
self.option_areaconstant = option_areaconstant
self.constantarea_years = constantarea_years
self.spinupyears = spinupyears
self.glac_idx_initial = [fl.thick.nonzero()[0] for fl in flowlines]
self.y0 = 0
self.y0 = y0
self.is_tidewater = is_tidewater
self.water_level = water_level

# widths_t0 = flowlines[0].widths_m
# area_v1 = widths_t0 * flowlines[0].dx_meter
# print('area v1:', area_v1.sum())
Expand Down Expand Up @@ -326,8 +323,13 @@ def run_until_and_store(self, y1, run_path=None, diag_path=None, store_monthly_s
def updategeometry(self, year, debug=False):
"""Update geometry for a given year"""

# get year index
year_idx = self.mb_model.get_year_index(year)
# get time step indices - note indexing should be [t_start:t_stop+1] to include final step in year
t_start, t_stop = self.mb_model.get_step_inds(year)
if debug:
print('year:', year)
print('year:', year, f'({year_idx})')
print('time steps:', f'[{t_start}, {t_stop}]')

# Loop over flowlines
for fl_id, fl in enumerate(self.fls):
Expand All @@ -338,8 +340,8 @@ def updategeometry(self, year, debug=False):
width_t0 = self.fls[fl_id].widths_m.copy()

# CONSTANT AREAS
# Mass redistribution ignored for calibration and spinup years (glacier properties constant)
if (self.option_areaconstant) or (year < self.spinupyears) or (year < self.constantarea_years):
# Mass redistribution ignored for calibration years (glacier properties constant)
if (self.option_areaconstant) or (year < self.y0 + self.constantarea_years):
# run mass balance
glac_bin_massbalclim_annual = self.mb_model.get_annual_mb(
heights, fls=self.fls, fl_id=fl_id, year=year, debug=False
Expand Down Expand Up @@ -380,7 +382,7 @@ def updategeometry(self, year, debug=False):
# If frontal ablation more than bin volume, remove entire bin
if fa_m3 > vol_last:
# Record frontal ablation (m3 w.e.) in mass balance model for output
self.mb_model.glac_bin_frontalablation[last_idx, int(12 * (year + 1) - 1)] = (
self.mb_model.glac_bin_frontalablation[last_idx, t_stop + 1] = (
vol_last
* pygem_prms['constants']['density_ice']
/ pygem_prms['constants']['density_water']
Expand All @@ -397,7 +399,7 @@ def updategeometry(self, year, debug=False):
section_t0[last_idx] = section_t0[last_idx] - fa_m3 / fl.dx_meter
self.fls[fl_id].section = section_t0
# Record frontal ablation(m3 w.e.)
self.mb_model.glac_bin_frontalablation[last_idx, int(12 * (year + 1) - 1)] = (
self.mb_model.glac_bin_frontalablation[last_idx, t_stop + 1] = (
fa_m3
* pygem_prms['constants']['density_ice']
/ pygem_prms['constants']['density_water']
Expand Down Expand Up @@ -433,11 +435,8 @@ def updategeometry(self, year, debug=False):
heights, fls=self.fls, fl_id=fl_id, year=year, debug=False
)
sec_in_year = (
self.mb_model.dates_table.loc[12 * year : 12 * (year + 1) - 1, 'daysinmonth'].values.sum()
* 24
* 3600
self.mb_model.dates_table.iloc[t_start : t_stop + 1]['days_in_step'].values.sum() * 24 * 3600
)

# print(' volume change [m3]:', (glac_bin_massbalclim_annual * sec_in_year *
# (width_t0 * fl.dx_meter)).sum())
# print(glac_bin_masssbalclim_annual)
Expand Down Expand Up @@ -469,14 +468,13 @@ def updategeometry(self, year, debug=False):
# Record glacier properties (volume [m3], area [m2], thickness [m], width [km])
# record the next year's properties as well
# 'year + 1' used so the glacier properties are consistent with mass balance computations
year = int(year) # required to ensure proper indexing with run_until_and_store (10/21/2020)
glacier_area = fl.widths_m * fl.dx_meter
glacier_area[fl.thick == 0] = 0
self.mb_model.glac_bin_area_annual[:, year + 1] = glacier_area
self.mb_model.glac_bin_icethickness_annual[:, year + 1] = fl.thick
self.mb_model.glac_bin_width_annual[:, year + 1] = fl.widths_m
self.mb_model.glac_wide_area_annual[year + 1] = glacier_area.sum()
self.mb_model.glac_wide_volume_annual[year + 1] = (fl.section * fl.dx_meter).sum()
self.mb_model.glac_bin_area_annual[:, year_idx + 1] = glacier_area
self.mb_model.glac_bin_icethickness_annual[:, year_idx + 1] = fl.thick
self.mb_model.glac_bin_width_annual[:, year_idx + 1] = fl.widths_m
self.mb_model.glac_wide_area_annual[year_idx + 1] = glacier_area.sum()
self.mb_model.glac_wide_volume_annual[year_idx + 1] = (fl.section * fl.dx_meter).sum()

# %% ----- FRONTAL ABLATION -----
def _get_annual_frontalablation(self, heights, year=None, fls=None, fl_id=None, calving_k=None, debug=False):
Expand Down
Loading