Skip to content

Clean up loading of flow parameters and add cfl_min_dt parameter to config #170

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 6 commits into from
Nov 13, 2025
Merged

Clean up loading of flow parameters and add cfl_min_dt parameter to config #170

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
2ac4508
Clean up loading of flow parameters and add cfl_min_dt parameter to c…
btobers Nov 13, 2025
6d24f26
int, float for cfl_min_dt
btobers Nov 13, 2025
fb5d26e
Update cfl params
btobers Nov 13, 2025
4c24a03
Remove nested try/except blocks for simpler debugging
btobers Nov 13, 2025
957988e
Update drive link for sample data in docs
btobers Nov 13, 2025
982a029
Update parser help string
btobers Nov 13, 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
2 changes: 1 addition & 1 deletion docs/install_pygem.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ Following installation, an initialization script should be executed.

The initialization script accomplishes two things:
1. Initializes the PyGEM configuration file *~/PyGEM/config.yaml*. If this file already exists, an overwrite prompt will appear.
2. Downloads and unzips a series of sample data files to *~/PyGEM/*, which can also be manually downloaded [here](https://drive.google.com/file/d/1Wu4ZqpOKxnc4EYhcRHQbwGq95FoOxMfZ/view?usp=drive_link).
2. Downloads and unzips a series of sample data files to *~/PyGEM/*, which can also be manually downloaded [here](https://drive.google.com/drive/folders/1jbQvp0jKXkBZYi47EYhnPCPBcn3We068?usp=sharing).

Run the initialization script by entering the following in the terminal:
```
Expand Down
3,677 changes: 1,818 additions & 1,859 deletions pygem/bin/run/run_calibration.py
View file
Open in desktop

Large diffs are not rendered by default.

33 changes: 19 additions & 14 deletions pygem/bin/run/run_calibration_frontalablation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -249,17 +249,22 @@ def reg_calving_flux(
pygem_prms['root'] + pygem_prms['sim']['oggm_dynamics']['glen_a_regional_relpath']
)
glena_idx = np.where(glena_df.O1Region == glacier_rgi_table.O1Region)[0][0]
glen_a_multiplier = glena_df.loc[glena_idx, 'glens_a_multiplier']
fs = glena_df.loc[glena_idx, 'fs']
# Check which columns exist
# Rounce et al. (2023) regional glen a calibration file has 'glens_a_multiplier' and 'fs'
# output of run_inversion has 'inversion_glen)a' and 'inversion_fs'
if {'glens_a_multiplier', 'fs'}.issubset(glena_df.columns):
glen_a = cfg.PARAMS['glen_a'] * glena_df.loc[glena_idx, 'glens_a_multiplier']
fs = glena_df.loc[glena_idx, 'fs']
elif {'inversion_glen_a', 'inversion_fs'}.issubset(glena_df.columns):
glen_a = glena_df.loc[glena_idx, 'inversion_glen_a']
fs = glena_df.loc[glena_idx, 'inversion_fs']
else:
fs = pygem_prms['sim']['oggm_dynamics']['fs']
glen_a_multiplier = pygem_prms['sim']['oggm_dynamics']['glen_a_multiplier']
glen_a = cfg.PARAMS['glen_a'] * pygem_prms['sim']['oggm_dynamics']['glen_a_multiplier']

# CFL number (may use different values for calving to prevent errors)
if glacier_rgi_table['TermType'] not in [1, 5] or not pygem_prms['setup']['include_frontalablation']:
cfg.PARAMS['cfl_number'] = pygem_prms['sim']['oggm_dynamics']['cfl_number']
else:
cfg.PARAMS['cfl_number'] = pygem_prms['sim']['oggm_dynamics']['cfl_number_calving']
# CFL params
cfg.PARAMS['cfl_number'] = pygem_prms['sim']['oggm_dynamics']['cfl_number']
cfg.PARAMS['cfl_min_dt'] = pygem_prms['sim']['oggm_dynamics']['cfl_min_dt']

# ----- Mass balance model for ice thickness inversion using OGGM -----
mbmod_inv = PyGEMMassBalance(
Expand All @@ -284,12 +289,12 @@ def reg_calving_flux(
if invert_standard:
apparent_mb_from_any_mb(gdir, mb_model=mbmod_inv)
tasks.prepare_for_inversion(gdir)
tasks.mass_conservation_inversion(gdir, glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier, fs=fs)
tasks.mass_conservation_inversion(gdir, glen_a=glen_a, fs=fs)
else:
tasks.find_inversion_calving_from_any_mb(
gdir,
mb_model=mbmod_inv,
glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier,
glen_a=glen_a,
fs=fs,
)

Expand All @@ -310,7 +315,7 @@ def reg_calving_flux(
nfls,
y0=args.ref_startyear,
mb_model=mbmod,
glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier,
glen_a=glen_a,
fs=fs,
is_tidewater=gdir.is_tidewater,
water_level=water_level,
Expand Down Expand Up @@ -374,7 +379,7 @@ def reg_calving_flux(
'OGGM dynamics, calving_k:',
np.round(calving_k, 4),
'glen_a:',
np.round(glen_a_multiplier, 2),
np.round(glen_a, 2),
)
print(
' calving front thickness [m]:',
Expand All @@ -394,7 +399,7 @@ def reg_calving_flux(
nfls,
mb_model=mbmod,
y0=args.ref_startyear,
glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier,
glen_a=glen_a,
fs=fs,
is_tidewater=gdir.is_tidewater,
water_level=water_level,
Expand Down Expand Up @@ -450,7 +455,7 @@ def reg_calving_flux(
'Mass Redistribution curve, calving_k:',
np.round(calving_k, 1),
'glen_a:',
np.round(glen_a_multiplier, 2),
np.round(glen_a, 2),
)
print(
' calving front thickness [m]:',
Expand Down
41 changes: 19 additions & 22 deletions pygem/bin/run/run_simulation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,7 @@
from pygem.output import calc_stats_array
from pygem.plot import graphics
from pygem.shop import debris
from pygem.utils._funcs import str2bool

cfg.PARAMS['hydro_month_nh'] = 1
cfg.PARAMS['hydro_month_sh'] = 1
Expand Down Expand Up @@ -272,9 +273,9 @@ def getparser():
parser.add_argument(
'-use_regional_glen_a',
action='store',
type=bool,
type=str2bool,
default=pygem_prms['sim']['oggm_dynamics']['use_regional_glen_a'],
help='Take the glacier flow parameterization from regionally calibrated priors (boolean: `0` or `1`, `True` or `False`)',
help='If True (False) take glacier flow parameterization from regionally calibrated priors (configuration file).',
)
parser.add_argument(
'-option_bias_adjustment',
Expand Down Expand Up @@ -774,17 +775,10 @@ def run(list_packed_vars):
if debug and gdir.is_tidewater:
print('calving_k:', calving_k)

# Load OGGM glacier dynamics parameters (if necessary)
# Load OGGM glacier dynamics parameters
if args.option_dynamics in ['OGGM', 'MassRedistributionCurves']:
# CFL number (may use different values for calving to prevent errors)
if (
glacier_rgi_table['TermType'] not in [1, 5]
or not pygem_prms['setup']['include_frontalablation']
):
cfg.PARAMS['cfl_number'] = pygem_prms['sim']['oggm_dynamics']['cfl_number']
else:
cfg.PARAMS['cfl_number'] = pygem_prms['sim']['oggm_dynamics']['cfl_number_calving']

cfg.PARAMS['cfl_number'] = pygem_prms['sim']['oggm_dynamics']['cfl_number']
cfg.PARAMS['cfl_min_dt'] = pygem_prms['sim']['oggm_dynamics']['cfl_min_dt']
if debug:
print('cfl number:', cfg.PARAMS['cfl_number'])

Expand All @@ -795,13 +789,18 @@ def run(list_packed_vars):
glena_O1regions = [int(x) for x in glena_df.O1Region.values]
assert glacier_rgi_table.O1Region in glena_O1regions, glacier_str + ' O1 region not in glena_df'
glena_idx = np.where(glena_O1regions == glacier_rgi_table.O1Region)[0][0]
glen_a_multiplier = glena_df.loc[glena_idx, 'glens_a_multiplier']
fs = glena_df.loc[glena_idx, 'fs']
# Check which columns exist
# Rounce et al. (2023) regional glen a calibration file has 'glens_a_multiplier' and 'fs'
# output of run_inversion has 'inversion_glen)a' and 'inversion_fs'
if {'glens_a_multiplier', 'fs'}.issubset(glena_df.columns):
glen_a = cfg.PARAMS['glen_a'] * glena_df.loc[glena_idx, 'glens_a_multiplier']
fs = glena_df.loc[glena_idx, 'fs']
elif {'inversion_glen_a', 'inversion_fs'}.issubset(glena_df.columns):
glen_a = glena_df.loc[glena_idx, 'inversion_glen_a']
fs = glena_df.loc[glena_idx, 'inversion_fs']
else:
args.option_dynamics = None
fs = pygem_prms['sim']['oggm_dynamics']['fs']
glen_a_multiplier = pygem_prms['sim']['oggm_dynamics']['glen_a_multiplier']
glen_a = cfg.PARAMS['glen_a'] * glen_a_multiplier
glen_a = cfg.PARAMS['glen_a'] * pygem_prms['sim']['oggm_dynamics']['glen_a_multiplier']

# spinup
if args.spinup:
Expand All @@ -818,8 +817,6 @@ def run(list_packed_vars):
nfls = get_spinup_flowlines(gdir, y0=args.sim_startyear)
except:
raise
glen_a = gdir.get_diagnostics()['inversion_glen_a']
fs = gdir.get_diagnostics()['inversion_fs']

# Time attributes and values
if pygem_prms['climate']['sim_wateryear'] == 'hydro':
Expand Down Expand Up @@ -914,7 +911,7 @@ def run(list_packed_vars):
tasks.prepare_for_inversion(gdir)
tasks.mass_conservation_inversion(
gdir,
glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier,
glen_a=glen_a,
fs=fs,
)

Expand All @@ -925,7 +922,7 @@ def run(list_packed_vars):
tasks.find_inversion_calving_from_any_mb(
gdir,
mb_model=mbmod_inv,
glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier,
glen_a=glen_a,
fs=fs,
)

Expand Down Expand Up @@ -1065,7 +1062,7 @@ def run(list_packed_vars):
nfls,
mb_model=mbmod,
y0=args.sim_startyear,
glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier,
glen_a=glen_a,
fs=fs,
is_tidewater=gdir.is_tidewater,
# water_level=gdir.get_diagnostics().get('calving_water_level', None)
Expand Down
2 changes: 1 addition & 1 deletion pygem/setup/config.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -352,7 +352,7 @@ def flatten_dict(d, parent_key=''):
'sim.out.export_binned_area_threshold': (int, float),
'sim.oggm_dynamics': dict,
'sim.oggm_dynamics.cfl_number': float,
'sim.oggm_dynamics.cfl_number_calving': float,
'sim.oggm_dynamics.cfl_min_dt': (int, float),
'sim.oggm_dynamics.glen_a_regional_relpath': str,
'sim.oggm_dynamics.use_regional_glen_a': bool,
'sim.oggm_dynamics.fs': (int, float),
Expand Down
4 changes: 2 additions & 2 deletions pygem/setup/config.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -253,8 +253,8 @@ sim:
export_binned_area_threshold: 0 # Area threshold for exporting binned ice thickness
# ===== OGGM DYNAMICS =====
oggm_dynamics:
cfl_number: 0.02 # Time step threshold (seconds)
cfl_number_calving: 0.01 # Time step threshold for marine-terimating glaciers (seconds)
cfl_number: 0.02 # Factor to to us in the CFL criterion to choose the time step
cfl_min_dt: 60 # Time step threshold (seconds)
use_regional_glen_a: true
glen_a_regional_relpath: /Output/calibration/glen_a_region.csv
# glen_a multiplier if not using regionally calibrated glens_a
Expand Down
2 changes: 1 addition & 1 deletion pygem/shop/loso25icebridge.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -486,7 +486,7 @@ def remove_outliers_zscore(self, zscore=3, inplace=False):

def save_elevchange1d(
self,
outdir=f'{pygem_prms["root"]}/{pygem_prms["calib"]["data"]["elev_change_1d"]["elev_change_1d_relpath"]}',
outdir=f'{pygem_prms["root"]}/{pygem_prms["calib"]["data"]["elev_change"]["dh_1d_relpath"]}',
csv=False,
):
"""
Expand Down
Loading