From f0e821fc4b6d8a68532b25a734f5f314f3147df6 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Sun, 26 Mar 2023 09:22:37 -0700 Subject: [PATCH 01/33] reorganization to prepare for branch runs --- .../{ensemble => }/albany_input.yaml | 0 .../__init__.py | 1 - .../ensemble_generator/ensemble_member.py | 19 +++++-------------- .../{ensemble => }/namelist.landice | 0 .../{ensemble => }/streams.landice | 0 5 files changed, 5 insertions(+), 15 deletions(-) rename compass/landice/tests/ensemble_generator/{ensemble => }/albany_input.yaml (100%) rename compass/landice/tests/ensemble_generator/{ensemble => control_ensemble}/__init__.py (98%) rename compass/landice/tests/ensemble_generator/{ensemble => }/namelist.landice (100%) rename compass/landice/tests/ensemble_generator/{ensemble => }/streams.landice (100%) diff --git a/compass/landice/tests/ensemble_generator/ensemble/albany_input.yaml b/compass/landice/tests/ensemble_generator/albany_input.yaml similarity index 100% rename from compass/landice/tests/ensemble_generator/ensemble/albany_input.yaml rename to compass/landice/tests/ensemble_generator/albany_input.yaml diff --git a/compass/landice/tests/ensemble_generator/ensemble/__init__.py b/compass/landice/tests/ensemble_generator/control_ensemble/__init__.py similarity index 98% rename from compass/landice/tests/ensemble_generator/ensemble/__init__.py rename to compass/landice/tests/ensemble_generator/control_ensemble/__init__.py index f8ace1f45a..e062e4934e 100644 --- a/compass/landice/tests/ensemble_generator/ensemble/__init__.py +++ b/compass/landice/tests/ensemble_generator/control_ensemble/__init__.py @@ -175,7 +175,6 @@ def configure(self): for run_num in range(self.start_run, self.end_run + 1): self.add_step(EnsembleMember( test_case=self, run_num=run_num, - test_resources_location='compass.landice.tests.ensemble_generator.ensemble', # noqa basal_fric_exp=param_dict['fric_exp']['vec'][run_num], mu_scale=param_dict['mu_scale']['vec'][run_num], stiff_scale=param_dict['stiff_scale']['vec'][run_num], diff --git a/compass/landice/tests/ensemble_generator/ensemble_member.py b/compass/landice/tests/ensemble_generator/ensemble_member.py index 8c68909f00..c20c503654 100644 --- a/compass/landice/tests/ensemble_generator/ensemble_member.py +++ b/compass/landice/tests/ensemble_generator/ensemble_member.py @@ -28,10 +28,6 @@ class EnsembleMember(Step): ntasks : integer the number of parallel (MPI) tasks the step would ideally use - test_resources_location : str - path to the python package that contains the resources to be - used for the test (namelist, streams, albany input file) - input_file_name : str name of the input file that was read from the config @@ -58,7 +54,6 @@ class EnsembleMember(Step): """ def __init__(self, test_case, run_num, - test_resources_location, basal_fric_exp=None, mu_scale=None, stiff_scale=None, @@ -78,10 +73,6 @@ def __init__(self, test_case, run_num, run_num : integer the run number for this ensemble member - test_resources_location : str - path to the python package that contains the resources to be - used for the test (namelist, streams, albany input file) - basal_fric_exp : float value of basal friction exponent to use @@ -105,7 +96,6 @@ def __init__(self, test_case, run_num, value of deltaT to use in ISMIP6 ice-shelf basal melt param. """ self.run_num = run_num - self.test_resources_location = test_resources_location # store assigned param values for this run self.basal_fric_exp = basal_fric_exp @@ -137,6 +127,8 @@ def setup(self): "'compass setup' again to set this experiment up.") return + module = self.__module__ + # Get config for info needed for setting up simulation config = self.config section = config['ensemble'] @@ -157,14 +149,13 @@ def setup(self): self.min_tasks = self.ntasks # Set up base run configuration - self.add_namelist_file(self.test_resources_location, - 'namelist.landice') + self.add_namelist_file(module, 'namelist.landice') # copy over albany yaml file # cannot use add_input functionality because need to modify the file # in this function, and inputs don't get processed until after this # function - with resources.path(self.test_resources_location, + with resources.path(module, 'albany_input.yaml') as package_path: target = str(package_path) shutil.copy(target, self.work_dir) @@ -260,7 +251,7 @@ def setup(self): # store accumulated namelist and streams options self.add_namelist_options(options=options, out_name='namelist.landice') - self.add_streams_file(self.test_resources_location, 'streams.landice', + self.add_streams_file(module, 'streams.landice', out_name='streams.landice', template_replacements=stream_replacements) diff --git a/compass/landice/tests/ensemble_generator/ensemble/namelist.landice b/compass/landice/tests/ensemble_generator/namelist.landice similarity index 100% rename from compass/landice/tests/ensemble_generator/ensemble/namelist.landice rename to compass/landice/tests/ensemble_generator/namelist.landice diff --git a/compass/landice/tests/ensemble_generator/ensemble/streams.landice b/compass/landice/tests/ensemble_generator/streams.landice similarity index 100% rename from compass/landice/tests/ensemble_generator/ensemble/streams.landice rename to compass/landice/tests/ensemble_generator/streams.landice From 170a8fcf299056a8c00179b607b15f030b2b46a8 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Sun, 26 Mar 2023 21:26:56 -0700 Subject: [PATCH 02/33] Add branch run test case for ensemble --- .../tests/ensemble_generator/__init__.py | 6 +- .../branch_ensemble/__init__.py | 85 +++++++++ .../branch_ensemble/branch_ensemble.cfg | 20 ++ .../branch_ensemble/branch_run.py | 173 ++++++++++++++++++ .../branch_ensemble/namelist.landice | 53 ++++++ .../branch_ensemble/streams.landice | 88 +++++++++ 6 files changed, 424 insertions(+), 1 deletion(-) create mode 100644 compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py create mode 100644 compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg create mode 100644 compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py create mode 100644 compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice create mode 100644 compass/landice/tests/ensemble_generator/branch_ensemble/streams.landice diff --git a/compass/landice/tests/ensemble_generator/__init__.py b/compass/landice/tests/ensemble_generator/__init__.py index 5fd53a238b..e9595be76f 100644 --- a/compass/landice/tests/ensemble_generator/__init__.py +++ b/compass/landice/tests/ensemble_generator/__init__.py @@ -1,4 +1,7 @@ -from compass.landice.tests.ensemble_generator.ensemble import Ensemble +from compass.landice.tests.ensemble_generator.branch_ensemble import ( + BranchEnsemble, +) +from compass.landice.tests.ensemble_generator.control_ensemble import Ensemble from compass.testgroup import TestGroup @@ -16,3 +19,4 @@ def __init__(self, mpas_core): name='ensemble_generator') self.add_test_case(Ensemble(test_group=self)) + self.add_test_case(BranchEnsemble(test_group=self)) diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py b/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py new file mode 100644 index 0000000000..79092fb3c5 --- /dev/null +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py @@ -0,0 +1,85 @@ +import os +import sys + +from compass.landice.tests.ensemble_generator.branch_ensemble.branch_run import ( # noqa + BranchRun, +) +from compass.landice.tests.ensemble_generator.ensemble_manager import ( + EnsembleManager, +) +from compass.testcase import TestCase + + +class BranchEnsemble(TestCase): + """ + A test case for performing an ensemble of + simulations for uncertainty quantification studies. + """ + + def __init__(self, test_group): + """ + Create the test case + + Parameters + ---------- + test_group : compass.landice.tests.ensemble_generator.EnsembleGenerator + The test group that this test case belongs to + + """ + name = 'branch_ensemble' + super().__init__(test_group=test_group, name=name) + + # We don't want to initialize all the individual runs + # So during init, we only add the run manager + self.add_step(EnsembleManager(test_case=self)) + + def configure(self): + """ + Configure a parameter ensemble of MALI simulations. + + Start by identifying the start and end run numbers to set up + from the config. + + Next, read a pre-defined unit parameter vector that can be used + for assigning parameter values to each ensemble member. + + The main work is using the unit parameter vector to set parameter + values for each parameter to be varied, over prescribed ranges. + + Then create the ensemble member as a step in the test case by calling + the EnsembleMember constructor. + + Finally, add this step to the test case's step_to_run. This normally + happens automatically if steps are added to the test case in the test + case constructor, but because we waited to add these steps until this + configure phase, we must explicitly add the steps to steps_to_run. + """ + + config = self.config + section = config['branch_ensemble'] + + control_test_dir = section.get('control_test_dir') + branch_year = section.getint('branch_year') + + # Determine start and end run numbers being requested + self.start_run = config.getint('ensemble', 'start_run') + self.end_run = config.getint('ensemble', 'end_run') + + for run_num in range(self.start_run, self.end_run + 1): + run_name = f'run{run_num:03}' + if os.path.isfile(os.path.join(control_test_dir, run_name, + f'rst.{branch_year}-01-01.nc')): + print(f"Adding {run_name}") + # use this run + self.add_step(BranchRun(test_case=self, run_num=run_num)) + # Note: do not add to steps_to_run, because ensemble_manager + # will handle submitting and running the runs + + # Have compass run only run the run_manager but not any actual runs. + # This is because the individual runs will be submitted as jobs + # by the ensemble manager. + self.steps_to_run = ['ensemble_manager',] + + # no run() method is needed + + # no validate() method is needed diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg new file mode 100644 index 0000000000..559b6374cd --- /dev/null +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg @@ -0,0 +1,20 @@ +# config options for branching an ensemble +[branch_ensemble] + +# start and end numbers for runs to set up and run +# Additional runs can be added and run to an existing ensemble +# without affecting existing runs, but trying to set up a run +# that already exists may result in unexpected behavior. +# Run numbers should be zero-based. +# If using uniform sampling, start_run should be 0 and end_run should be +# equal to max_samples, otherwise unexpected behavior may result. +# These values do not affect viz/analysis, which will include any +# runs it finds. +start_run = 0 +end_run = 3 + +# Path to thermal forcing file for the mesh to be used +TF_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Thwaites_UQ/Thwaites_4to20km_r02_20230126/forcing/ocean_thermal_forcing/obs/Thwaites_4to20km_r02_20230126_obs_TF_1995-2017_8km_x_60m_no_xtime.nc + +# Path to SMB forcing file for the mesh to be used +SMB_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Thwaites_UQ/Thwaites_4to20km_r02_20230126/forcing/atmosphere_forcing/RACMO_climatology_1995-2017/Thwaites_4to20km_r02_20230126_RACMO2.3p2_ANT27_smb_climatology_1995-2017.nc diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py new file mode 100644 index 0000000000..0f3f44f57f --- /dev/null +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py @@ -0,0 +1,173 @@ +import os +import shutil + +import netCDF4 + +import compass.namelist +from compass.io import symlink +from compass.job import write_job_script +from compass.model import run_model +from compass.step import Step + + +class BranchRun(Step): + """ + A step for setting up a single ensemble member + + Attributes + ---------- + run_num : integer + the run number for this ensemble member + + name : str + the name of the run being set up in this step + + ntasks : integer + the number of parallel (MPI) tasks the step would ideally use + + input_file_name : str + name of the input file that was read from the config + + basal_fric_exp : float + value of basal friction exponent to use + + mu_scale : float + value to scale muFriction by + + stiff_scale : float + value to scale stiffnessFactor by + + von_mises_threshold : float + value of von Mises stress threshold to use + + calv_spd_lim : float + value of calving speed limit to use + + gamma0 : float + value of gamma0 to use in ISMIP6 ice-shelf basal melt param. + + deltaT : float + value of deltaT to use in ISMIP6 ice-shelf basal melt param. + """ + + def __init__(self, test_case, run_num, + basal_fric_exp=None, + mu_scale=None, + stiff_scale=None, + von_mises_threshold=None, + calv_spd_lim=None, + gamma0=None, + deltaT=None): + """ + Creates a new run within an ensemble + + Parameters + ---------- + test_case : compass.TestCase + The test case this step belongs to + + run_num : integer + the run number for this ensemble member + """ + self.run_num = run_num + + # define step (run) name + self.name = f'run{run_num:03}' + + super().__init__(test_case=test_case, name=self.name) + + def setup(self): + """ + Set up this run by setting up a baseline run configuration + and then modifying parameters for this ensemble member based on + an externally provided unit parameter vector + """ + + print(f'Setting up run number {self.run_num}') + + config = self.config + section = config['branch_ensemble'] + + control_test_dir = section.get('control_test_dir') + branch_year = section.getint('branch_year') + + ctrl_dir = os.path.join(os.path.join(control_test_dir, self.name)) + + # copy over the following: + # restart file - but change year + rst_file = os.path.join(ctrl_dir, f'rst.{branch_year:04}-01-01.nc') + shutil.copy(rst_file, os.path.join(self.work_dir, + 'rst.2015-01-01.nc')) + f = netCDF4.Dataset(os.path.join(self.work_dir, + 'rst.2015-01-01.nc'), 'r+') + xtime = f.variables['xtime'] + xtime[0, :] = list('2015-01-01_00:00:00'.ljust(64)) + f.close() + + # create restart_timestamp + with open(os.path.join(self.work_dir, 'restart_timestamp'), 'w') as f: + f.write('2015-01-01_00:00:00') + + # yaml file + shutil.copy(os.path.join(ctrl_dir, 'albany_input.yaml'), + self.work_dir) + + # set up namelist + options = {'config_do_restart': '.true.', + 'config_start_time': "'file'", + 'config_stop_time': "'2300-01-01_00:00:00'"} + namelist = compass.namelist.ingest(os.path.join(ctrl_dir, + 'namelist.landice')) + namelist = compass.namelist.replace(namelist, options) + compass.namelist.write(namelist, os.path.join(self.work_dir, + 'namelist.landice')) + + # set up streams + stream_replacements = {} + TF_file_path = section.get('TF_file_path') + stream_replacements['TF_file_path'] = TF_file_path + SMB_file_path = section.get('SMB_file_path') + stream_replacements['SMB_file_path'] = SMB_file_path + strm_src = 'compass.landice.tests.ensemble_generator.branch_ensemble' + self.add_streams_file(strm_src, + 'streams.landice', + out_name='streams.landice', + template_replacements=stream_replacements) + + # copy run_info file + shutil.copy(os.path.join(ctrl_dir, 'run_info.cfg'), self.work_dir) + + # copy graph files + shutil.copy(os.path.join(ctrl_dir, 'graph.info'), self.work_dir) + + # save number of tasks to use + # eventually compass could determine this, but for now we want + # explicit control + self.ntasks = self.config.getint('ensemble', 'ntasks') + self.min_tasks = self.ntasks + + self.add_model_as_input() + + # set job name to run number so it will get set in batch script + self.config.set('job', 'job_name', f'uq_{self.name}') + machine = self.config.get('deploy', 'machine') + write_job_script(self.config, machine, + target_cores=self.ntasks, min_cores=self.min_tasks, + work_dir=self.work_dir) + + # COMPASS does not create symlinks for the load script in step dirs, + # so use the standard approach for creating that symlink in each + # step dir. + if 'LOAD_COMPASS_ENV' in os.environ: + script_filename = os.environ['LOAD_COMPASS_ENV'] + # make a symlink to the script for loading the compass conda env. + symlink(script_filename, os.path.join(self.work_dir, + 'load_compass_env.sh')) + + def run(self): + """ + Run this member of the ensemble. + Eventually we want this function to handle restarts. + """ + + run_model(self) diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice b/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice new file mode 100644 index 0000000000..0d5fec76b9 --- /dev/null +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice @@ -0,0 +1,53 @@ + config_velocity_solver = 'FO' + config_unrealistic_velocity = 0.00317 ! 100 km/yr + config_nonconvergence_error = .false. + config_flowParamA_calculation = 'PB1982' ! required for VM calving + + config_thickness_advection = 'fo' + config_tracer_advection = 'fo' + + config_calving = 'von_Mises_stress' + config_grounded_von_Mises_threshold_stress = 250.0e3 + config_floating_von_Mises_threshold_stress = 250.0e3 + config_calving_speed_limit = 0.000952 ! 30 km/yr + config_damage_calving_threshold = 0.95 + config_calculate_damage = .true. + config_damage_calving_method = 'none' + config_restore_calving_front = .false. + config_remove_icebergs = .true. + config_remove_small_islands = .true. + config_distribute_unablatedVolumeDynCell = .true. + + config_thermal_solver = 'temperature' + config_thermal_calculate_bmb = .true. + config_temperature_init = 'file' + config_thermal_thickness = 0.0 + config_surface_air_temperature_source = 'file' + config_basal_heat_flux_source = 'file' + + config_basal_mass_bal_float = 'ismip6' + config_front_mass_bal_grounded = 'none' + + config_ice_density = 910.0 + config_ocean_density = 1028.0 + config_dynamic_thickness = 10.0 + + config_adaptive_timestep = .true. + config_adaptive_timestep_calvingCFL_fraction = 0.8 + config_adaptive_timestep_include_calving = .true. + config_adaptive_timestep_CFL_fraction = 0.2 + config_adaptive_timestep_force_interval = '0001-00-00_00:00:00' + + config_do_restart = .false. + config_restart_timestamp_name = 'restart_timestamp' + config_start_time = '2000-01-01_00:00:00' + config_stop_time = '2200-01-01_00:00:00' + + config_pio_num_iotasks = 1 + config_pio_stride = 32 + + config_AM_globalStats_enable = .true. + config_AM_globalStats_compute_interval = 'output_interval' + config_AM_globalStats_stream_name = 'globalStats' + config_AM_globalStats_compute_on_startup = .true. + config_AM_globalStats_write_on_startup = .true. diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/streams.landice b/compass/landice/tests/ensemble_generator/branch_ensemble/streams.landice new file mode 100644 index 0000000000..720a9bf80c --- /dev/null +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/streams.landice @@ -0,0 +1,88 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 2325ce405e793d5877ffbe2686d720fe0677c7f9 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 25 Apr 2023 10:15:46 -0700 Subject: [PATCH 03/33] Update test group cfg for Amery ensemble This is just a template, but better to have it match the Amery settings instead of Thwaites, now that we've switched to that. --- .../ensemble_generator/ensemble_generator.cfg | 37 +++++++++---------- 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/ensemble_generator.cfg b/compass/landice/tests/ensemble_generator/ensemble_generator.cfg index e2148133d8..a740675e36 100644 --- a/compass/landice/tests/ensemble_generator/ensemble_generator.cfg +++ b/compass/landice/tests/ensemble_generator/ensemble_generator.cfg @@ -18,7 +18,7 @@ end_run = 3 # for a single parameter sensitivity study. It will sample uniformly across # all dimensions simultaneously, thus sampling only a small fraction of # parameter space -sampling_method = uniform +sampling_method = sobol # maximum number of sample considered. # max_samples needs to be greater or equal to (end_run + 1) @@ -27,7 +27,7 @@ sampling_method = uniform # max_samples should not be changed after the first set of ensemble. # So, when using Sobol sequence, max_samples might be set larger than # (end_run + 1) if you plan to add more samples to the ensemble later. -max_samples = 4 +max_samples = 1024 # basin for comparing model results with observational estimates in # visualization script. @@ -35,7 +35,7 @@ max_samples = 4 # If desired basin does not exist, it can be added to that dataset. # (They need not be mutually exclusive.) # If a basin is not provided, observational comparisons will not be made. -basin = None +basin = ISMIP6BasinBC # fraction of CFL-limited time step to be used by the adaptive timestepper # This value is explicitly included here to force the user to consciously @@ -57,29 +57,28 @@ cfl_fraction = 0.7 # Eventually this could be hard-coded to use files on the input data # server, but initially we want flexibility to experiment with different # inputs and forcings -input_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Thwaites_UQ/Thwaites_4to20km_r02_20230126/relaxation/Thwaites_4to20km_r02_20230126_withStiffness_10yrRelax.nc +input_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Amery_UQ/Amery_4to20km_from_whole_AIS/Amery.nc # the value of the friction exponent used for the calculation of muFriction # in the input file orig_fric_exp = 0.2 # Path to ISMIP6 ice-shelf basal melt parameter input file. -basal_melt_param_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Thwaites_UQ/Thwaites_4to20km_r02_20230126/forcing/basal_melt/parameterizations/Thwaites_4to20km_r02_20230126_basin_and_coeff_gamma0_DeltaT_quadratic_non_local_median.nc +basal_melt_param_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Amery_UQ/Amery_4to20km_from_whole_AIS/forcing/basal_melt/parameterizations/Amery_4to20km_basin_and_coeff_gamma0_DeltaT_quadratic_non_local_median_allBasin2.nc # Path to thermal forcing file for the mesh to be used -TF_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Thwaites_UQ/Thwaites_4to20km_r02_20230126/forcing/ocean_thermal_forcing/obs/Thwaites_4to20km_r02_20230126_obs_TF_1995-2017_8km_x_60m_no_xtime.nc +TF_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Amery_UQ/Amery_4to20km_from_whole_AIS/forcing/ocean_thermal_forcing/obs/Amery_4to20km_obs_TF_1995-2017_8km_x_60m.nc # Path to SMB forcing file for the mesh to be used -SMB_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Thwaites_UQ/Thwaites_4to20km_r02_20230126/forcing/atmosphere_forcing/RACMO_climatology_1995-2017/Thwaites_4to20km_r02_20230126_RACMO2.3p2_ANT27_smb_climatology_1995-2017.nc +SMB_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Amery_UQ/Amery_4to20km_from_whole_AIS/forcing/atmosphere_forcing/RACMO_climatology_1995-2017/Amery_4to20km_RACMO2.3p2_ANT27_smb_climatology_1995-2017_no_xtime_noBareLandAdvance.nc # number of tasks that each ensemble member should be run with # Eventually, compass could determine this, but we want explicit control for now -# ntasks=32 for cori ntasks = 128 # whether basal friction exponent is being varied # [unitless] -use_fric_exp = False +use_fric_exp = True # min value to vary over fric_exp_min = 0.1 # max value to vary over @@ -97,17 +96,17 @@ mu_scale_max = 1.2 # [unitless: 1.0=no scaling] use_stiff_scale = True # min value to vary over -stiff_scale_min = 0.5 +stiff_scale_min = 0.8 # max value to vary over -stiff_scale_max = 1.5 +stiff_scale_max = 1.2 # whether the von Mises threshold stress (sigma_max) is being varied # [units: Pa] -use_von_mises_threshold = False +use_von_mises_threshold = True # min value to vary over -von_mises_threshold_min = 100.0e3 +von_mises_threshold_min = 80.0e3 # max value to vary over -von_mises_threshold_max = 300.0e3 +von_mises_threshold_max = 180.0e3 # whether the calving speed limit is being varied # [units: km/yr] @@ -119,7 +118,7 @@ calv_limit_max = 50.0 # whether ocean melt parameterization coefficient is being varied # [units: m/yr] -use_gamma0 = False +use_gamma0 = True # min value to vary over gamma0_min = 9620.0 # max value to vary over @@ -127,11 +126,11 @@ gamma0_max = 471000.0 # whether target ice-shelf basal melt flux is being varied # [units: Gt/yr] -use_meltflux = False +use_meltflux = True # min value to vary over -meltflux_min = 90.5 +meltflux_min = 12. # max value to vary over -meltflux_max = 114.5 +meltflux_max = 58. # ice-shelf area associated with target melt rates # [units: m^2] -iceshelf_area_obs = 4411.0e6 +iceshelf_area_obs = 60654.e6 From 1604b7fad5b49266128ebb26670f11e11e1bbe63 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 24 Apr 2023 10:31:56 -0700 Subject: [PATCH 04/33] Add grounded volume; adjust GL flux --- .../tests/ensemble_generator/plot_ensemble.py | 24 +++++++++++++++---- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/plot_ensemble.py b/compass/landice/tests/ensemble_generator/plot_ensemble.py index a6d810466f..3962c09e3e 100644 --- a/compass/landice/tests/ensemble_generator/plot_ensemble.py +++ b/compass/landice/tests/ensemble_generator/plot_ensemble.py @@ -73,6 +73,11 @@ 'units': 'km$^2$', 'values': np.zeros((nRuns,)) * np.nan, 'obs': None}, + 'grd vol': { + 'title': f'Grounded vol change at year {targetYear}', + 'units': 'Gt', + 'values': np.zeros((nRuns,)) * np.nan, + 'obs': None}, 'GL flux': { 'title': f'Grounding line flux at year {targetYear}', 'units': 'Gt/yr', @@ -210,8 +215,17 @@ totalArea = f.variables['totalIceArea'][:] / 1.0e6 # in km2 fltArea = f.variables['floatingIceArea'][:] / 1.0e6 # in km2 grdArea = f.variables['groundedIceArea'][:] / 1.0e6 # in km2 - groundingLineFlux = f.variables['groundingLineFlux'][:] / 1.0e12 # Gt - BMB = f.variables['totalFloatingBasalMassBal'][:] / -1.0e12 # in Gt + iceArea = f.variables['totalIceArea'][:] / 1000.0**2 # in km^2 + grdVol = f.variables['groundedIceVolume'][:] / (1.0e12 / rhoi) # in Gt + BMB = f.variables['totalFloatingBasalMassBal'][:] / -1.0e12 # in Gt/yr + SMB = f.variables['totalGroundedSfcMassBal'][:] / -1.0e12 # in Gt/yr + groundingLineFlux = f.variables['groundingLineFlux'][:] \ + / 1.0e12 # Gt/yr + GLMigFlux = f.variables['groundingLineMigrationFlux'][:] \ + / 1.0e12 # Gt/yr + w = 50 + GLMigFlux2 = np.convolve(GLMigFlux, np.ones(w), 'same') / w + GLflux2 = groundingLineFlux + GLMigFlux2 # find target year index indices = np.nonzero(years >= targetYear)[0] @@ -235,7 +249,7 @@ color=col, alpha=alph) axFAts.plot(years, fltArea, linewidth=lw, color=col, alpha=alph) # ignore first entry which is 0 - axGLFts.plot(years[1:], groundingLineFlux[1:], linewidth=lw, + axGLFts.plot(years[1:], GLflux2[1:], linewidth=lw, color=col, alpha=alph) # ignore first entry which is 0 axBMBts.plot(years[1:], BMB[1:], linewidth=lw, @@ -248,10 +262,10 @@ print(f'{run} using year {years[ii]}') qoi_info['SLR']['values'][idx] = SLR[ii] - grdArea = f.variables['groundedIceArea'][:] / 1000.0**2 # in km^2 qoi_info['grd area']['values'][idx] = grdArea[ii] - grdArea[0] - iceArea = f.variables['totalIceArea'][:] / 1000.0**2 # in km^2 + qoi_info['grd vol']['values'][idx] = grdVol[ii] - grdVol[0] + qoi_info['total area']['values'][idx] = iceArea[ii] - iceArea[0] qoi_info['GL flux']['values'][idx] = groundingLineFlux[ii] From e873d4b9d85930cb1bdb7c74157e66d63fc32ae9 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 24 Apr 2023 10:38:08 -0700 Subject: [PATCH 05/33] standardize naming; make all plots optional --- .../tests/ensemble_generator/plot_ensemble.py | 331 +++++++++--------- 1 file changed, 172 insertions(+), 159 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/plot_ensemble.py b/compass/landice/tests/ensemble_generator/plot_ensemble.py index 3962c09e3e..aa9e705f6b 100644 --- a/compass/landice/tests/ensemble_generator/plot_ensemble.py +++ b/compass/landice/tests/ensemble_generator/plot_ensemble.py @@ -16,8 +16,11 @@ # -------------- # general settings # -------------- -targetYear = 100.0 # model year from start at which to calculate statistics -labelRuns = False +target_year = 100.0 # model year from start at which to calculate statistics +label_runs = False +plot_time_series = True +plot_single_param_sensitivies = False +plot_pairwise_param_sensitivities = False plot_maps = False lw = 0.5 # linewidth for ensemble plots @@ -59,32 +62,32 @@ # The values array is 1d array of values from each run qoi_info = { 'SLR': { - 'title': f'SLR at year {targetYear}', + 'title': f'SLR at year {target_year}', 'units': 'mm', 'values': np.zeros((nRuns,)) * np.nan, 'obs': None}, 'total area': { - 'title': f'Total area change at year {targetYear}', + 'title': f'Total area change at year {target_year}', 'units': 'km$^2$', 'values': np.zeros((nRuns,)) * np.nan, 'obs': None}, 'grd area': { - 'title': f'Grounded area change at year {targetYear}', + 'title': f'Grounded area change at year {target_year}', 'units': 'km$^2$', 'values': np.zeros((nRuns,)) * np.nan, 'obs': None}, 'grd vol': { - 'title': f'Grounded vol change at year {targetYear}', + 'title': f'Grounded vol change at year {target_year}', 'units': 'Gt', 'values': np.zeros((nRuns,)) * np.nan, 'obs': None}, 'GL flux': { - 'title': f'Grounding line flux at year {targetYear}', + 'title': f'Grounding line flux at year {target_year}', 'units': 'Gt/yr', 'values': np.zeros((nRuns,)) * np.nan, 'obs': None}, 'melt flux': { - 'title': f'Ice-shelf basal melt flux at year {targetYear}', + 'title': f'Ice-shelf basal melt flux at year {target_year}', 'units': 'Gt/yr', 'values': np.zeros((nRuns,)) * np.nan, 'obs': None}} @@ -131,43 +134,44 @@ # Set up time series plots # -------------- -# Set up axes for time series plots before reading data. -# Time series are plotted as they are read. -figTS = plt.figure(1, figsize=(8, 12), facecolor='w') -nrow = 6 -ncol = 1 -axSLRts = figTS.add_subplot(nrow, ncol, 1) -plt.ylabel('SLR\ncontribution\n(mm)') -plt.grid() - -axTAts = figTS.add_subplot(nrow, ncol, 2, sharex=axSLRts) -plt.ylabel('Total area\nchange (km2)') -plt.grid() - -axGAts = figTS.add_subplot(nrow, ncol, 3, sharex=axSLRts) -plt.ylabel('Grounded area\nchange (km2)') -plt.grid() - -axFAts = figTS.add_subplot(nrow, ncol, 4, sharex=axSLRts) -plt.ylabel('Floating\narea (km2)') -plt.grid() - -axBMBts = figTS.add_subplot(nrow, ncol, 5, sharex=axSLRts) -plt.ylabel('Ice-shelf\nbasal melt\nflux (Gt/yr)') -plt.grid() -axBMBts.fill_between(obs_melt_yrs, - obs_melt - obs_melt_unc, - obs_melt + obs_melt_unc, - color='b', alpha=0.2, label='melt obs') - -axGLFts = figTS.add_subplot(nrow, ncol, 6, sharex=axSLRts) -plt.xlabel('Year') -plt.ylabel('GL flux\n(Gt/yr)') -plt.grid() -axGLFts.fill_between(obs_discharge_yrs, - obs_discharge - obs_discharge_unc, - obs_discharge + obs_discharge_unc, - color='b', alpha=0.2, label='D obs') +if plot_time_series: + # Set up axes for time series plots before reading data. + # Time series are plotted as they are read. + figTS = plt.figure(1, figsize=(8, 12), facecolor='w') + nrow = 6 + ncol = 1 + axSLRts = figTS.add_subplot(nrow, ncol, 1) + plt.ylabel('SLR\ncontribution\n(mm)') + plt.grid() + + axTAts = figTS.add_subplot(nrow, ncol, 2, sharex=axSLRts) + plt.ylabel('Total area\nchange (km2)') + plt.grid() + + axGAts = figTS.add_subplot(nrow, ncol, 3, sharex=axSLRts) + plt.ylabel('Grounded area\nchange (km2)') + plt.grid() + + axFAts = figTS.add_subplot(nrow, ncol, 4, sharex=axSLRts) + plt.ylabel('Floating\narea (km2)') + plt.grid() + + axBMBts = figTS.add_subplot(nrow, ncol, 5, sharex=axSLRts) + plt.ylabel('Ice-shelf\nbasal melt\nflux (Gt/yr)') + plt.grid() + axBMBts.fill_between(obs_melt_yrs, + obs_melt - obs_melt_unc, + obs_melt + obs_melt_unc, + color='b', alpha=0.2, label='melt obs') + + axGLFts = figTS.add_subplot(nrow, ncol, 6, sharex=axSLRts) + plt.xlabel('Year') + plt.ylabel('GL flux\n(Gt/yr)') + plt.grid() + axGLFts.fill_between(obs_discharge_yrs, + obs_discharge - obs_discharge_unc, + obs_discharge + obs_discharge_unc, + color='b', alpha=0.2, label='D obs') # -------------- # maps plotting setup @@ -228,35 +232,36 @@ GLflux2 = groundingLineFlux + GLMigFlux2 # find target year index - indices = np.nonzero(years >= targetYear)[0] - - # color lines depending on if they match obs or not - col = 'k' - alph = 0.2 - GLobs = qoi_info['GL flux']['obs'] - if GLobs is not None and len(indices) > 0: - ii = indices[0] - if groundingLineFlux[ii] > (GLobs[0] - GLobs[1]) and \ - groundingLineFlux[ii] < (GLobs[0] + GLobs[1]): - col = 'r' - alph = 0.7 - - # plot time series - axSLRts.plot(years, SLR, linewidth=lw, color=col, alpha=alph) - axTAts.plot(years, totalArea - totalArea[0], linewidth=lw, - color=col, alpha=alph) - axGAts.plot(years, grdArea - grdArea[0], linewidth=lw, - color=col, alpha=alph) - axFAts.plot(years, fltArea, linewidth=lw, color=col, alpha=alph) - # ignore first entry which is 0 - axGLFts.plot(years[1:], GLflux2[1:], linewidth=lw, - color=col, alpha=alph) - # ignore first entry which is 0 - axBMBts.plot(years[1:], BMB[1:], linewidth=lw, - color=col, alpha=alph) + indices = np.nonzero(years >= target_year)[0] + + if plot_time_series: + # color lines depending on if they match obs or not + col = 'k' + alph = 0.2 + GLobs = qoi_info['GL flux']['obs'] + if GLobs is not None and len(indices) > 0: + ii = indices[0] + if groundingLineFlux[ii] > (GLobs[0] - GLobs[1]) and \ + groundingLineFlux[ii] < (GLobs[0] + GLobs[1]): + col = 'r' + alph = 0.7 + + # plot time series + axSLRts.plot(years, SLR, linewidth=lw, color=col, alpha=alph) + axTAts.plot(years, totalArea - totalArea[0], linewidth=lw, + color=col, alpha=alph) + axGAts.plot(years, grdArea - grdArea[0], linewidth=lw, + color=col, alpha=alph) + axFAts.plot(years, fltArea, linewidth=lw, color=col, alpha=alph) + # ignore first entry which is 0 + axGLFts.plot(years[1:], GLflux2[1:], linewidth=lw, + color=col, alpha=alph) + # ignore first entry which is 0 + axBMBts.plot(years[1:], BMB[1:], linewidth=lw, + color=col, alpha=alph) # Only process runs that have reached target year - indices = np.nonzero(years >= targetYear)[0] + indices = np.nonzero(years >= target_year)[0] if len(indices) > 0: ii = indices[0] print(f'{run} using year {years[ii]}') @@ -279,7 +284,7 @@ decode_timedelta=False, chunks={"Time": 10}) yearsOutput = DS['daysSinceStart'].values[:] / 365.0 - indices = np.nonzero(yearsOutput >= targetYear)[0] + indices = np.nonzero(yearsOutput >= target_year)[0] if len(indices) > 0: ii = indices[0] @@ -329,101 +334,109 @@ axMaps2.hist2d(GLX, GLY, (50, 50), cmap=plt.cm.jet) figMaps.savefig('figure_maps.png') -figTS.tight_layout() -figTS.savefig('figure_time_series.png') +if plot_time_series: + figTS.tight_layout() + figTS.savefig('figure_time_series.png') # -------------- # single parameter plots # -------------- -fig_num = 0 -fig_offset = 100 -for param in param_info: - if param_info[param]['active']: - fig = plt.figure(fig_offset + fig_num, figsize=(13, 8), facecolor='w') - nrow = 2 - ncol = 3 - fig.suptitle(f'{param} sensitivities') - # create subplot for each QOI - n_sub = 1 - for qoi in qoi_info: - ax = fig.add_subplot(nrow, ncol, n_sub) - plt.title(qoi_info[qoi]['title']) - plt.xlabel(f'{param} ({param_info[param]["units"]})') - plt.ylabel(f'{qoi} ({qoi_info[qoi]["units"]})') - pvalues = param_info[param]['values'] - qvalues = qoi_info[qoi]['values'] - obs = qoi_info[qoi]['obs'] - if obs is not None: - plt.fill_between([pvalues.min(), pvalues.max()], - np.array([1., 1.]) * (obs[0] - obs[1]), - np.array([1., 1.]) * (obs[0] + obs[1]), - color='k', alpha=0.2, label='melt obs') - plt.plot(pvalues, qvalues, '.') - if labelRuns: - for i in range(nRuns): - plt.annotate(f'{runs[i][3:]}', (pvalues[i], qvalues[i])) - n_sub += 1 - fig.tight_layout() - fig.savefig(f'figure_sensitivity_{param}.png') - fig_num += 1 +if plot_single_param_sensitivies: + fig_num = 0 + fig_offset = 100 + for param in param_info: + if param_info[param]['active']: + fig = plt.figure(fig_offset + fig_num, figsize=(13, 8), + facecolor='w') + nrow = 2 + ncol = 3 + fig.suptitle(f'{param} sensitivities') + # create subplot for each QOI + n_sub = 1 + for qoi in qoi_info: + ax = fig.add_subplot(nrow, ncol, n_sub) + plt.title(qoi_info[qoi]['title']) + plt.xlabel(f'{param} ({param_info[param]["units"]})') + plt.ylabel(f'{qoi} ({qoi_info[qoi]["units"]})') + pvalues = param_info[param]['values'] + qvalues = qoi_info[qoi]['values'] + obs = qoi_info[qoi]['obs'] + if obs is not None: + plt.fill_between([pvalues.min(), pvalues.max()], + np.array([1., 1.]) * (obs[0] - obs[1]), + np.array([1., 1.]) * (obs[0] + obs[1]), + color='k', alpha=0.2, label='melt obs') + plt.plot(pvalues, qvalues, '.') + if label_runs: + for i in range(nRuns): + plt.annotate(f'{runs[i][3:]}', + (pvalues[i], qvalues[i])) + n_sub += 1 + fig.tight_layout() + fig.savefig(f'figure_sensitivity_{param}.png') + fig_num += 1 # -------------- # pairwise parameter plots # -------------- -fig_num = 0 -fig_offset = 200 -markerSize = 100 -for count1, param1 in enumerate(param_info): - if param_info[param1]['active']: - p1_cnt = count1 - for count2, param2 in enumerate(param_info): - if count2 > count1 and param_info[param2]['active']: - fig = plt.figure(fig_offset + fig_num, figsize=(13, 8), - facecolor='w') - nrow = 2 - ncol = 3 - fig.suptitle(f'{param1} vs. {param2} sensitivities') - # create subplot for each QOI - n_sub = 1 - for qoi in qoi_info: - ax = fig.add_subplot(nrow, ncol, n_sub) - plt.title(f'{qoi_info[qoi]["title"]} ' - f'({qoi_info[qoi]["units"]})') - plt.xlabel(f'{param1} ({param_info[param1]["units"]})') - plt.ylabel(f'{param2} ({param_info[param2]["units"]})') - xdata = param_info[param1]['values'] - ydata = param_info[param2]['values'] - zdata = qoi_info[qoi]['values'] - if np.isfinite(zdata).sum() == 0: - print(f"No valid data for {param1} vs. {param2} " - f"sensitivity plot for {qoi}, skipping") - continue - plt.scatter(xdata, ydata, s=markerSize, c=zdata, - plotnonfinite=False) - badIdx = np.nonzero(np.isnan(zdata))[0] - goodIdx = np.nonzero(np.logical_not(np.isnan(zdata)))[0] - plt.plot(xdata[badIdx], ydata[badIdx], 'kx') - obs = qoi_info[qoi]['obs'] - plt.colorbar() - if obs is not None: - try: - plt.tricontour(xdata[goodIdx], ydata[goodIdx], - zdata[goodIdx], - [obs[0] - obs[1], obs[0] + obs[1]], - colors='k') - except ValueError: - print(f"Skipping obs contour for {param1} vs. " - f"{param2}, because outside model range") - if labelRuns: - for i in range(nRuns): - plt.annotate(f'{runs[i][3:]}', - (xdata[i], ydata[i])) - n_sub += 1 - fig.tight_layout() - fig.savefig( - f'figure_pairwise_sensitivity_{param1}_{param2}.png') - fig_num += 1 +if plot_pairwise_param_sensitivities: + fig_num = 0 + fig_offset = 200 + markerSize = 100 + for count1, param1 in enumerate(param_info): + if param_info[param1]['active']: + p1_cnt = count1 + for count2, param2 in enumerate(param_info): + if count2 > count1 and param_info[param2]['active']: + fig = plt.figure(fig_offset + fig_num, figsize=(13, 8), + facecolor='w') + nrow = 2 + ncol = 3 + fig.suptitle(f'{param1} vs. {param2} sensitivities') + # create subplot for each QOI + n_sub = 1 + for qoi in qoi_info: + ax = fig.add_subplot(nrow, ncol, n_sub) + plt.title(f'{qoi_info[qoi]["title"]} ' + f'({qoi_info[qoi]["units"]})') + plt.xlabel(f'{param1} ({param_info[param1]["units"]})') + plt.ylabel(f'{param2} ({param_info[param2]["units"]})') + xdata = param_info[param1]['values'] + ydata = param_info[param2]['values'] + zdata = qoi_info[qoi]['values'] + if np.isfinite(zdata).sum() == 0: + print(f"No valid data for {param1} vs. {param2} " + f"sensitivity plot for {qoi}, skipping") + continue + plt.scatter(xdata, ydata, s=markerSize, c=zdata, + plotnonfinite=False) + badIdx = np.nonzero(np.isnan(zdata))[0] + goodIdx = np.nonzero(np.logical_not(np.isnan( + zdata)))[0] + plt.plot(xdata[badIdx], ydata[badIdx], 'kx') + obs = qoi_info[qoi]['obs'] + plt.colorbar() + if obs is not None: + try: + plt.tricontour(xdata[goodIdx], ydata[goodIdx], + zdata[goodIdx], + [obs[0] - obs[1], + obs[0] + obs[1]], + colors='k') + except ValueError: + print(f"Skipping obs contour for {param1} " + f"vs. {param2}, because outside model " + "range") + if label_runs: + for i in range(nRuns): + plt.annotate(f'{runs[i][3:]}', + (xdata[i], ydata[i])) + n_sub += 1 + fig.tight_layout() + fig.savefig( + f'figure_pairwise_sensitivity_{param1}_{param2}.png') + fig_num += 1 plt.show() From 77882408b474ab91371bff3c06b1b874001b2073 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 24 Apr 2023 10:55:37 -0700 Subject: [PATCH 06/33] Replace SLR with VAF change everywhere --- .../tests/ensemble_generator/plot_ensemble.py | 47 +++++++++---------- 1 file changed, 22 insertions(+), 25 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/plot_ensemble.py b/compass/landice/tests/ensemble_generator/plot_ensemble.py index aa9e705f6b..9c1aa20027 100644 --- a/compass/landice/tests/ensemble_generator/plot_ensemble.py +++ b/compass/landice/tests/ensemble_generator/plot_ensemble.py @@ -19,7 +19,7 @@ target_year = 100.0 # model year from start at which to calculate statistics label_runs = False plot_time_series = True -plot_single_param_sensitivies = False +plot_single_param_sensitivies = True plot_pairwise_param_sensitivities = False plot_maps = False lw = 0.5 # linewidth for ensemble plots @@ -61,22 +61,22 @@ # Set up nested dictionary for possible quantities of interest. # The values array is 1d array of values from each run qoi_info = { - 'SLR': { - 'title': f'SLR at year {target_year}', - 'units': 'mm', + 'VAF change': { + 'title': f'VAF change at year {target_year}', + 'units': 'Gt', 'values': np.zeros((nRuns,)) * np.nan, 'obs': None}, - 'total area': { + 'total area change': { 'title': f'Total area change at year {target_year}', 'units': 'km$^2$', 'values': np.zeros((nRuns,)) * np.nan, 'obs': None}, - 'grd area': { + 'grd area change': { 'title': f'Grounded area change at year {target_year}', 'units': 'km$^2$', 'values': np.zeros((nRuns,)) * np.nan, 'obs': None}, - 'grd vol': { + 'grd vol change': { 'title': f'Grounded vol change at year {target_year}', 'units': 'Gt', 'values': np.zeros((nRuns,)) * np.nan, @@ -140,23 +140,23 @@ figTS = plt.figure(1, figsize=(8, 12), facecolor='w') nrow = 6 ncol = 1 - axSLRts = figTS.add_subplot(nrow, ncol, 1) - plt.ylabel('SLR\ncontribution\n(mm)') + axVAFts = figTS.add_subplot(nrow, ncol, 1) + plt.ylabel('VAF change\n(mm)') plt.grid() - axTAts = figTS.add_subplot(nrow, ncol, 2, sharex=axSLRts) + axTAts = figTS.add_subplot(nrow, ncol, 2, sharex=axVAFts) plt.ylabel('Total area\nchange (km2)') plt.grid() - axGAts = figTS.add_subplot(nrow, ncol, 3, sharex=axSLRts) + axGAts = figTS.add_subplot(nrow, ncol, 3, sharex=axVAFts) plt.ylabel('Grounded area\nchange (km2)') plt.grid() - axFAts = figTS.add_subplot(nrow, ncol, 4, sharex=axSLRts) + axFAts = figTS.add_subplot(nrow, ncol, 4, sharex=axVAFts) plt.ylabel('Floating\narea (km2)') plt.grid() - axBMBts = figTS.add_subplot(nrow, ncol, 5, sharex=axSLRts) + axBMBts = figTS.add_subplot(nrow, ncol, 5, sharex=axVAFts) plt.ylabel('Ice-shelf\nbasal melt\nflux (Gt/yr)') plt.grid() axBMBts.fill_between(obs_melt_yrs, @@ -164,7 +164,7 @@ obs_melt + obs_melt_unc, color='b', alpha=0.2, label='melt obs') - axGLFts = figTS.add_subplot(nrow, ncol, 6, sharex=axSLRts) + axGLFts = figTS.add_subplot(nrow, ncol, 6, sharex=axVAFts) plt.xlabel('Year') plt.ylabel('GL flux\n(Gt/yr)') plt.grid() @@ -215,7 +215,6 @@ years = f.variables['daysSinceStart'][:] / 365.0 VAF = f.variables['volumeAboveFloatation'][:] - SLR = (VAF[0] - VAF) / 3.62e14 * rhoi / rhosw * 1000. totalArea = f.variables['totalIceArea'][:] / 1.0e6 # in km2 fltArea = f.variables['floatingIceArea'][:] / 1.0e6 # in km2 grdArea = f.variables['groundedIceArea'][:] / 1.0e6 # in km2 @@ -247,7 +246,8 @@ alph = 0.7 # plot time series - axSLRts.plot(years, SLR, linewidth=lw, color=col, alpha=alph) + axVAFts.plot(years, VAF - VAF[0], linewidth=lw, color=col, + alpha=alph) axTAts.plot(years, totalArea - totalArea[0], linewidth=lw, color=col, alpha=alph) axGAts.plot(years, grdArea - grdArea[0], linewidth=lw, @@ -265,16 +265,13 @@ if len(indices) > 0: ii = indices[0] print(f'{run} using year {years[ii]}') - qoi_info['SLR']['values'][idx] = SLR[ii] - - qoi_info['grd area']['values'][idx] = grdArea[ii] - grdArea[0] - - qoi_info['grd vol']['values'][idx] = grdVol[ii] - grdVol[0] - - qoi_info['total area']['values'][idx] = iceArea[ii] - iceArea[0] - + qoi_info['VAF change']['values'][idx] = VAF[ii] - VAF[0] + qoi_info['grd vol change']['values'][idx] = grdVol[ii] - grdVol[0] + qoi_info['grd area change']['values'][idx] = (grdArea[ii] - + grdArea[0]) + qoi_info['total area change']['values'][idx] = (iceArea[ii] - + iceArea[0]) qoi_info['GL flux']['values'][idx] = groundingLineFlux[ii] - qoi_info['melt flux']['values'][idx] = BMB[ii] # plot map From 9ac020e5671d4eb4ea669061c1e2e2d72b14e1d8 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 24 Apr 2023 11:21:54 -0700 Subject: [PATCH 07/33] Updates to time series plots Add more panels, expand across 2 plots --- .../tests/ensemble_generator/plot_ensemble.py | 111 ++++++++++++------ 1 file changed, 74 insertions(+), 37 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/plot_ensemble.py b/compass/landice/tests/ensemble_generator/plot_ensemble.py index 9c1aa20027..fc92181f21 100644 --- a/compass/landice/tests/ensemble_generator/plot_ensemble.py +++ b/compass/landice/tests/ensemble_generator/plot_ensemble.py @@ -19,7 +19,7 @@ target_year = 100.0 # model year from start at which to calculate statistics label_runs = False plot_time_series = True -plot_single_param_sensitivies = True +plot_single_param_sensitivies = False plot_pairwise_param_sensitivities = False plot_maps = False lw = 0.5 # linewidth for ensemble plots @@ -137,41 +137,64 @@ if plot_time_series: # Set up axes for time series plots before reading data. # Time series are plotted as they are read. - figTS = plt.figure(1, figsize=(8, 12), facecolor='w') - nrow = 6 + fig_ts_mb = plt.figure(1, figsize=(8, 12), facecolor='w') + nrow = 5 ncol = 1 - axVAFts = figTS.add_subplot(nrow, ncol, 1) - plt.ylabel('VAF change\n(mm)') + + ax_ts_vaf = fig_ts_mb.add_subplot(nrow, ncol, 1) + plt.ylabel('VAF change\n(Gt)') + plt.grid() + + ax_ts_grvol = fig_ts_mb.add_subplot(nrow, ncol, 2, sharex=ax_ts_vaf) + plt.ylabel('Grounded volume\nchange (Gt)') + plt.grid() + + ax_ts_glf = fig_ts_mb.add_subplot(nrow, ncol, 3, sharex=ax_ts_vaf) + plt.xlabel('Year') + plt.ylabel('GL flux\n(Gt/yr)') + plt.grid() + ax_ts_glf.fill_between(obs_discharge_yrs, + obs_discharge - obs_discharge_unc, + obs_discharge + obs_discharge_unc, + color='b', alpha=0.2, label='D obs') + + ax_ts_glf2 = fig_ts_mb.add_subplot(nrow, ncol, 4, sharex=ax_ts_vaf) + plt.xlabel('Year') + plt.ylabel('GL flux+\nGL mig. flux\n(Gt/yr)') plt.grid() + ax_ts_glf2.fill_between(obs_discharge_yrs, + obs_discharge - obs_discharge_unc, + obs_discharge + obs_discharge_unc, + color='b', alpha=0.2, label='D obs') - axTAts = figTS.add_subplot(nrow, ncol, 2, sharex=axVAFts) + ax_ts_smb = fig_ts_mb.add_subplot(nrow, ncol, 5, sharex=ax_ts_vaf) + plt.ylabel('Grounded SMB\n(Gt/yr)') + plt.grid() + + # second plot to avoid crowding + fig_ts_area = plt.figure(2, figsize=(8, 12), facecolor='w') + nrow = 4 + ncol = 1 + + ax_ts_ta = fig_ts_area.add_subplot(nrow, ncol, 1, sharex=ax_ts_vaf) plt.ylabel('Total area\nchange (km2)') plt.grid() - axGAts = figTS.add_subplot(nrow, ncol, 3, sharex=axVAFts) + ax_ts_ga = fig_ts_area.add_subplot(nrow, ncol, 2, sharex=ax_ts_vaf) plt.ylabel('Grounded area\nchange (km2)') plt.grid() - axFAts = figTS.add_subplot(nrow, ncol, 4, sharex=axVAFts) + ax_ts_fa = fig_ts_area.add_subplot(nrow, ncol, 3, sharex=ax_ts_vaf) plt.ylabel('Floating\narea (km2)') plt.grid() - axBMBts = figTS.add_subplot(nrow, ncol, 5, sharex=axVAFts) + ax_ts_bmb = fig_ts_area.add_subplot(nrow, ncol, 4, sharex=ax_ts_vaf) plt.ylabel('Ice-shelf\nbasal melt\nflux (Gt/yr)') plt.grid() - axBMBts.fill_between(obs_melt_yrs, - obs_melt - obs_melt_unc, - obs_melt + obs_melt_unc, - color='b', alpha=0.2, label='melt obs') - - axGLFts = figTS.add_subplot(nrow, ncol, 6, sharex=axVAFts) - plt.xlabel('Year') - plt.ylabel('GL flux\n(Gt/yr)') - plt.grid() - axGLFts.fill_between(obs_discharge_yrs, - obs_discharge - obs_discharge_unc, - obs_discharge + obs_discharge_unc, - color='b', alpha=0.2, label='D obs') + ax_ts_bmb.fill_between(obs_melt_yrs, + obs_melt - obs_melt_unc, + obs_melt + obs_melt_unc, + color='b', alpha=0.2, label='melt obs') # -------------- # maps plotting setup @@ -214,14 +237,15 @@ f = netCDF4.Dataset(fpath, 'r') years = f.variables['daysSinceStart'][:] / 365.0 - VAF = f.variables['volumeAboveFloatation'][:] + VAF = f.variables['volumeAboveFloatation'][:] / \ + (1.0e12 / rhoi) # in Gt totalArea = f.variables['totalIceArea'][:] / 1.0e6 # in km2 fltArea = f.variables['floatingIceArea'][:] / 1.0e6 # in km2 grdArea = f.variables['groundedIceArea'][:] / 1.0e6 # in km2 iceArea = f.variables['totalIceArea'][:] / 1000.0**2 # in km^2 grdVol = f.variables['groundedIceVolume'][:] / (1.0e12 / rhoi) # in Gt BMB = f.variables['totalFloatingBasalMassBal'][:] / -1.0e12 # in Gt/yr - SMB = f.variables['totalGroundedSfcMassBal'][:] / -1.0e12 # in Gt/yr + grSMB = f.variables['totalGroundedSfcMassBal'][:] / -1.0e12 # in Gt/yr groundingLineFlux = f.variables['groundingLineFlux'][:] \ / 1.0e12 # Gt/yr GLMigFlux = f.variables['groundingLineMigrationFlux'][:] \ @@ -246,19 +270,30 @@ alph = 0.7 # plot time series - axVAFts.plot(years, VAF - VAF[0], linewidth=lw, color=col, - alpha=alph) - axTAts.plot(years, totalArea - totalArea[0], linewidth=lw, - color=col, alpha=alph) - axGAts.plot(years, grdArea - grdArea[0], linewidth=lw, - color=col, alpha=alph) - axFAts.plot(years, fltArea, linewidth=lw, color=col, alpha=alph) + # plot 1 + ax_ts_vaf.plot(years, VAF - VAF[0], linewidth=lw, color=col, + alpha=alph) + ax_ts_grvol.plot(years, grdVol - grdVol[0], linewidth=lw, + color=col, alpha=alph) + # ignore first entry which is 0 + ax_ts_glf.plot(years[1:], groundingLineFlux[1:], linewidth=lw, + color=col, alpha=alph) + # ignore first entry which is 0 + ax_ts_glf2.plot(years[1:], GLflux2[1:], linewidth=lw, + color=col, alpha=alph) # ignore first entry which is 0 - axGLFts.plot(years[1:], GLflux2[1:], linewidth=lw, - color=col, alpha=alph) + ax_ts_smb.plot(years[1:], grSMB[1:], linewidth=lw, + color=col, alpha=alph) + + # plot 2 + ax_ts_ta.plot(years, totalArea - totalArea[0], linewidth=lw, + color=col, alpha=alph) + ax_ts_ga.plot(years, grdArea - grdArea[0], linewidth=lw, + color=col, alpha=alph) + ax_ts_fa.plot(years, fltArea, linewidth=lw, color=col, alpha=alph) # ignore first entry which is 0 - axBMBts.plot(years[1:], BMB[1:], linewidth=lw, - color=col, alpha=alph) + ax_ts_bmb.plot(years[1:], BMB[1:], linewidth=lw, + color=col, alpha=alph) # Only process runs that have reached target year indices = np.nonzero(years >= target_year)[0] @@ -332,8 +367,10 @@ figMaps.savefig('figure_maps.png') if plot_time_series: - figTS.tight_layout() - figTS.savefig('figure_time_series.png') + fig_ts_mb.tight_layout() + fig_ts_mb.savefig('figure_time_series1.png') + fig_ts_area.tight_layout() + fig_ts_area.savefig('figure_time_series2.png') # -------------- # single parameter plots From 6570b243962088fe4ffb0a02fef091b3df2ffec9 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 24 Apr 2023 13:28:59 -0700 Subject: [PATCH 08/33] Add filtering ability to plot script To do so, also add optional filtering criteria to the ais_observations module. --- compass/landice/ais_observations.py | 12 +- .../tests/ensemble_generator/plot_ensemble.py | 127 ++++++++++++------ 2 files changed, 96 insertions(+), 43 deletions(-) diff --git a/compass/landice/ais_observations.py b/compass/landice/ais_observations.py index 7f41fae80f..97ea787e5a 100644 --- a/compass/landice/ais_observations.py +++ b/compass/landice/ais_observations.py @@ -43,7 +43,17 @@ 'input': [73, 10], 'outflow': [77, 4], 'net': [-4, 11], - 'shelf_melt': [35.5, 23.0]}, + 'shelf_melt': [35.5, 23.0], + 'filter criteria': { + 'total area change': { + 'rate': False, + 'values': [-10000.0, 10000.0]}, # km2 + 'grd area change': { + 'rate': True, + 'values': [-15.0, 14.56]}, # km2/yr + 'grd vol change': { + 'rate': True, + 'values': [-18.9, 15.5]}}}, # Gt/yr 'ISMIP6BasinCCp': { 'name': 'Phillipi, Denman', 'input': [81, 13], diff --git a/compass/landice/tests/ensemble_generator/plot_ensemble.py b/compass/landice/tests/ensemble_generator/plot_ensemble.py index fc92181f21..97f7489d27 100644 --- a/compass/landice/tests/ensemble_generator/plot_ensemble.py +++ b/compass/landice/tests/ensemble_generator/plot_ensemble.py @@ -18,9 +18,10 @@ # -------------- target_year = 100.0 # model year from start at which to calculate statistics label_runs = False +filter_runs = True plot_time_series = True -plot_single_param_sensitivies = False -plot_pairwise_param_sensitivities = False +plot_single_param_sensitivies = True +plot_pairwise_param_sensitivities = True plot_maps = False lw = 0.5 # linewidth for ensemble plots @@ -63,34 +64,26 @@ qoi_info = { 'VAF change': { 'title': f'VAF change at year {target_year}', - 'units': 'Gt', - 'values': np.zeros((nRuns,)) * np.nan, - 'obs': None}, + 'units': 'Gt'}, 'total area change': { 'title': f'Total area change at year {target_year}', - 'units': 'km$^2$', - 'values': np.zeros((nRuns,)) * np.nan, - 'obs': None}, + 'units': 'km$^2$'}, 'grd area change': { 'title': f'Grounded area change at year {target_year}', - 'units': 'km$^2$', - 'values': np.zeros((nRuns,)) * np.nan, - 'obs': None}, + 'units': 'km$^2$'}, 'grd vol change': { 'title': f'Grounded vol change at year {target_year}', - 'units': 'Gt', - 'values': np.zeros((nRuns,)) * np.nan, - 'obs': None}, + 'units': 'Gt'}, 'GL flux': { 'title': f'Grounding line flux at year {target_year}', - 'units': 'Gt/yr', - 'values': np.zeros((nRuns,)) * np.nan, - 'obs': None}, + 'units': 'Gt/yr'}, 'melt flux': { 'title': f'Ice-shelf basal melt flux at year {target_year}', - 'units': 'Gt/yr', - 'values': np.zeros((nRuns,)) * np.nan, - 'obs': None}} + 'units': 'Gt/yr'}} +# Initialize some common attributes to be set later +for qoi in qoi_info: + qoi_info[qoi]['values'] = np.zeros((nRuns,)) * np.nan + qoi_info[qoi]['obs'] = None # Get ensemble-wide information basin = None @@ -109,6 +102,29 @@ print(f"Using observations from basin {basin} " f"({ais_basin_info[basin]['name']}).") + +def filter_run(): + valid_run = True + if filter_runs is True and 'filter criteria' in ais_basin_info[basin]: + for criterion in ais_basin_info[basin]['filter criteria']: + # calculate as annual rate + qoi_val = qoi_info[criterion]['values'][idx] + filter_info = ais_basin_info[basin]['filter criteria'][criterion] + min_val = filter_info['values'][0] + max_val = filter_info['values'][1] + if filter_info['rate'] is True: + min_val *= target_year + max_val *= target_year + if qoi_val < min_val or qoi_val > max_val: + valid_run = False + # if run marked invalid, set all qoi to nan + if valid_run is False: + print(f"Marking run {idx} as invalid due to filtering criteria") + for qoi in qoi_info: + qoi_info[qoi]['values'][idx] = np.nan + return valid_run + + # -------------- # Observations information # -------------- @@ -254,20 +270,32 @@ GLMigFlux2 = np.convolve(GLMigFlux, np.ones(w), 'same') / w GLflux2 = groundingLineFlux + GLMigFlux2 - # find target year index + # Only process qois for runs that have reached target year indices = np.nonzero(years >= target_year)[0] + if len(indices) > 0: + ii = indices[0] + print(f'{run} using year {years[ii]}') + qoi_info['VAF change']['values'][idx] = VAF[ii] - VAF[0] + qoi_info['grd vol change']['values'][idx] = grdVol[ii] - grdVol[0] + qoi_info['grd area change']['values'][idx] = (grdArea[ii] - + grdArea[0]) + qoi_info['total area change']['values'][idx] = (iceArea[ii] - + iceArea[0]) + qoi_info['GL flux']['values'][idx] = groundingLineFlux[ii] + qoi_info['melt flux']['values'][idx] = BMB[ii] + + # filter run + valid_run = filter_run() + else: + valid_run = False if plot_time_series: # color lines depending on if they match obs or not col = 'k' alph = 0.2 - GLobs = qoi_info['GL flux']['obs'] - if GLobs is not None and len(indices) > 0: - ii = indices[0] - if groundingLineFlux[ii] > (GLobs[0] - GLobs[1]) and \ - groundingLineFlux[ii] < (GLobs[0] + GLobs[1]): - col = 'r' - alph = 0.7 + if valid_run: + col = 'r' + alph = 0.7 # plot time series # plot 1 @@ -295,20 +323,6 @@ ax_ts_bmb.plot(years[1:], BMB[1:], linewidth=lw, color=col, alpha=alph) - # Only process runs that have reached target year - indices = np.nonzero(years >= target_year)[0] - if len(indices) > 0: - ii = indices[0] - print(f'{run} using year {years[ii]}') - qoi_info['VAF change']['values'][idx] = VAF[ii] - VAF[0] - qoi_info['grd vol change']['values'][idx] = grdVol[ii] - grdVol[0] - qoi_info['grd area change']['values'][idx] = (grdArea[ii] - - grdArea[0]) - qoi_info['total area change']['values'][idx] = (iceArea[ii] - - iceArea[0]) - qoi_info['GL flux']['values'][idx] = groundingLineFlux[ii] - qoi_info['melt flux']['values'][idx] = BMB[ii] - # plot map if plot_maps: DS = xr.open_mfdataset(run + '/output/' + 'output_*.nc', @@ -372,6 +386,33 @@ fig_ts_area.tight_layout() fig_ts_area.savefig('figure_time_series2.png') + # optionally add some filtering indicators on the time series + if filter_runs is True and 'filter criteria' in ais_basin_info[basin]: + filter_info = ais_basin_info[basin]['filter criteria'] + + criterion = 'total area change' + if criterion in filter_info: + min_val = filter_info[criterion]['values'][0] + max_val = filter_info[criterion]['values'][1] + ax_ts_ta.plot([0, target_year], [min_val, min_val], 'b:') + ax_ts_ta.plot([0, target_year], [max_val, max_val], 'b:') + + criterion = 'grd area change' + if criterion in ais_basin_info[basin]['filter criteria']: + min_val = filter_info[criterion]['values'][0] + max_val = filter_info[criterion]['values'][1] + ax_ts_ga.plot([0, target_year], [0, min_val * target_year], 'b:') + ax_ts_ga.plot([0, target_year], [0, max_val * target_year], 'b:') + + criterion = 'grd vol change' + if criterion in ais_basin_info[basin]['filter criteria']: + min_val = filter_info[criterion]['values'][0] + max_val = filter_info[criterion]['values'][1] + ax_ts_grvol.plot([0, target_year], [0, min_val * target_year], + 'b:') + ax_ts_grvol.plot([0, target_year], [0, max_val * target_year], + 'b:') + # -------------- # single parameter plots # -------------- @@ -402,6 +443,8 @@ np.array([1., 1.]) * (obs[0] + obs[1]), color='k', alpha=0.2, label='melt obs') plt.plot(pvalues, qvalues, '.') + badIdx = np.nonzero(np.isnan(qvalues))[0] + plt.plot(pvalues[badIdx], np.zeros(len(badIdx),), 'k.') if label_runs: for i in range(nRuns): plt.annotate(f'{runs[i][3:]}', From e70cf243e42ca6d79303ff206d3e34930e4ecbe1 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 25 Apr 2023 11:36:21 -0700 Subject: [PATCH 09/33] Add info about run statuses in plot script --- .../tests/ensemble_generator/plot_ensemble.py | 20 +++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/plot_ensemble.py b/compass/landice/tests/ensemble_generator/plot_ensemble.py index 97f7489d27..821c25db3d 100644 --- a/compass/landice/tests/ensemble_generator/plot_ensemble.py +++ b/compass/landice/tests/ensemble_generator/plot_ensemble.py @@ -166,7 +166,6 @@ def filter_run(): plt.grid() ax_ts_glf = fig_ts_mb.add_subplot(nrow, ncol, 3, sharex=ax_ts_vaf) - plt.xlabel('Year') plt.ylabel('GL flux\n(Gt/yr)') plt.grid() ax_ts_glf.fill_between(obs_discharge_yrs, @@ -175,7 +174,6 @@ def filter_run(): color='b', alpha=0.2, label='D obs') ax_ts_glf2 = fig_ts_mb.add_subplot(nrow, ncol, 4, sharex=ax_ts_vaf) - plt.xlabel('Year') plt.ylabel('GL flux+\nGL mig. flux\n(Gt/yr)') plt.grid() ax_ts_glf2.fill_between(obs_discharge_yrs, @@ -185,6 +183,7 @@ def filter_run(): ax_ts_smb = fig_ts_mb.add_subplot(nrow, ncol, 5, sharex=ax_ts_vaf) plt.ylabel('Grounded SMB\n(Gt/yr)') + plt.xlabel('Year') plt.grid() # second plot to avoid crowding @@ -206,6 +205,7 @@ def filter_run(): ax_ts_bmb = fig_ts_area.add_subplot(nrow, ncol, 4, sharex=ax_ts_vaf) plt.ylabel('Ice-shelf\nbasal melt\nflux (Gt/yr)') + plt.xlabel('Year') plt.grid() ax_ts_bmb.fill_between(obs_melt_yrs, obs_melt - obs_melt_unc, @@ -233,8 +233,14 @@ def filter_run(): # -------------- # Loop through runs and gather data # -------------- +n_dirs = 0 +n_with_output = 0 +n_at_target_yr = 0 +n_filtered = 0 for idx, run in enumerate(runs): print(f'Analyzing {run}') + n_dirs += 1 + # get param values for this run run_cfg = configparser.ConfigParser() run_cfg.read(os.path.join(run, 'run_info.cfg')) @@ -250,6 +256,7 @@ def filter_run(): fpath = run + "/output/globalStats.nc" if os.path.exists(fpath): + n_with_output += 1 f = netCDF4.Dataset(fpath, 'r') years = f.variables['daysSinceStart'][:] / 365.0 @@ -273,6 +280,7 @@ def filter_run(): # Only process qois for runs that have reached target year indices = np.nonzero(years >= target_year)[0] if len(indices) > 0: + n_at_target_yr += 1 ii = indices[0] print(f'{run} using year {years[ii]}') qoi_info['VAF change']['values'][idx] = VAF[ii] - VAF[0] @@ -286,6 +294,8 @@ def filter_run(): # filter run valid_run = filter_run() + if valid_run: + n_filtered += 1 else: valid_run = False @@ -363,6 +373,12 @@ def filter_run(): f.close() +# Print information about runs +print(f'# runs with directories: {n_dirs}') +print(f'# runs with output files: {n_with_output}') +print(f'# runs reaching target year: {n_at_target_yr}') +print(f'# runs passing filter: {n_filtered}') + # -------------- # save qoi structure # -------------- From 9a032c56835d075ed8856b447256a8e156dadf67 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 25 Apr 2023 12:31:11 -0700 Subject: [PATCH 10/33] Standardize some naming conventions --- compass/landice/tests/ensemble_generator/__init__.py | 6 ++++-- .../tests/ensemble_generator/control_ensemble/__init__.py | 4 ++-- docs/developers_guide/landice/api.rst | 7 +++++-- 3 files changed, 11 insertions(+), 6 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/__init__.py b/compass/landice/tests/ensemble_generator/__init__.py index e9595be76f..7dbb694330 100644 --- a/compass/landice/tests/ensemble_generator/__init__.py +++ b/compass/landice/tests/ensemble_generator/__init__.py @@ -1,7 +1,9 @@ from compass.landice.tests.ensemble_generator.branch_ensemble import ( BranchEnsemble, ) -from compass.landice.tests.ensemble_generator.control_ensemble import Ensemble +from compass.landice.tests.ensemble_generator.control_ensemble import ( + ControlEnsemble, +) from compass.testgroup import TestGroup @@ -18,5 +20,5 @@ def __init__(self, mpas_core): super().__init__(mpas_core=mpas_core, name='ensemble_generator') - self.add_test_case(Ensemble(test_group=self)) + self.add_test_case(ControlEnsemble(test_group=self)) self.add_test_case(BranchEnsemble(test_group=self)) diff --git a/compass/landice/tests/ensemble_generator/control_ensemble/__init__.py b/compass/landice/tests/ensemble_generator/control_ensemble/__init__.py index e062e4934e..6698790665 100644 --- a/compass/landice/tests/ensemble_generator/control_ensemble/__init__.py +++ b/compass/landice/tests/ensemble_generator/control_ensemble/__init__.py @@ -14,7 +14,7 @@ from compass.validate import compare_variables -class Ensemble(TestCase): +class ControlEnsemble(TestCase): """ A test case for performing an ensemble of simulations for uncertainty quantification studies. @@ -30,7 +30,7 @@ def __init__(self, test_group): The test group that this test case belongs to """ - name = 'ensemble' + name = 'control_ensemble' super().__init__(test_group=test_group, name=name) # We don't want to initialize all the individual runs diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst index 02f598c416..8565422f9f 100644 --- a/docs/developers_guide/landice/api.rst +++ b/docs/developers_guide/landice/api.rst @@ -175,8 +175,11 @@ ensemble_generator ensemble_member.EnsembleMember.setup ensemble_member.EnsembleMember.run - ensemble.Ensemble - ensemble.Ensemble.configure + control_ensemble.ControlEnsemble + control_ensemble.ControlEnsemble.configure + + branch_ensemble.BranchEnsemble + branch_ensemble.BranchEnsemble.configure greenland ~~~~~~~~~ From 46476b9447467146ec942fcddb0c46d2d11299fd Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 26 Apr 2023 12:19:08 -0700 Subject: [PATCH 11/33] Fixup to 489f31f38 Correct resource location module name --- .../landice/tests/ensemble_generator/ensemble_member.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/ensemble_member.py b/compass/landice/tests/ensemble_generator/ensemble_member.py index c20c503654..ca08833cff 100644 --- a/compass/landice/tests/ensemble_generator/ensemble_member.py +++ b/compass/landice/tests/ensemble_generator/ensemble_member.py @@ -127,7 +127,7 @@ def setup(self): "'compass setup' again to set this experiment up.") return - module = self.__module__ + resource_module = 'compass.landice.tests.ensemble_generator' # Get config for info needed for setting up simulation config = self.config @@ -149,13 +149,13 @@ def setup(self): self.min_tasks = self.ntasks # Set up base run configuration - self.add_namelist_file(module, 'namelist.landice') + self.add_namelist_file(resource_module, 'namelist.landice') # copy over albany yaml file # cannot use add_input functionality because need to modify the file # in this function, and inputs don't get processed until after this # function - with resources.path(module, + with resources.path(resource_module, 'albany_input.yaml') as package_path: target = str(package_path) shutil.copy(target, self.work_dir) @@ -251,7 +251,7 @@ def setup(self): # store accumulated namelist and streams options self.add_namelist_options(options=options, out_name='namelist.landice') - self.add_streams_file(module, 'streams.landice', + self.add_streams_file(resource_module, 'streams.landice', out_name='streams.landice', template_replacements=stream_replacements) From 8c4287c8104b4d4f769b87931fe798b2347effde Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 28 Apr 2023 10:44:32 -0700 Subject: [PATCH 12/33] Add missing info from branch_ensemble cfg file Also switch to using start_run, end_run from the [branch_ensemble] section --- .../branch_ensemble/__init__.py | 4 +-- .../branch_ensemble/branch_ensemble.cfg | 25 ++++++++++--------- 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py b/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py index 79092fb3c5..1d59e38ca3 100644 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py @@ -62,8 +62,8 @@ def configure(self): branch_year = section.getint('branch_year') # Determine start and end run numbers being requested - self.start_run = config.getint('ensemble', 'start_run') - self.end_run = config.getint('ensemble', 'end_run') + self.start_run = section.getint('start_run') + self.end_run = section.getint('end_run') for run_num in range(self.start_run, self.end_run + 1): run_name = f'run{run_num:03}' diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg index 559b6374cd..0930986e49 100644 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg @@ -2,19 +2,20 @@ [branch_ensemble] # start and end numbers for runs to set up and run -# Additional runs can be added and run to an existing ensemble -# without affecting existing runs, but trying to set up a run -# that already exists may result in unexpected behavior. -# Run numbers should be zero-based. -# If using uniform sampling, start_run should be 0 and end_run should be -# equal to max_samples, otherwise unexpected behavior may result. -# These values do not affect viz/analysis, which will include any -# runs it finds. +# branch runs. +# It is assumed that control runs have already been +# conducted for these runs. start_run = 0 end_run = 3 -# Path to thermal forcing file for the mesh to be used -TF_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Thwaites_UQ/Thwaites_4to20km_r02_20230126/forcing/ocean_thermal_forcing/obs/Thwaites_4to20km_r02_20230126_obs_TF_1995-2017_8km_x_60m_no_xtime.nc +# Path to thermal forcing file for the mesh to be used in the branch run +TF_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Amery_UQ/Amery_4to20km_from_whole_AIS/forcing/ocean_thermal_forcing/UKESM1-0-LL_SSP585/1995-2300/Amery_4to20km_TF_UKESM1-0-LL_SSP585_2300.nc -# Path to SMB forcing file for the mesh to be used -SMB_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Thwaites_UQ/Thwaites_4to20km_r02_20230126/forcing/atmosphere_forcing/RACMO_climatology_1995-2017/Thwaites_4to20km_r02_20230126_RACMO2.3p2_ANT27_smb_climatology_1995-2017.nc +# Path to SMB forcing file for the mesh to be used in the branch run +SMB_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Amery_UQ/Amery_4to20km_from_whole_AIS/forcing/atmosphere_forcing/UKESM1-0-LL_SSP585/1995-2300/Amery_4to20km_SMB_UKESM1-0-LL_SSP585_2300_noBareLandAdvance.nc + +# location of control ensemble to branch from +control_test_dir = /pscratch/sd/h/hoffman2/AMERY_corrected_forcing_6param_ensemble_2023-03-18/landice/ensemble_generator/ensemble + +# year of control simulation from which to branch runs +branch_year = 2050 From 50110109755ad3faf7cad91f50769a5afe410297 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 12 Jun 2023 15:28:52 -0700 Subject: [PATCH 13/33] Add speed error QOI to plot script This commit adds a new quantity of interest for area-averaged error in ice surface speed. Because this QOI requires evaluating spatial information instead of using pre-generated globalStats metrics, this commit also handles opening the output.nc dataset for each run, as well as gathering fields from the input netcdf file. Note that this makes the script run significantly slower. --- .../tests/ensemble_generator/plot_ensemble.py | 47 +++++++++++++++++-- 1 file changed, 43 insertions(+), 4 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/plot_ensemble.py b/compass/landice/tests/ensemble_generator/plot_ensemble.py index 821c25db3d..efadb8172c 100644 --- a/compass/landice/tests/ensemble_generator/plot_ensemble.py +++ b/compass/landice/tests/ensemble_generator/plot_ensemble.py @@ -18,7 +18,7 @@ # -------------- target_year = 100.0 # model year from start at which to calculate statistics label_runs = False -filter_runs = True +filter_runs = False plot_time_series = True plot_single_param_sensitivies = True plot_pairwise_param_sensitivities = True @@ -79,7 +79,12 @@ 'units': 'Gt/yr'}, 'melt flux': { 'title': f'Ice-shelf basal melt flux at year {target_year}', - 'units': 'Gt/yr'}} + 'units': 'Gt/yr'}, + 'speed error': { + 'title': f'Normalized error in modeled ice surface speed at ' + f'year {target_year}', + 'units': 'std. devs.'}} + # Initialize some common attributes to be set later for qoi in qoi_info: qoi_info[qoi]['values'] = np.zeros((nRuns,)) * np.nan @@ -96,6 +101,7 @@ basin = ens_info['basin'] if basin == 'None': basin = None + input_file_path = ens_info['input_file_path'] if basin is None: print("No basin found. Not using observational data.") else: @@ -230,6 +236,17 @@ def filter_run(): GLX = np.array([]) GLY = np.array([]) +# -------------- +# Get needed fields from input file +# -------------- + +DSinput = xr.open_dataset(input_file_path) +observedSurfaceVelocityX = DSinput['observedSurfaceVelocityX'].values +observedSurfaceVelocityY = DSinput['observedSurfaceVelocityY'].values +obsSpdUnc = DSinput['observedSurfaceVelocityUncertainty'].values +obsSpd = (observedSurfaceVelocityX**2 + observedSurfaceVelocityY**2)**0.5 +DSinput.close() + # -------------- # Loop through runs and gather data # -------------- @@ -299,6 +316,28 @@ def filter_run(): else: valid_run = False + # calculate qoi's requiring spatial output + DS = xr.open_mfdataset(run + '/output/' + 'output_*.nc', + combine='nested', concat_dim='Time', + decode_timedelta=False, + chunks={"Time": 10}) + yearsOutput = DS['daysSinceStart'].values[:] / 365.0 + indices = np.nonzero(yearsOutput >= target_year)[0] + if len(indices) > 0: + ii = indices[0] + + surfaceSpeed = DS['surfaceSpeed'].values[ii, :] + thickness = DS['thickness'].values[ii, :] + areaCell = DS['areaCell'].values[0, :] + + spdError = (surfaceSpeed - obsSpd) / obsSpdUnc + # only evaluate where modeled ice exists and observed speed + # uncertainy is a reasonable value (<~1e6 m/yr) + mask = (thickness > 0.0) * (obsSpdUnc < 0.1) + qoi_info['speed error']['values'][idx] = \ + (spdError * mask * areaCell).sum() / (mask * areaCell).sum() + DS.close() + if plot_time_series: # color lines depending on if they match obs or not col = 'k' @@ -440,7 +479,7 @@ def filter_run(): if param_info[param]['active']: fig = plt.figure(fig_offset + fig_num, figsize=(13, 8), facecolor='w') - nrow = 2 + nrow = 3 ncol = 3 fig.suptitle(f'{param} sensitivities') # create subplot for each QOI @@ -485,7 +524,7 @@ def filter_run(): if count2 > count1 and param_info[param2]['active']: fig = plt.figure(fig_offset + fig_num, figsize=(13, 8), facecolor='w') - nrow = 2 + nrow = 3 ncol = 3 fig.suptitle(f'{param1} vs. {param2} sensitivities') # create subplot for each QOI From fd58fc0b0d1a22a204a4559d3c2188688fea69ad Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 12 Jun 2023 15:38:42 -0700 Subject: [PATCH 14/33] Switch previous commit to use log of speed --- compass/landice/tests/ensemble_generator/plot_ensemble.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/plot_ensemble.py b/compass/landice/tests/ensemble_generator/plot_ensemble.py index efadb8172c..17faf3507e 100644 --- a/compass/landice/tests/ensemble_generator/plot_ensemble.py +++ b/compass/landice/tests/ensemble_generator/plot_ensemble.py @@ -16,7 +16,7 @@ # -------------- # general settings # -------------- -target_year = 100.0 # model year from start at which to calculate statistics +target_year = 50.0 # model year from start at which to calculate statistics label_runs = False filter_runs = False plot_time_series = True @@ -330,7 +330,8 @@ def filter_run(): thickness = DS['thickness'].values[ii, :] areaCell = DS['areaCell'].values[0, :] - spdError = (surfaceSpeed - obsSpd) / obsSpdUnc + spdError = ((np.log10(surfaceSpeed) - np.log10(obsSpd)) / + np.log10(obsSpdUnc)) # only evaluate where modeled ice exists and observed speed # uncertainy is a reasonable value (<~1e6 m/yr) mask = (thickness > 0.0) * (obsSpdUnc < 0.1) @@ -499,7 +500,7 @@ def filter_run(): color='k', alpha=0.2, label='melt obs') plt.plot(pvalues, qvalues, '.') badIdx = np.nonzero(np.isnan(qvalues))[0] - plt.plot(pvalues[badIdx], np.zeros(len(badIdx),), 'k.') + plt.plot(pvalues[badIdx], np.zeros(len(badIdx),), 'kx') if label_runs: for i in range(nRuns): plt.annotate(f'{runs[i][3:]}', From 955af8088c3a1b233bc7d1aa57a106bfbb7e8988 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 16 Jun 2023 13:50:43 -0700 Subject: [PATCH 15/33] Update speed error calc; add flt and grd versions --- .../tests/ensemble_generator/plot_ensemble.py | 41 +++++++++++++++---- 1 file changed, 32 insertions(+), 9 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/plot_ensemble.py b/compass/landice/tests/ensemble_generator/plot_ensemble.py index 17faf3507e..fcf13fc75a 100644 --- a/compass/landice/tests/ensemble_generator/plot_ensemble.py +++ b/compass/landice/tests/ensemble_generator/plot_ensemble.py @@ -16,7 +16,7 @@ # -------------- # general settings # -------------- -target_year = 50.0 # model year from start at which to calculate statistics +target_year = 1.0 # model year from start at which to calculate statistics label_runs = False filter_runs = False plot_time_series = True @@ -81,8 +81,16 @@ 'title': f'Ice-shelf basal melt flux at year {target_year}', 'units': 'Gt/yr'}, 'speed error': { - 'title': f'Normalized error in modeled ice surface speed at ' + 'title': f'Speed error at ' f'year {target_year}', + 'units': 'std. devs.'}, + 'flt speed error': { + 'title': f'Speed error over ' + f'floating ice at year {target_year}', + 'units': 'std. devs.'}, + 'grd speed error': { + 'title': f'Speed error over ' + f'grounded ice at year {target_year}', 'units': 'std. devs.'}} # Initialize some common attributes to be set later @@ -241,10 +249,11 @@ def filter_run(): # -------------- DSinput = xr.open_dataset(input_file_path) -observedSurfaceVelocityX = DSinput['observedSurfaceVelocityX'].values -observedSurfaceVelocityY = DSinput['observedSurfaceVelocityY'].values -obsSpdUnc = DSinput['observedSurfaceVelocityUncertainty'].values +observedSurfaceVelocityX = DSinput['observedSurfaceVelocityX'].values[0, :] +observedSurfaceVelocityY = DSinput['observedSurfaceVelocityY'].values[0, :] +obsSpdUnc = DSinput['observedSurfaceVelocityUncertainty'].values[0, :] obsSpd = (observedSurfaceVelocityX**2 + observedSurfaceVelocityY**2)**0.5 +obsSpdUnc += obsSpd * 0.25 DSinput.close() # -------------- @@ -328,15 +337,29 @@ def filter_run(): surfaceSpeed = DS['surfaceSpeed'].values[ii, :] thickness = DS['thickness'].values[ii, :] + cellMask = DS['cellMask'].values[ii, :] areaCell = DS['areaCell'].values[0, :] - spdError = ((np.log10(surfaceSpeed) - np.log10(obsSpd)) / - np.log10(obsSpdUnc)) # only evaluate where modeled ice exists and observed speed # uncertainy is a reasonable value (<~1e6 m/yr) - mask = (thickness > 0.0) * (obsSpdUnc < 0.1) + mask = ((thickness > 0.0) * + (obsSpdUnc < 0.1) * + (obsSpdUnc > 0.0)) + ind = np.nonzero(mask)[0] qoi_info['speed error']['values'][idx] = \ - (spdError * mask * areaCell).sum() / (mask * areaCell).sum() + ((areaCell[ind] * (surfaceSpeed[ind] - obsSpd[ind])**2 / + obsSpdUnc[ind]**2).sum() / areaCell[ind].sum())**0.5 + # float speed error + floatMask = (cellMask & 4) // 4 + ind = np.nonzero(mask * floatMask)[0] + qoi_info['flt speed error']['values'][idx] = \ + ((areaCell[ind] * (surfaceSpeed[ind] - obsSpd[ind])**2 / + obsSpdUnc[ind]**2).sum() / areaCell[ind].sum())**0.5 + # grd speed error + ind = np.nonzero(mask * np.logical_not(floatMask))[0] + qoi_info['grd speed error']['values'][idx] = \ + ((areaCell[ind] * (surfaceSpeed[ind] - obsSpd[ind])**2 / + obsSpdUnc[ind]**2).sum() / areaCell[ind].sum())**0.5 DS.close() if plot_time_series: From e23a4fdc9f9ac95a96ad923e4a5d126830899e22 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 27 Jul 2023 21:51:05 -0700 Subject: [PATCH 16/33] Add scalar calving metrics to output files --- .../ensemble_generator/branch_ensemble/streams.landice | 9 +++++++++ compass/landice/tests/ensemble_generator/streams.landice | 9 +++++++++ 2 files changed, 18 insertions(+) diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/streams.landice b/compass/landice/tests/ensemble_generator/branch_ensemble/streams.landice index 720a9bf80c..d121a5d873 100644 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/streams.landice +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/streams.landice @@ -65,6 +65,12 @@ + + + + + + + + + diff --git a/compass/landice/tests/ensemble_generator/streams.landice b/compass/landice/tests/ensemble_generator/streams.landice index fb3f93c827..2f10d57361 100644 --- a/compass/landice/tests/ensemble_generator/streams.landice +++ b/compass/landice/tests/ensemble_generator/streams.landice @@ -83,6 +83,12 @@ + + + + + + + + + From 1e83993ae3f9ed64c2d074af877aeeb07716d92f Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 27 Jul 2023 21:54:56 -0700 Subject: [PATCH 17/33] Add ability to only branch filtered runs --- .../branch_ensemble/__init__.py | 18 ++++++++++++++++-- .../branch_ensemble/branch_ensemble.cfg | 6 ++++++ 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py b/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py index 1d59e38ca3..cc088a6dda 100644 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py @@ -1,6 +1,9 @@ import os +import pickle import sys +import numpy as np + from compass.landice.tests.ensemble_generator.branch_ensemble.branch_run import ( # noqa BranchRun, ) @@ -65,10 +68,21 @@ def configure(self): self.start_run = section.getint('start_run') self.end_run = section.getint('end_run') + # Determine whether to only set up filtered runs + self.set_up_filtered_only = section.getboolean('set_up_filtered_only') + self.ensemble_pickle_file = section.get('ensemble_pickle_file') + if self.set_up_filtered_only: + with open(self.ensemble_pickle_file, 'rb') as f: + [param_info, qoi_info] = pickle.load(f) + filtered_runs = np.isfinite(qoi_info['VAF change']['values']) + else: + filtered_runs = np.ones((self.end_run + 1,)) + for run_num in range(self.start_run, self.end_run + 1): run_name = f'run{run_num:03}' - if os.path.isfile(os.path.join(control_test_dir, run_name, - f'rst.{branch_year}-01-01.nc')): + if (filtered_runs[run_num] and + os.path.isfile(os.path.join(control_test_dir, run_name, + f'rst.{branch_year}-01-01.nc'))): print(f"Adding {run_name}") # use this run self.add_step(BranchRun(test_case=self, run_num=run_num)) diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg index 0930986e49..e84b0f0679 100644 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg @@ -19,3 +19,9 @@ control_test_dir = /pscratch/sd/h/hoffman2/AMERY_corrected_forcing_6param_ensemb # year of control simulation from which to branch runs branch_year = 2050 + +# whether to only set up branch runs for filtered runs or all runs +set_up_filtered_only = True + +# path to pickle file containing filtering information generated by plot_ensemble.py +ensemble_pickle_file = None From cf52a181c501d71e73e86303cc74b3783138670a Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 31 Aug 2023 18:21:11 -0700 Subject: [PATCH 18/33] Update namelist settings for branch run and spinup --- .../ensemble_generator/branch_ensemble/namelist.landice | 6 ++++-- compass/landice/tests/ensemble_generator/namelist.landice | 3 ++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice b/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice index 0d5fec76b9..265ab387e0 100644 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice @@ -7,9 +7,10 @@ config_tracer_advection = 'fo' config_calving = 'von_Mises_stress' - config_grounded_von_Mises_threshold_stress = 250.0e3 + config_grounded_von_Mises_threshold_stress = 1.0e9 config_floating_von_Mises_threshold_stress = 250.0e3 config_calving_speed_limit = 0.000952 ! 30 km/yr + config_calving_error_threshold = 1.0e9 config_damage_calving_threshold = 0.95 config_calculate_damage = .true. config_damage_calving_method = 'none' @@ -26,7 +27,8 @@ config_basal_heat_flux_source = 'file' config_basal_mass_bal_float = 'ismip6' - config_front_mass_bal_grounded = 'none' + config_front_mass_bal_grounded = 'ismip6' + config_use_3d_thermal_forcing_for_face_melt = .true. config_ice_density = 910.0 config_ocean_density = 1028.0 diff --git a/compass/landice/tests/ensemble_generator/namelist.landice b/compass/landice/tests/ensemble_generator/namelist.landice index 0d5fec76b9..f93af513c1 100644 --- a/compass/landice/tests/ensemble_generator/namelist.landice +++ b/compass/landice/tests/ensemble_generator/namelist.landice @@ -26,7 +26,8 @@ config_basal_heat_flux_source = 'file' config_basal_mass_bal_float = 'ismip6' - config_front_mass_bal_grounded = 'none' + config_front_mass_bal_grounded = 'ismip6' + config_use_3d_thermal_forcing_for_face_melt = .true. config_ice_density = 910.0 config_ocean_density = 1028.0 From 3c2153d0d3c0dc61257980155fe5f43987561d36 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 31 Aug 2023 21:57:34 -0700 Subject: [PATCH 19/33] Add helper script for copy existing runs to new workdir --- .../branch_ensemble/copy_existing_runs.sh | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100755 compass/landice/tests/ensemble_generator/branch_ensemble/copy_existing_runs.sh diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/copy_existing_runs.sh b/compass/landice/tests/ensemble_generator/branch_ensemble/copy_existing_runs.sh new file mode 100755 index 0000000000..2d301d7509 --- /dev/null +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/copy_existing_runs.sh @@ -0,0 +1,18 @@ +#!/bin/bash +# This script can be used to copy over all of the contents of an existing ensemble +# to a fresh generation of the ensemble directory structure. +# This is useful if your compass environment gets corrupted or similar +# situation requiring you to start over with creating the environment. + +SRC_PATH=/pscratch/sd/h/hoffman2/AMERY_corrected_forcing_6param_ensemble_2023-03-18_branch_runs_yr2050_2023-07-28_noGroundedCalving_noCalvingError_calvingMetrics_200runsFiltered/landice/ensemble_generator/branch_ensemble +DEST_PATH=/pscratch/sd/h/hoffman2/AMERY_corrected_forcing_6param_ensemble_2023-03-18_branch_runs_yr2050_2023-08-31/landice/ensemble_generator/branch_ensemble + +for dir in ${SRC_PATH}/run* ; do + echo $dir + run=`basename $dir` + cp ${dir}/log.landice* ${DEST_PATH}/${run}/ + cp ${dir}/restart_timestamp ${DEST_PATH}/${run}/ + cp ${dir}/rst*nc ${DEST_PATH}/${run}/ + cp ${dir}/uq_run*.o* ${DEST_PATH}/${run}/ + cp -R ${dir}/output ${DEST_PATH}/${run}/ +done From 8971c8c73cd07a903653a6f8fa28afe3809e2ea7 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 18 Sep 2023 08:06:30 -0700 Subject: [PATCH 20/33] Update plot script to plot info about run durations --- .../tests/ensemble_generator/plot_ensemble.py | 36 +++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/compass/landice/tests/ensemble_generator/plot_ensemble.py b/compass/landice/tests/ensemble_generator/plot_ensemble.py index fcf13fc75a..6b2482ed05 100644 --- a/compass/landice/tests/ensemble_generator/plot_ensemble.py +++ b/compass/landice/tests/ensemble_generator/plot_ensemble.py @@ -98,6 +98,9 @@ qoi_info[qoi]['values'] = np.zeros((nRuns,)) * np.nan qoi_info[qoi]['obs'] = None +final_time = np.zeros((nRuns,)) * np.nan +run_nums = np.ones((nRuns,), dtype=int) * -1 + # Get ensemble-wide information basin = None ens_cfg = configparser.ConfigParser() @@ -267,6 +270,8 @@ def filter_run(): print(f'Analyzing {run}') n_dirs += 1 + run_nums[idx] = int(run[3:]) + # get param values for this run run_cfg = configparser.ConfigParser() run_cfg.read(os.path.join(run, 'run_info.cfg')) @@ -286,6 +291,8 @@ def filter_run(): f = netCDF4.Dataset(fpath, 'r') years = f.variables['daysSinceStart'][:] / 365.0 + final_time[idx] = years[-1] + VAF = f.variables['volumeAboveFloatation'][:] / \ (1.0e12 / rhoi) # in Gt totalArea = f.variables['totalIceArea'][:] / 1.0e6 # in km2 @@ -492,6 +499,35 @@ def filter_run(): ax_ts_grvol.plot([0, target_year], [0, max_val * target_year], 'b:') +# -------------- +# run duration plots +# -------------- + +fig_duration = plt.figure(99, figsize=(8, 8), facecolor='w') +nrow = 3 +ncol = 1 + +ax_yr_histo = fig_duration.add_subplot(nrow, ncol, 1) +plt.hist(final_time, bins=np.arange(final_time.min(), final_time.max() + 1)) +plt.xlabel('final simulated year') +plt.ylabel('count') +plt.grid() + +ax_yr_histo_cum = fig_duration.add_subplot(nrow, ncol, 2) +plt.hist(final_time, bins=np.arange(final_time.min(), final_time.max() + 1), + cumulative=True) +plt.xlabel('final simulated year') +plt.ylabel('cumulative count') +plt.grid() + +ax_yr_by_run = fig_duration.add_subplot(nrow, ncol, 3) +plt.bar(run_nums, final_time, width=0.9) +plt.xlabel('run number') +plt.ylabel('final simulated year') +plt.grid(axis='y') + +fig_duration.tight_layout() + # -------------- # single parameter plots # -------------- From 5dcf8f944deadf53cd5a164bd2a0f15165e62f4d Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 18 Sep 2023 08:12:45 -0700 Subject: [PATCH 21/33] Allow plot script to use 2 different cfg file names This allows the script to work with both spinup and projection ensembles --- .../tests/ensemble_generator/plot_ensemble.py | 27 ++++++++++++------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/plot_ensemble.py b/compass/landice/tests/ensemble_generator/plot_ensemble.py index 6b2482ed05..5f889eee47 100644 --- a/compass/landice/tests/ensemble_generator/plot_ensemble.py +++ b/compass/landice/tests/ensemble_generator/plot_ensemble.py @@ -4,6 +4,7 @@ import glob import os import pickle +import sys import matplotlib.tri as tri import netCDF4 @@ -104,15 +105,23 @@ # Get ensemble-wide information basin = None ens_cfg = configparser.ConfigParser() -ens_cfg_file = 'ensemble.cfg' -if os.path.isfile(ens_cfg_file): - ens_cfg.read(ens_cfg_file) - ens_info = ens_cfg['ensemble'] - if 'basin' in ens_info: - basin = ens_info['basin'] - if basin == 'None': - basin = None - input_file_path = ens_info['input_file_path'] +# Check for presence of two possible cfg file names +ens_cfg_file1 = 'ensemble.cfg' +ens_cfg_file2 = 'branch_ensemble.cfg' +if os.path.isfile(ens_cfg_file1): + ens_cfg_file = ens_cfg_file1 +elif os.path.isfile(ens_cfg_file2): + ens_cfg_file = ens_cfg_file2 +else: + sys.exit("A usable cfg file for the ensemble was not found. " + "Please correct the configuration or disable this check.") +ens_cfg.read(ens_cfg_file) +ens_info = ens_cfg['ensemble'] +if 'basin' in ens_info: + basin = ens_info['basin'] + if basin == 'None': + basin = None +input_file_path = ens_info['input_file_path'] if basin is None: print("No basin found. Not using observational data.") else: From dfba73c3686c67da85a0b249323327b781998369 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 2 Feb 2024 16:41:50 -0800 Subject: [PATCH 22/33] Skip setting up runs that already exist This is useful for adding additional ensemeble members after an ensemble has already been set up. --- .../ensemble_generator/branch_ensemble/__init__.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py b/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py index cc088a6dda..0503d0fd34 100644 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py @@ -83,11 +83,15 @@ def configure(self): if (filtered_runs[run_num] and os.path.isfile(os.path.join(control_test_dir, run_name, f'rst.{branch_year}-01-01.nc'))): - print(f"Adding {run_name}") - # use this run - self.add_step(BranchRun(test_case=self, run_num=run_num)) - # Note: do not add to steps_to_run, because ensemble_manager - # will handle submitting and running the runs + if os.path.exists(os.path.join(self.work_dir, run_name)): + print(f"WARNING: {run_name} path already exists; " + "skipping. ") + else: + print(f"Adding {run_name}") + # use this run + self.add_step(BranchRun(test_case=self, run_num=run_num)) + # Note: do not add to steps_to_run; ensemble_manager + # will handle submitting and running the runs # Have compass run only run the run_manager but not any actual runs. # This is because the individual runs will be submitted as jobs From 0c7e6cea6505c43506dbeadf71548c435c7d0c63 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 2 Feb 2024 16:44:58 -0800 Subject: [PATCH 23/33] Handle errors and completed runs on submitting a restart In compass run, check for runs that have produced error log files and run that have reached the stop time, and skip submitting them. --- .../ensemble_generator/ensemble_manager.py | 24 ++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/ensemble_manager.py b/compass/landice/tests/ensemble_generator/ensemble_manager.py index a44ab3592b..9e89fbf22e 100644 --- a/compass/landice/tests/ensemble_generator/ensemble_manager.py +++ b/compass/landice/tests/ensemble_generator/ensemble_manager.py @@ -1,8 +1,10 @@ +import glob import os from importlib.resources import path from mpas_tools.logging import check_call +import compass.namelist from compass.io import symlink from compass.step import Step @@ -63,8 +65,24 @@ def run(self): # Get step object from 'steps' dictionary runStep = self.test_case.steps[run] os.chdir(runStep.work_dir) - # TODO: assess if this run is unrun, partially run, or complete, - # and adjust accordingly + + # Skip run if an error is detected + err_files = glob.glob('log.landice.*.err') + if len(err_files) > 0: + print(f'{run} created an error log file in the ' + 'previous job and will not be restarted.') + continue + # Skip run if completed + if os.path.isfile('restart_timestamp'): + with open('restart_timestamp') as f: + curr_time = f.readline().strip() + namelist = compass.namelist.ingest('namelist.landice') + stop_time = (namelist['time_management'] + ['config_stop_time'].strip().strip("'")) + if curr_time == stop_time: + print(f'{run} has reached stop time. Skipping.') + continue + # If we didn't skip this run, submit it check_call(['sbatch', 'job_script.sh'], logger) - logger.info(f'Run {run} submitted.') + logger.info(f'{run} submitted.') logger.info('All runs submitted.') From 20b7c37ed1acf0af5acf20f80a091144abf8fb01 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 2 Feb 2024 17:00:21 -0800 Subject: [PATCH 24/33] Update namelist, streams, submit script for branch runs In the future, the hardcoded namelist updates in branch_run.py should be moved to a separate namelist file. I'm deleting the existing namelist.landice, because it is not actually used. --- .../branch_ensemble/branch_run.py | 14 ++++- .../branch_ensemble/namelist.landice | 55 ------------------- .../branch_ensemble/streams.landice | 4 ++ 3 files changed, 16 insertions(+), 57 deletions(-) delete mode 100644 compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py index 0f3f44f57f..4f148cef6b 100644 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py @@ -115,7 +115,12 @@ def setup(self): # set up namelist options = {'config_do_restart': '.true.', 'config_start_time': "'file'", - 'config_stop_time': "'2300-01-01_00:00:00'"} + 'config_stop_time': "'2300-01-01_00:00:00'", + 'config_grounded_von_Mises_threshold_stress': '1.0e9', + 'config_min_adaptive_timestep': '21600', + 'config_calving_error_threshold': '1.0e9', + 'config_front_mass_bal_grounded': "'ismip6'", + 'config_use_3d_thermal_forcing_for_face_melt': '.true.'} namelist = compass.namelist.ingest(os.path.join(ctrl_dir, 'namelist.landice')) namelist = compass.namelist.replace(namelist, options) @@ -151,9 +156,14 @@ def setup(self): # set job name to run number so it will get set in batch script self.config.set('job', 'job_name', f'uq_{self.name}') machine = self.config.get('deploy', 'machine') + pre_run_cmd = ('LOGDIR=previous_logs_`date +"%Y-%m-%d_%H-%M-%S"`;' + 'mkdir $LOGDIR; cp log* $LOGDIR; date') + post_run_cmd = "date" write_job_script(self.config, machine, target_cores=self.ntasks, min_cores=self.min_tasks, - work_dir=self.work_dir) + work_dir=self.work_dir, + pre_run_commands=pre_run_cmd, + post_run_commands=post_run_cmd) # COMPASS does not create symlinks for the load script in step dirs, # so use the standard approach for creating that symlink in each diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice b/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice deleted file mode 100644 index 265ab387e0..0000000000 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice +++ /dev/null @@ -1,55 +0,0 @@ - config_velocity_solver = 'FO' - config_unrealistic_velocity = 0.00317 ! 100 km/yr - config_nonconvergence_error = .false. - config_flowParamA_calculation = 'PB1982' ! required for VM calving - - config_thickness_advection = 'fo' - config_tracer_advection = 'fo' - - config_calving = 'von_Mises_stress' - config_grounded_von_Mises_threshold_stress = 1.0e9 - config_floating_von_Mises_threshold_stress = 250.0e3 - config_calving_speed_limit = 0.000952 ! 30 km/yr - config_calving_error_threshold = 1.0e9 - config_damage_calving_threshold = 0.95 - config_calculate_damage = .true. - config_damage_calving_method = 'none' - config_restore_calving_front = .false. - config_remove_icebergs = .true. - config_remove_small_islands = .true. - config_distribute_unablatedVolumeDynCell = .true. - - config_thermal_solver = 'temperature' - config_thermal_calculate_bmb = .true. - config_temperature_init = 'file' - config_thermal_thickness = 0.0 - config_surface_air_temperature_source = 'file' - config_basal_heat_flux_source = 'file' - - config_basal_mass_bal_float = 'ismip6' - config_front_mass_bal_grounded = 'ismip6' - config_use_3d_thermal_forcing_for_face_melt = .true. - - config_ice_density = 910.0 - config_ocean_density = 1028.0 - config_dynamic_thickness = 10.0 - - config_adaptive_timestep = .true. - config_adaptive_timestep_calvingCFL_fraction = 0.8 - config_adaptive_timestep_include_calving = .true. - config_adaptive_timestep_CFL_fraction = 0.2 - config_adaptive_timestep_force_interval = '0001-00-00_00:00:00' - - config_do_restart = .false. - config_restart_timestamp_name = 'restart_timestamp' - config_start_time = '2000-01-01_00:00:00' - config_stop_time = '2200-01-01_00:00:00' - - config_pio_num_iotasks = 1 - config_pio_stride = 32 - - config_AM_globalStats_enable = .true. - config_AM_globalStats_compute_interval = 'output_interval' - config_AM_globalStats_stream_name = 'globalStats' - config_AM_globalStats_compute_on_startup = .true. - config_AM_globalStats_write_on_startup = .true. diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/streams.landice b/compass/landice/tests/ensemble_generator/branch_ensemble/streams.landice index d121a5d873..f1cb2af75f 100644 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/streams.landice +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/streams.landice @@ -58,6 +58,10 @@ + + + + From ba9554480999e4ac19917e078579b85157a4db79 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 2 Feb 2024 17:21:15 -0800 Subject: [PATCH 25/33] updates to plotting script MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * add histograms for all QOIs * add rate versions of key QOIs * add labeling to time series plots Note: the new rate versions of QOIs replace the change QOIs in the ais_observations.py module. The ‘rate’ attribute there should be removed and the associated logic in the plotting script should be removed in a future commit. --- compass/landice/ais_observations.py | 16 ++-- .../tests/ensemble_generator/plot_ensemble.py | 85 +++++++++++++++++-- 2 files changed, 84 insertions(+), 17 deletions(-) diff --git a/compass/landice/ais_observations.py b/compass/landice/ais_observations.py index 97ea787e5a..3d76476269 100644 --- a/compass/landice/ais_observations.py +++ b/compass/landice/ais_observations.py @@ -45,15 +45,15 @@ 'net': [-4, 11], 'shelf_melt': [35.5, 23.0], 'filter criteria': { - 'total area change': { + 'total area change rate': { 'rate': False, - 'values': [-10000.0, 10000.0]}, # km2 - 'grd area change': { - 'rate': True, - 'values': [-15.0, 14.56]}, # km2/yr - 'grd vol change': { - 'rate': True, - 'values': [-18.9, 15.5]}}}, # Gt/yr + 'values': [-99.7361194, 75.85552239]}, # km2 + 'grd area change rate': { + 'rate': False, + 'values': [-29.1723, 31.1477]}, # km2/yr + 'grd vol change rate': { + 'rate': False, + 'values': [-24.54, 21.228]}}}, # Gt/yr 'ISMIP6BasinCCp': { 'name': 'Phillipi, Denman', 'input': [81, 13], diff --git a/compass/landice/tests/ensemble_generator/plot_ensemble.py b/compass/landice/tests/ensemble_generator/plot_ensemble.py index 5f889eee47..fb59dee868 100644 --- a/compass/landice/tests/ensemble_generator/plot_ensemble.py +++ b/compass/landice/tests/ensemble_generator/plot_ensemble.py @@ -17,12 +17,16 @@ # -------------- # general settings # -------------- -target_year = 1.0 # model year from start at which to calculate statistics -label_runs = False +target_year = 50.0 # model year from start at which to calculate statistics +# model year from start defining start of time interval over which to +# calculate rates of change +start_year_rate = 0.0 +label_runs = True filter_runs = False plot_time_series = True -plot_single_param_sensitivies = True -plot_pairwise_param_sensitivities = True +plot_single_param_sensitivies = False +plot_pairwise_param_sensitivities = False +plot_qoi_histograms = True plot_maps = False lw = 0.5 # linewidth for ensemble plots @@ -66,15 +70,27 @@ 'VAF change': { 'title': f'VAF change at year {target_year}', 'units': 'Gt'}, + 'VAF change rate': { + 'title': f'VAF change rate between year {start_year_rate} and {target_year}', # noqa + 'units': 'Gt/yr'}, 'total area change': { 'title': f'Total area change at year {target_year}', 'units': 'km$^2$'}, + 'total area change rate': { + 'title': f'Total area change rate between year {start_year_rate} and {target_year}', # noqa + 'units': 'km$^2$/yr'}, 'grd area change': { 'title': f'Grounded area change at year {target_year}', 'units': 'km$^2$'}, + 'grd area change rate': { + 'title': f'Grounded area change rate between year {start_year_rate} and {target_year}', # noqa + 'units': 'km$^2$/yr'}, 'grd vol change': { 'title': f'Grounded vol change at year {target_year}', 'units': 'Gt'}, + 'grd vol change rate': { + 'title': f'Grounded vol change rate between year {start_year_rate} and {target_year}', # noqa + 'units': 'Gt/yr'}, 'GL flux': { 'title': f'Grounding line flux at year {target_year}', 'units': 'Gt/yr'}, @@ -93,6 +109,7 @@ 'title': f'Speed error over ' f'grounded ice at year {target_year}', 'units': 'std. devs.'}} +n_qoi = len(qoi_info) # Initialize some common attributes to be set later for qoi in qoi_info: @@ -324,13 +341,28 @@ def filter_run(): if len(indices) > 0: n_at_target_yr += 1 ii = indices[0] - print(f'{run} using year {years[ii]}') + jj = np.nonzero(years >= start_year_rate)[0][0] + print(f'{run} using year {years[ii]}. ' + f'Using start year for rates of {years[jj]}. ' + f'Max year is {years[-1]}.') + + # find index to start year of rate calculations + rate_dt = years[ii] - years[jj] + qoi_info['VAF change']['values'][idx] = VAF[ii] - VAF[0] + qoi_info['VAF change rate']['values'][idx] = \ + (VAF[ii] - VAF[jj]) / rate_dt qoi_info['grd vol change']['values'][idx] = grdVol[ii] - grdVol[0] + qoi_info['grd vol change rate']['values'][idx] = \ + (grdVol[ii] - grdVol[jj]) / rate_dt qoi_info['grd area change']['values'][idx] = (grdArea[ii] - grdArea[0]) + qoi_info['grd area change rate']['values'][idx] = \ + (grdArea[ii] - grdArea[jj]) / rate_dt qoi_info['total area change']['values'][idx] = (iceArea[ii] - iceArea[0]) + qoi_info['total area change rate']['values'][idx] = \ + (iceArea[ii] - iceArea[jj]) / rate_dt qoi_info['GL flux']['values'][idx] = groundingLineFlux[ii] qoi_info['melt flux']['values'][idx] = BMB[ii] @@ -390,8 +422,14 @@ def filter_run(): # plot 1 ax_ts_vaf.plot(years, VAF - VAF[0], linewidth=lw, color=col, alpha=alph) + if label_runs: + ax_ts_vaf.annotate(run, (years[-1], VAF[-1] - VAF[0]), + fontsize=6) ax_ts_grvol.plot(years, grdVol - grdVol[0], linewidth=lw, color=col, alpha=alph) + if label_runs: + ax_ts_grvol.annotate(run, (years[-1], grdVol[-1] - grdVol[0]), + fontsize=6) # ignore first entry which is 0 ax_ts_glf.plot(years[1:], groundingLineFlux[1:], linewidth=lw, color=col, alpha=alph) @@ -405,8 +443,15 @@ def filter_run(): # plot 2 ax_ts_ta.plot(years, totalArea - totalArea[0], linewidth=lw, color=col, alpha=alph) + if label_runs: + ax_ts_ta.annotate(run, + (years[-1], totalArea[-1] - totalArea[0]), + fontsize=6) ax_ts_ga.plot(years, grdArea - grdArea[0], linewidth=lw, color=col, alpha=alph) + if label_runs: + ax_ts_ga.annotate(run, (years[-1], grdArea[-1] - grdArea[0]), + fontsize=6) ax_ts_fa.plot(years, fltArea, linewidth=lw, color=col, alpha=alph) # ignore first entry which is 0 ax_ts_bmb.plot(years[1:], BMB[1:], linewidth=lw, @@ -548,8 +593,8 @@ def filter_run(): if param_info[param]['active']: fig = plt.figure(fig_offset + fig_num, figsize=(13, 8), facecolor='w') - nrow = 3 - ncol = 3 + ncol = int(np.ceil(n_qoi**0.5)) + nrow = int(np.ceil(n_qoi / ncol)) fig.suptitle(f'{param} sensitivities') # create subplot for each QOI n_sub = 1 @@ -593,8 +638,8 @@ def filter_run(): if count2 > count1 and param_info[param2]['active']: fig = plt.figure(fig_offset + fig_num, figsize=(13, 8), facecolor='w') - nrow = 3 - ncol = 3 + ncol = int(np.ceil(n_qoi**0.5)) + nrow = int(np.ceil(n_qoi / ncol)) fig.suptitle(f'{param1} vs. {param2} sensitivities') # create subplot for each QOI n_sub = 1 @@ -640,4 +685,26 @@ def filter_run(): f'figure_pairwise_sensitivity_{param1}_{param2}.png') fig_num += 1 +# -------------- +# QOI histograms +# -------------- + +if plot_qoi_histograms: + print("Plotting QOI histograms") + fig = plt.figure(300, figsize=(13, 8), facecolor='w') + fig.suptitle('QOI histograms') + ncol = int(np.ceil(n_qoi**0.5)) + nrow = int(np.ceil(n_qoi / ncol)) + n_sub = 0 + for qoi in qoi_info: + n_sub += 1 + ax = fig.add_subplot(nrow, ncol, n_sub) + plt.title(qoi_info[qoi]['title']) + plt.xlabel(f'{qoi} ({qoi_info[qoi]["units"]})') + plt.ylabel('count') + qvalues = qoi_info[qoi]['values'] + ax.hist(qvalues) + fig.tight_layout() + fig.savefig('figure_qoi_histograms.png') + plt.show() From e2eca1b6beab1473cb445c34a1d39e18d222eea3 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 2 Feb 2024 18:36:21 -0800 Subject: [PATCH 26/33] rename control ensemble to spinup ensemble --- .../landice/tests/ensemble_generator/__init__.py | 6 +++--- .../ensemble_generator/branch_ensemble/__init__.py | 4 ++-- .../branch_ensemble/branch_ensemble.cfg | 8 ++++---- .../branch_ensemble/branch_run.py | 14 +++++++------- .../__init__.py | 4 ++-- 5 files changed, 18 insertions(+), 18 deletions(-) rename compass/landice/tests/ensemble_generator/{control_ensemble => spinup_ensemble}/__init__.py (99%) diff --git a/compass/landice/tests/ensemble_generator/__init__.py b/compass/landice/tests/ensemble_generator/__init__.py index 7dbb694330..6dc9af8843 100644 --- a/compass/landice/tests/ensemble_generator/__init__.py +++ b/compass/landice/tests/ensemble_generator/__init__.py @@ -1,8 +1,8 @@ from compass.landice.tests.ensemble_generator.branch_ensemble import ( BranchEnsemble, ) -from compass.landice.tests.ensemble_generator.control_ensemble import ( - ControlEnsemble, +from compass.landice.tests.ensemble_generator.spinup_ensemble import ( + SpinupEnsemble, ) from compass.testgroup import TestGroup @@ -20,5 +20,5 @@ def __init__(self, mpas_core): super().__init__(mpas_core=mpas_core, name='ensemble_generator') - self.add_test_case(ControlEnsemble(test_group=self)) + self.add_test_case(SpinupEnsemble(test_group=self)) self.add_test_case(BranchEnsemble(test_group=self)) diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py b/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py index 0503d0fd34..07c6a4ebec 100644 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/__init__.py @@ -61,7 +61,7 @@ def configure(self): config = self.config section = config['branch_ensemble'] - control_test_dir = section.get('control_test_dir') + spinup_test_dir = section.get('spinup_test_dir') branch_year = section.getint('branch_year') # Determine start and end run numbers being requested @@ -81,7 +81,7 @@ def configure(self): for run_num in range(self.start_run, self.end_run + 1): run_name = f'run{run_num:03}' if (filtered_runs[run_num] and - os.path.isfile(os.path.join(control_test_dir, run_name, + os.path.isfile(os.path.join(spinup_test_dir, run_name, f'rst.{branch_year}-01-01.nc'))): if os.path.exists(os.path.join(self.work_dir, run_name)): print(f"WARNING: {run_name} path already exists; " diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg index e84b0f0679..78953eda17 100644 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_ensemble.cfg @@ -3,7 +3,7 @@ # start and end numbers for runs to set up and run # branch runs. -# It is assumed that control runs have already been +# It is assumed that spinup runs have already been # conducted for these runs. start_run = 0 end_run = 3 @@ -14,10 +14,10 @@ TF_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Amery_UQ/Amery_4to20km_fr # Path to SMB forcing file for the mesh to be used in the branch run SMB_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Amery_UQ/Amery_4to20km_from_whole_AIS/forcing/atmosphere_forcing/UKESM1-0-LL_SSP585/1995-2300/Amery_4to20km_SMB_UKESM1-0-LL_SSP585_2300_noBareLandAdvance.nc -# location of control ensemble to branch from -control_test_dir = /pscratch/sd/h/hoffman2/AMERY_corrected_forcing_6param_ensemble_2023-03-18/landice/ensemble_generator/ensemble +# location of spinup ensemble to branch from +spinup_test_dir = /pscratch/sd/h/hoffman2/AMERY_corrected_forcing_6param_ensemble_2023-03-18/landice/ensemble_generator/ensemble -# year of control simulation from which to branch runs +# year of spinup simulation from which to branch runs branch_year = 2050 # whether to only set up branch runs for filtered runs or all runs diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py index 4f148cef6b..72b2d42efd 100644 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py @@ -88,14 +88,14 @@ def setup(self): config = self.config section = config['branch_ensemble'] - control_test_dir = section.get('control_test_dir') + spinup_test_dir = section.get('spinup_test_dir') branch_year = section.getint('branch_year') - ctrl_dir = os.path.join(os.path.join(control_test_dir, self.name)) + spinup_dir = os.path.join(os.path.join(spinup_test_dir, self.name)) # copy over the following: # restart file - but change year - rst_file = os.path.join(ctrl_dir, f'rst.{branch_year:04}-01-01.nc') + rst_file = os.path.join(spinup_dir, f'rst.{branch_year:04}-01-01.nc') shutil.copy(rst_file, os.path.join(self.work_dir, 'rst.2015-01-01.nc')) f = netCDF4.Dataset(os.path.join(self.work_dir, @@ -109,7 +109,7 @@ def setup(self): f.write('2015-01-01_00:00:00') # yaml file - shutil.copy(os.path.join(ctrl_dir, 'albany_input.yaml'), + shutil.copy(os.path.join(spinup_dir, 'albany_input.yaml'), self.work_dir) # set up namelist @@ -121,7 +121,7 @@ def setup(self): 'config_calving_error_threshold': '1.0e9', 'config_front_mass_bal_grounded': "'ismip6'", 'config_use_3d_thermal_forcing_for_face_melt': '.true.'} - namelist = compass.namelist.ingest(os.path.join(ctrl_dir, + namelist = compass.namelist.ingest(os.path.join(spinup_dir, 'namelist.landice')) namelist = compass.namelist.replace(namelist, options) compass.namelist.write(namelist, os.path.join(self.work_dir, @@ -140,10 +140,10 @@ def setup(self): template_replacements=stream_replacements) # copy run_info file - shutil.copy(os.path.join(ctrl_dir, 'run_info.cfg'), self.work_dir) + shutil.copy(os.path.join(spinup_dir, 'run_info.cfg'), self.work_dir) # copy graph files - shutil.copy(os.path.join(ctrl_dir, 'graph.info'), self.work_dir) + shutil.copy(os.path.join(spinup_dir, 'graph.info'), self.work_dir) # save number of tasks to use # eventually compass could determine this, but for now we want diff --git a/compass/landice/tests/ensemble_generator/control_ensemble/__init__.py b/compass/landice/tests/ensemble_generator/spinup_ensemble/__init__.py similarity index 99% rename from compass/landice/tests/ensemble_generator/control_ensemble/__init__.py rename to compass/landice/tests/ensemble_generator/spinup_ensemble/__init__.py index 6698790665..1b6aae8a80 100644 --- a/compass/landice/tests/ensemble_generator/control_ensemble/__init__.py +++ b/compass/landice/tests/ensemble_generator/spinup_ensemble/__init__.py @@ -14,7 +14,7 @@ from compass.validate import compare_variables -class ControlEnsemble(TestCase): +class SpinupEnsemble(TestCase): """ A test case for performing an ensemble of simulations for uncertainty quantification studies. @@ -30,7 +30,7 @@ def __init__(self, test_group): The test group that this test case belongs to """ - name = 'control_ensemble' + name = 'spinup_ensemble' super().__init__(test_group=test_group, name=name) # We don't want to initialize all the individual runs From 495015260b3c1e527d227e96857a6e3074b0ff37 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 2 Feb 2024 19:02:08 -0800 Subject: [PATCH 27/33] remove 'rate' attribute that is no longer used --- compass/landice/ais_observations.py | 3 --- compass/landice/tests/ensemble_generator/plot_ensemble.py | 3 --- 2 files changed, 6 deletions(-) diff --git a/compass/landice/ais_observations.py b/compass/landice/ais_observations.py index 3d76476269..335285ad2c 100644 --- a/compass/landice/ais_observations.py +++ b/compass/landice/ais_observations.py @@ -46,13 +46,10 @@ 'shelf_melt': [35.5, 23.0], 'filter criteria': { 'total area change rate': { - 'rate': False, 'values': [-99.7361194, 75.85552239]}, # km2 'grd area change rate': { - 'rate': False, 'values': [-29.1723, 31.1477]}, # km2/yr 'grd vol change rate': { - 'rate': False, 'values': [-24.54, 21.228]}}}, # Gt/yr 'ISMIP6BasinCCp': { 'name': 'Phillipi, Denman', diff --git a/compass/landice/tests/ensemble_generator/plot_ensemble.py b/compass/landice/tests/ensemble_generator/plot_ensemble.py index fb59dee868..f6300060d3 100644 --- a/compass/landice/tests/ensemble_generator/plot_ensemble.py +++ b/compass/landice/tests/ensemble_generator/plot_ensemble.py @@ -155,9 +155,6 @@ def filter_run(): filter_info = ais_basin_info[basin]['filter criteria'][criterion] min_val = filter_info['values'][0] max_val = filter_info['values'][1] - if filter_info['rate'] is True: - min_val *= target_year - max_val *= target_year if qoi_val < min_val or qoi_val > max_val: valid_run = False # if run marked invalid, set all qoi to nan From 6f0fef25aa03b7e4f1fc6d2b72b71ac1f95ae6d8 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 2 Feb 2024 19:43:43 -0800 Subject: [PATCH 28/33] parse nl replacements from branch_ensemble module this replaces them being hardcoded in the .py file --- .../branch_ensemble/branch_run.py | 15 +++++++-------- .../branch_ensemble/namelist.landice | 8 ++++++++ 2 files changed, 15 insertions(+), 8 deletions(-) create mode 100644 compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py index 72b2d42efd..cf0a30e76c 100644 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py @@ -113,16 +113,15 @@ def setup(self): self.work_dir) # set up namelist - options = {'config_do_restart': '.true.', - 'config_start_time': "'file'", - 'config_stop_time': "'2300-01-01_00:00:00'", - 'config_grounded_von_Mises_threshold_stress': '1.0e9', - 'config_min_adaptive_timestep': '21600', - 'config_calving_error_threshold': '1.0e9', - 'config_front_mass_bal_grounded': "'ismip6'", - 'config_use_3d_thermal_forcing_for_face_melt': '.true.'} + # start with the namelist from the spinup + # Note: this differs from most compass tests, which would start with + # the default namelist from the mpas build dir namelist = compass.namelist.ingest(os.path.join(spinup_dir, 'namelist.landice')) + # use the namelist in this module to update the spinup namelist + options = compass.namelist.parse_replacements( + 'compass.landice.tests.ensemble_generator.branch_ensemble', + 'namelist.landice') namelist = compass.namelist.replace(namelist, options) compass.namelist.write(namelist, os.path.join(self.work_dir, 'namelist.landice')) diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice b/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice new file mode 100644 index 0000000000..d09cb4290f --- /dev/null +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice @@ -0,0 +1,8 @@ +config_do_restart = '.true.' +config_start_time = "'file'" +config_stop_time = "'2300-01-01_00:00:00'" +config_grounded_von_Mises_threshold_stress = 1.0e9 +config_min_adaptive_timestep = 21600 +config_calving_error_threshold = 1.0e9 +config_front_mass_bal_grounded = "'ismip6'" +config_use_3d_thermal_forcing_for_face_melt = '.true.' From 17e6f4e521c182697621444d40c3d0dde038e6fb Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 2 Feb 2024 19:54:57 -0800 Subject: [PATCH 29/33] add note about getting customized job script to be written --- .../tests/ensemble_generator/branch_ensemble/branch_run.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py index cf0a30e76c..864a751ff0 100644 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/branch_run.py @@ -153,6 +153,8 @@ def setup(self): self.add_model_as_input() # set job name to run number so it will get set in batch script + # Note: currently, for this to work right, one has to delete/comment + # the call to write_job_script at line 316-7 in compass/setup.py self.config.set('job', 'job_name', f'uq_{self.name}') machine = self.config.get('deploy', 'machine') pre_run_cmd = ('LOGDIR=previous_logs_`date +"%Y-%m-%d_%H-%M-%S"`;' From 3c85883977898dc792bed90a0f592e8cd8f36f29 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Sat, 3 Feb 2024 15:01:11 -0800 Subject: [PATCH 30/33] update docs --- docs/developers_guide/landice/api.rst | 4 +- .../test_groups/ensemble_generator.rst | 36 +++++-- .../test_groups/ensemble_generator.rst | 93 ++++++++++++++++--- 3 files changed, 111 insertions(+), 22 deletions(-) diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst index 8565422f9f..4bbbd37376 100644 --- a/docs/developers_guide/landice/api.rst +++ b/docs/developers_guide/landice/api.rst @@ -175,8 +175,8 @@ ensemble_generator ensemble_member.EnsembleMember.setup ensemble_member.EnsembleMember.run - control_ensemble.ControlEnsemble - control_ensemble.ControlEnsemble.configure + spinup_ensemble.SpinupEnsemble + spinup_ensemble.SpinupEnsemble.configure branch_ensemble.BranchEnsemble branch_ensemble.BranchEnsemble.configure diff --git a/docs/developers_guide/landice/test_groups/ensemble_generator.rst b/docs/developers_guide/landice/test_groups/ensemble_generator.rst index 127ca4b499..9f62d40fc5 100644 --- a/docs/developers_guide/landice/test_groups/ensemble_generator.rst +++ b/docs/developers_guide/landice/test_groups/ensemble_generator.rst @@ -90,12 +90,13 @@ methods perform minimal operations. The ``run`` method submits each run in the ensemble as a slurm job. Eventually the ``ensemble_manager`` will be able to assess if runs need restarts and modify them to be submitted as such. -ensemble --------- +spinup_ensemble +--------------- -The :py:class:`compass.landice.tests.ensemble_generator.ensemble.Ensemble` -uses the framework described above to set up an ensemble of -simulations. The constructor simply adds the ensemble manager as the only step. +The :py:class:`compass.landice.tests.ensemble_generator.spinup_ensemble.SpinupEnsemble` +uses the framework described above to set up an ensemble of spinup or historical +simulations from a common initial condition. +The constructor simply adds the ensemble manager as the only step. This allows the test case to be listed by ``compass list`` without having all ensemble members listed in a verbose listing. Because there may be dozens of ensemble members, it is better to wait to have them added until the setup @@ -112,7 +113,9 @@ Finally, each run is now added to the test case as a step to run, because they were not automatically added by compass during the test case constructor phase. -There are no ``run`` or ``validate`` steps required. The ensemble manager +The ``run`` step simply sets up a graph file and runs the model. + +The ensemble manager handles "running" ensemble members by submitting them as slurm jobs. This is a major difference in how this test case functions from most compass test cases. @@ -120,3 +123,24 @@ The visualization script ``plot_ensemble.py`` is symlinked in the test case work directory and can be run manually to assess the status of the ensemble, but there is not a formal analysis step that can be run through compass. + +branch_ensemble +--------------- + +The :py:class:`compass.landice.tests.ensemble_generator.branch_ensemble.BranchEnsemble` +sets up an ensemble of runs each of which are branched from an ensemble +member of a previously run spinup ensemble. +The constructor adds the ensemble_manager as a step, as with the spinup_ensemble. + +The ``configure`` method searches over the range of runs requested and assesses if +the corresponding spinup_ensemble member reached the requested branch time. +If so, and if the branch_ensemble memebr directory does not already exist, that +run is added as a step. Within each run (step), the restart file from the branch +year is copied to the branch run directory. The time stamp is reassigned to +2015 (this could be made a cfg option in the future). Also copied over are +the namelist and albany_input.yamlm files. The namelist is updated with +settings specific to the branch ensemble, and a streams file specific to the +branch run is added. Finally, details for managing runs are set up, including +a job script. + +As in the spinup_ensemble, the ``run`` step just runs the model. diff --git a/docs/users_guide/landice/test_groups/ensemble_generator.rst b/docs/users_guide/landice/test_groups/ensemble_generator.rst index 4dadc126d3..06dbac22ea 100644 --- a/docs/users_guide/landice/test_groups/ensemble_generator.rst +++ b/docs/users_guide/landice/test_groups/ensemble_generator.rst @@ -3,7 +3,7 @@ ensemble_generator ================== -The ``landice/ensemble_generator`` test group creates ensemble of MALI +The ``landice/ensemble_generator`` test group creates ensembles of MALI simulations with different parameter values. The ensemble framework sets up a user-defined number of simulations with parameter values selected from either uniform sampling or a space-filling Sobol sequence. @@ -49,7 +49,8 @@ Additional parameters can be easily added in the future. ``compass setup`` will set up the simulations and the ensemble manager. ``compass run`` from the test case work directory will submit each run as a separate slurm job. -Individual runs can be run independently through ``compass run`` executed in the +Individual runs can be run independently through the job script or +``compass run`` executed in the run directory. (E.g., if you want to test or debug a run without running the entire ensemble.) @@ -60,17 +61,26 @@ target year. The visualization script plots a small number of quantities of interest as a function of each active parameter. It also plots pairwise parameter sensitivities for each pair of parameters being varied. Finally, it plots time-series plots for the quantities of interest for all runs in the -ensemble. +ensemble. There are a number of options in the visualization script that +can be manually set near the top of script. Future improvements may include: * enabling the ensemble manager to identify runs that need to be restarted - so the restarts do not need to be managed manually + so the restarts in the spinup_ensemble do not need to be managed manually * safety checks or warnings before submitting ensembles that will use large amounts of computing resources -The test group includes a single test case for creating an ensemble. +The test group includes two test cases: + +* ``spinup_ensemble``: a set of simulations from the same initial condition + but with different parameter values. This could either be fixed climate + relaxation spinup or forced by time-evolving historical conditions. + +* ``branch_ensemble``: a set of simulations branched from each member of the + spinup_ensemble in a specified year with a different forcing. Multiple + branch ensembles can be branched from one spinup_ensemble config options -------------- @@ -239,16 +249,24 @@ jobs for each ensemble member. [job] wall_time = 1:30:00 -ensemble --------- +Note that currently there is not functionality +to automatically enable restart settings if runs in the spinup_ensemble +do not reach the desired year. This could be added in the future, but to +date it has been practical to set ``wall_time`` long enough to ensure this +is not a problem. Runs in a branch_ensemble are set as restarts from the +spinup_ensemble runs, so there is no need to change settings if runs +need to be continued beyond the first job. + +spinup_ensemble +--------------- -``landice/ensemble_generator/ensemble`` uses the ensemble framework to create -and ensemble of simulations integrated from 2000 to 2100. The test case +``landice/ensemble_generator/spinup_ensemble`` uses the ensemble framework to create +an ensemble of simulations integrated over a specified time range. The test case can be applied to any domain and set of input files. If the default namelist and streams settings are not appropriate, they can be adjusted or a new test case can be set up mirroring the existing one. -The model configuration uses: +The default model configuration uses: * first-order velocity solver @@ -264,11 +282,53 @@ The model configuration uses: The initial condition and forcing files are specified in the ``ensemble_generator.cfg`` file or a user modification of it. +branch_ensemble +--------------- + +``landice/ensemble_generator/branch_ensemble`` uses the ensemble framework to create +an ensemble of simulations that are branched from corresponding runs of the +``spinup_ensemble`` at a specified year with a different forcing. In general, +any namelist or streams modifications can be applied to the branch runs. + +The branch_ensemble test-case-specific config options are: + +.. code-block:: cfg + + # config options for setting up an ensemble + + # config options for branching an ensemble + [branch_ensemble] + + # start and end numbers for runs to set up and run + # branch runs. + # It is assumed that spinup runs have already been + # conducted for these runs. + start_run = 0 + end_run = 3 + + # Path to thermal forcing file for the mesh to be used in the branch run + TF_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Amery_UQ/Amery_4to20km_from_whole_AIS/forcing/ocean_thermal_forcing/UKESM1-0-LL_SSP585/1995-2300/Amery_4to20km_TF_UKESM1-0-LL_SSP585_2300.nc + + # Path to SMB forcing file for the mesh to be used in the branch run + SMB_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Amery_UQ/Amery_4to20km_from_whole_AIS/forcing/atmosphere_forcing/UKESM1-0-LL_SSP585/1995-2300/Amery_4to20km_SMB_UKESM1-0-LL_SSP585_2300_noBareLandAdvance.nc + + # location of spinup ensemble to branch from + spinup_test_dir = /pscratch/sd/h/hoffman2/AMERY_corrected_forcing_6param_ensemble_2023-03-18/landice/ensemble_generator/ensemble + + # year of spinup simulation from which to branch runs + branch_year = 2050 + + # whether to only set up branch runs for filtered runs or all runs + set_up_filtered_only = True + + # path to pickle file containing filtering information generated by plot_ensemble.py + ensemble_pickle_file = None + Steps for setting up and running an ensmble ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1. With a compass conda environment set up, run, e.g., - ``compass setup -t landice/ensemble_generator/ensemble -w WORK_DIR_PATH -f USER.cfg`` + ``compass setup -t landice/ensemble_generator/spinup_ensemble -w WORK_DIR_PATH -f USER.cfg`` where ``WORK_DIR_PATH`` is a location that can store the whole ensemble (typically a scratch drive) and ``USER.cfg`` is the user-defined config described in the previous section that includes @@ -279,7 +339,7 @@ Steps for setting up and running an ensmble 2. After ``compass setup`` completes and all runs are set up, go to the ``WORK_DIR_PATH`` and change to the - ``landice/ensemble_generator/ensemble`` subdirectory. + ``landice/ensemble_generator/spinup_ensemble`` subdirectory. From there you will see subdirectories for each run, a subdirectory for the ``ensemble_manager`` and symlink to the visualization script. @@ -290,12 +350,13 @@ Steps for setting up and running an ensmble 4. Each run will have its own batch job that can be monitored with ``squeue`` or similar commands. -5. When the ensemble has completed, you can assess the result through the +5. When the ensemble has completed, or as it is progressing, + you can assess the result through the basic visualization script ``plot_ensemble.py``. The script will skip runs that are incomplete or failed, so you can run it while an ensemble is still running to assess progress. -6. If you want to add additional ensemble members, adjust +6. If you want to run add additional ensemble members, adjust ``start_run`` and ``end_run`` in your config file and redo steps 1-5. The ensemble_manager will always be set to run the most recent run numbers defined in the config when ``compass setup`` was run. @@ -304,3 +365,7 @@ Steps for setting up and running an ensmble It is also possible to run an individual run manually by changing to the run directory and submitting the job script yourself with ``sbatch``. + +Setting up and running a branch ensemble follows the same steps. Multiple +branch ensembles (e.g., with different climate forcing scenarios) can be +conducted from one spinup ensemble. From 75827dd6d30e509347d6db49710d2091c1d692f6 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 5 Feb 2024 10:04:58 -0800 Subject: [PATCH 31/33] Fix lil bug in plot script and update default output --- compass/landice/tests/ensemble_generator/plot_ensemble.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/plot_ensemble.py b/compass/landice/tests/ensemble_generator/plot_ensemble.py index f6300060d3..bd131a422c 100644 --- a/compass/landice/tests/ensemble_generator/plot_ensemble.py +++ b/compass/landice/tests/ensemble_generator/plot_ensemble.py @@ -26,7 +26,7 @@ plot_time_series = True plot_single_param_sensitivies = False plot_pairwise_param_sensitivities = False -plot_qoi_histograms = True +plot_qoi_histograms = False plot_maps = False lw = 0.5 # linewidth for ensemble plots @@ -123,7 +123,7 @@ basin = None ens_cfg = configparser.ConfigParser() # Check for presence of two possible cfg file names -ens_cfg_file1 = 'ensemble.cfg' +ens_cfg_file1 = 'spinup_ensemble.cfg' ens_cfg_file2 = 'branch_ensemble.cfg' if os.path.isfile(ens_cfg_file1): ens_cfg_file = ens_cfg_file1 From ef6ae921f79a18c52623211f800b653d210a1b5f Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 5 Feb 2024 10:26:41 -0800 Subject: [PATCH 32/33] fix syntax errors in namelist recently added --- .../branch_ensemble/namelist.landice | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice b/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice index d09cb4290f..9564610e57 100644 --- a/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice +++ b/compass/landice/tests/ensemble_generator/branch_ensemble/namelist.landice @@ -1,8 +1,8 @@ -config_do_restart = '.true.' -config_start_time = "'file'" -config_stop_time = "'2300-01-01_00:00:00'" +config_do_restart = .true. +config_start_time = 'file' +config_stop_time = '2300-01-01_00:00:00' config_grounded_von_Mises_threshold_stress = 1.0e9 config_min_adaptive_timestep = 21600 config_calving_error_threshold = 1.0e9 -config_front_mass_bal_grounded = "'ismip6'" -config_use_3d_thermal_forcing_for_face_melt = '.true.' +config_front_mass_bal_grounded = 'ismip6' +config_use_3d_thermal_forcing_for_face_melt = .true. From ded772bcec7856b4a4ba3e73a396473f011bf101 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 28 Feb 2024 14:18:41 -0800 Subject: [PATCH 33/33] Revisions from code review --- .../ensemble_generator/ensemble_generator.cfg | 4 +- .../tests/ensemble_generator/plot_ensemble.py | 10 +++- .../test_groups/ensemble_generator.rst | 51 ++++++++++--------- 3 files changed, 39 insertions(+), 26 deletions(-) diff --git a/compass/landice/tests/ensemble_generator/ensemble_generator.cfg b/compass/landice/tests/ensemble_generator/ensemble_generator.cfg index a740675e36..4cbab8b830 100644 --- a/compass/landice/tests/ensemble_generator/ensemble_generator.cfg +++ b/compass/landice/tests/ensemble_generator/ensemble_generator.cfg @@ -20,10 +20,10 @@ end_run = 3 # parameter space sampling_method = sobol -# maximum number of sample considered. +# maximum number of samples to be considered. # max_samples needs to be greater or equal to (end_run + 1) # When using uniform sampling, max_samples should equal (end_run + 1). -# When using Sobol sequence, max_samples ought to be a multiple of 2. +# When using Sobol sequence, max_samples ought to be a power of 2. # max_samples should not be changed after the first set of ensemble. # So, when using Sobol sequence, max_samples might be set larger than # (end_run + 1) if you plan to add more samples to the ensemble later. diff --git a/compass/landice/tests/ensemble_generator/plot_ensemble.py b/compass/landice/tests/ensemble_generator/plot_ensemble.py index bd131a422c..2a2cdb8924 100644 --- a/compass/landice/tests/ensemble_generator/plot_ensemble.py +++ b/compass/landice/tests/ensemble_generator/plot_ensemble.py @@ -279,7 +279,8 @@ def filter_run(): observedSurfaceVelocityY = DSinput['observedSurfaceVelocityY'].values[0, :] obsSpdUnc = DSinput['observedSurfaceVelocityUncertainty'].values[0, :] obsSpd = (observedSurfaceVelocityX**2 + observedSurfaceVelocityY**2)**0.5 -obsSpdUnc += obsSpd * 0.25 +# May want to set speed uncertainty a function of speed, e.g., +# obsSpdUnc += obsSpd * 0.25 DSinput.close() # -------------- @@ -329,8 +330,15 @@ def filter_run(): / 1.0e12 # Gt/yr GLMigFlux = f.variables['groundingLineMigrationFlux'][:] \ / 1.0e12 # Gt/yr + # Apply smoothing to GL migration flux, because it is very noisy. + # w is the width of a window over time levels to use for smoothing. + # May need to play with this to get nice looking GL flux plots. + # The convolution is a boxcar filter with width w. + # (Note: if w is larger than the length of your time series, an error + # will occur.) w = 50 GLMigFlux2 = np.convolve(GLMigFlux, np.ones(w), 'same') / w + # GLflux2 includes both the groundingLineFlux and the migration flux GLflux2 = groundingLineFlux + GLMigFlux2 # Only process qois for runs that have reached target year diff --git a/docs/users_guide/landice/test_groups/ensemble_generator.rst b/docs/users_guide/landice/test_groups/ensemble_generator.rst index 06dbac22ea..d8f77e4a4c 100644 --- a/docs/users_guide/landice/test_groups/ensemble_generator.rst +++ b/docs/users_guide/landice/test_groups/ensemble_generator.rst @@ -72,6 +72,10 @@ Future improvements may include: * safety checks or warnings before submitting ensembles that will use large amounts of computing resources +* a method for maintaining namelist, streams, and albany_input.yaml files for + different ensembles. Currently, these input files are specific to the Amery + Ice Shelf ensemble run in 2023. + The test group includes two test cases: * ``spinup_ensemble``: a set of simulations from the same initial condition @@ -95,7 +99,6 @@ The test-case-specific config options are: .. code-block:: cfg - # config options for setting up an ensemble [ensemble] # start and end numbers for runs to set up and run @@ -115,16 +118,16 @@ The test-case-specific config options are: # for a single parameter sensitivity study. It will sample uniformly across # all dimensions simultaneously, thus sampling only a small fraction of # parameter space - sampling_method = uniform + sampling_method = sobol - # maximum number of sample considered. + # maximum number of samples to be considered. # max_samples needs to be greater or equal to (end_run + 1) # When using uniform sampling, max_samples should equal (end_run + 1). - # When using Sobol sequence, max_samples ought to be a multiple of 2. + # When using Sobol sequence, max_samples ought to be a power of 2. # max_samples should not be changed after the first set of ensemble. # So, when using Sobol sequence, max_samples might be set larger than # (end_run + 1) if you plan to add more samples to the ensemble later. - max_samples = 4 + max_samples = 1024 # basin for comparing model results with observational estimates in # visualization script. @@ -132,7 +135,7 @@ The test-case-specific config options are: # If desired basin does not exist, it can be added to that dataset. # (They need not be mutually exclusive.) # If a basin is not provided, observational comparisons will not be made. - basin = None + basin = ISMIP6BasinBC # fraction of CFL-limited time step to be used by the adaptive timestepper # This value is explicitly included here to force the user to consciously @@ -154,30 +157,28 @@ The test-case-specific config options are: # Eventually this could be hard-coded to use files on the input data # server, but initially we want flexibility to experiment with different # inputs and forcings - input_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Thwaites_UQ/Thwaites_4to20km_r02_20230126/relaxation/Thwaites_4to20km_r02_20230126_withStiffness_10yrRelax.nc + input_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Amery_UQ/Amery_4to20km_from_whole_AIS/Amery.nc # the value of the friction exponent used for the calculation of muFriction # in the input file orig_fric_exp = 0.2 # Path to ISMIP6 ice-shelf basal melt parameter input file. - basal_melt_param_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Thwaites_UQ/Thwaites_4to20km_r02_20230126/forcing/basal_melt/parameterizations/Thwaites_4to20km_r02_20230126_basin_and_coeff_gamma0_DeltaT_quadratic_non_local_median.nc + basal_melt_param_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Amery_UQ/Amery_4to20km_from_whole_AIS/forcing/basal_melt/parameterizations/Amery_4to20km_basin_and_coeff_gamma0_DeltaT_quadratic_non_local_median_allBasin2.nc # Path to thermal forcing file for the mesh to be used - TF_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Thwaites_UQ/Thwaites_4to20km_r02_20230126/forcing/ocean_thermal_forcing/obs/Thwaites_4to20km_r02_20230126_obs_TF_1995-2017_8km_x_60m_no_xtime.nc + TF_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Amery_UQ/Amery_4to20km_from_whole_AIS/forcing/ocean_thermal_forcing/obs/Amery_4to20km_obs_TF_1995-2017_8km_x_60m.nc # Path to SMB forcing file for the mesh to be used - SMB_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Thwaites_UQ/Thwaites_4to20km_r02_20230126/forcing/atmosphere_forcing/RACMO_climatology_1995-2017/Thwaites_4to20km_r02_202 - 30126_RACMO2.3p2_ANT27_smb_climatology_1995-2017.nc + SMB_file_path = /global/cfs/cdirs/fanssie/MALI_projects/Amery_UQ/Amery_4to20km_from_whole_AIS/forcing/atmosphere_forcing/RACMO_climatology_1995-2017/Amery_4to20km_RACMO2.3p2_ANT27_smb_climatology_1995-2017_no_xtime_noBareLandAdvance.nc # number of tasks that each ensemble member should be run with # Eventually, compass could determine this, but we want explicit control for now - # ntasks=32 for cori ntasks = 128 # whether basal friction exponent is being varied # [unitless] - use_fric_exp = False + use_fric_exp = True # min value to vary over fric_exp_min = 0.1 # max value to vary over @@ -195,17 +196,17 @@ The test-case-specific config options are: # [unitless: 1.0=no scaling] use_stiff_scale = True # min value to vary over - stiff_scale_min = 0.5 + stiff_scale_min = 0.8 # max value to vary over - stiff_scale_max = 1.5 + stiff_scale_max = 1.2 # whether the von Mises threshold stress (sigma_max) is being varied # [units: Pa] - use_von_mises_threshold = False + use_von_mises_threshold = True # min value to vary over - von_mises_threshold_min = 100.0e3 + von_mises_threshold_min = 80.0e3 # max value to vary over - von_mises_threshold_max = 300.0e3 + von_mises_threshold_max = 180.0e3 # whether the calving speed limit is being varied # [units: km/yr] @@ -217,7 +218,7 @@ The test-case-specific config options are: # whether ocean melt parameterization coefficient is being varied # [units: m/yr] - use_gamma0 = False + use_gamma0 = True # min value to vary over gamma0_min = 9620.0 # max value to vary over @@ -225,14 +226,14 @@ The test-case-specific config options are: # whether target ice-shelf basal melt flux is being varied # [units: Gt/yr] - use_meltflux = False + use_meltflux = True # min value to vary over - meltflux_min = 90.5 + meltflux_min = 12. # max value to vary over - meltflux_max = 114.5 + meltflux_max = 58. # ice-shelf area associated with target melt rates # [units: m^2] - iceshelf_area_obs = 4411.0e6 + iceshelf_area_obs = 60654.e6 A user should copy the default config file to a user-defined config file before setting up the test case and any necessary adjustments made. @@ -369,3 +370,7 @@ directory and submitting the job script yourself with ``sbatch``. Setting up and running a branch ensemble follows the same steps. Multiple branch ensembles (e.g., with different climate forcing scenarios) can be conducted from one spinup ensemble. + +It is intended that a single MALI executable is used for all simulations in an +ensemble. If the version of the code is updated, it is up to the user to anticipate +undesirable consequences.