add new spinup option for idealized experiments

combine1d/core/cost_function.py

@@ -56,7 +56,8 @@ def get_known_parameters(data_logger):
 
     # save guard if new control vars are added
     potential_control_vars = ['bed_h', 'surface_h', 'lambdas', 'w0_m',
-                              'height_shift_spinup', 'area_bed_h']
+                              'height_shift_spinup', 'area_bed_h',
+                              'section']
 
     # extract the ice free parts of variables of the control_vars which are
     # assumed to be known
@@ -76,6 +77,12 @@ def get_known_parameters(data_logger):
         elif con_var in ['surface_h']:
             prefix_var = ''
             known_index = np.full(ice_mask.shape, False)
+        elif con_var in ['section']:
+            prefix_var = ''
+            known_index = np.full(ice_mask.shape, True)
+            last_point = sum(ice_mask) + \
+                data_logger.spinup_options['section']['extra_grid_points']
+            known_index[:last_point] = False
         elif con_var in ['height_shift_spinup']:
             # no known variable to save
             continue
@@ -102,6 +109,15 @@ def get_indices_for_unknown_parameters(data_logger):
             parameter_length = ice_grid_points
         elif control_var in ['surface_h']:
             parameter_length = grid_points
+        elif control_var in ['section']:
+            extra_grid_points = \
+                data_logger.spinup_options['section']['extra_grid_points']
+            if extra_grid_points is None:
+                parameter_length = grid_points
+            else:
+                tmp_length = ice_grid_points + extra_grid_points
+                parameter_length = grid_points if tmp_length > grid_points \
+                    else tmp_length
         elif control_var in ['w0_m', 'lambdas']:
             parameter_length = trapezoid_grid_points
         elif control_var in ['height_shift_spinup']:
@@ -181,6 +197,15 @@ def define_scaling_terms(data_logger):
             data_logger.reg_scaling_terms['smoothed_bed'] = \
                 data_logger.regularisation_terms['smoothed_bed'] / \
                 smoothing_scale
+        elif regularisation_term == 'smoothed_flux':
+            data_logger.reg_scaling_terms['smoothed_flux'] = \
+                data_logger.regularisation_terms['smoothed_flux']
+        elif regularisation_term == 'distance_from_fg':
+            tmp_reg_scaling = np.ones(data_logger.len_unknown_parameter)
+            for var in data_logger.parameter_indices.keys():
+                tmp_reg_scaling[data_logger.parameter_indices[var]] = \
+                    data_logger.regularisation_terms['distance_from_fg'][var]
+            data_logger.reg_scaling_terms['distance_from_fg'] = tmp_reg_scaling
         else:
             raise NotImplementedError(f'{regularisation_term}')
 
@@ -227,6 +252,11 @@ def cost_fct(unknown_parameters, data_logger):
 
     mb_models = initialise_mb_models(unknown_parameters_descaled, data_logger)
 
+    if 'smoothed_flux' in data_logger.regularisation_terms.keys():
+        initial_flux = dynamic_model(flowline).flux_stag[0]
+    else:
+        initial_flux = None
+
     if data_logger.spinup_type == 'height_shift_spinup':
         try:
             # Here a spinup run is conducted using the control variable
@@ -242,13 +272,16 @@ def cost_fct(unknown_parameters, data_logger):
             cost = np.array(np.Inf)
             grad = np.empty(len(unknown_parameters)) * np.nan
             return cost, grad
-    elif data_logger.spinup_type not in [None, 'surface_h']:
+    elif data_logger.spinup_type not in [None, 'surface_h', 'perfect_sfc_h',
+                                         'perfect_thickness',
+                                         'perfect_section', 'section']:
         raise NotImplementedError(f'The spinup option {data_logger.spinup_type} '
                                   'is not implemented!')
 
     # sfc_h saved to be able to save the past glacier evolution after the
     # minimisation and for the evaluation of idealized experiments
     sfc_h_start = flowline.surface_h.detach().clone()
+    section_start = flowline.section.detach().clone()
 
     observations = data_logger.observations
 
@@ -274,6 +307,8 @@ def cost_fct(unknown_parameters, data_logger):
     c_terms, reg_terms = get_cost_terms(
         observations_mdl,
         final_fl,  # for regularisation term 'smoothed_bed'
+        initial_flux,  # for regularisation term 'smoothed_flux'
+        unknown_parameters,  # for distance from fg regularisation
         data_logger  # for reg_parameters and observations
         )
 
@@ -299,6 +334,7 @@ def cost_fct(unknown_parameters, data_logger):
     data_logger.save_data_in_datalogger(
         'observations_mdl', detach_observations_mdl(observations_mdl))
     data_logger.save_data_in_datalogger('sfc_h_start', sfc_h_start)
+    data_logger.save_data_in_datalogger('section_start', section_start)
     data_logger.save_data_in_datalogger('flowlines', detach_flowline(final_fl))
     data_logger.save_data_in_datalogger('costs', cost)
     data_logger.save_data_in_datalogger('grads', grad)
@@ -385,7 +421,8 @@ def initialise_flowline(unknown_parameters, data_logger):
 
     # initialise all potential control variables for flowline as empty tensor
     # and fill them
-    all_potential_control_vars = ['surface_h', 'lambdas', 'w0_m']
+    all_potential_control_vars = ['surface_h', 'lambdas', 'w0_m',
+                                  'section']
     if not any(k in parameter_indices.keys() for k in ['bed_h', 'area_bed_h']):
         all_potential_control_vars.insert(0, 'bed_h')
     elif all(k in parameter_indices.keys() for k in ['bed_h', 'area_bed_h']):
@@ -396,6 +433,8 @@ def initialise_flowline(unknown_parameters, data_logger):
         all_potential_control_vars.insert(0, 'bed_h')
     elif 'area_bed_h' in parameter_indices.keys():
         all_potential_control_vars.insert(0, 'area_bed_h')
+    elif 'perfect_bed_h' in data_logger.spinup_options.keys():
+        all_potential_control_vars.insert(0, 'bed_h')
     else:
         raise NotImplementedError('I have not expected to come to this point!')
     fl_vars_total = {}  # they are used for actual initialisation
@@ -418,6 +457,11 @@ def initialise_flowline(unknown_parameters, data_logger):
                 var_index = np.full(ice_mask.shape, True)
             elif var in ['lambdas', 'w0_m']:
                 var_index = (trap_index & ice_mask)
+            elif var in ['section']:
+                var_index = np.full(ice_mask.shape, False)
+                last_point = sum(ice_mask) + \
+                    data_logger.spinup_options['section']['extra_grid_points']
+                var_index[:last_point] = True
             else:
                 raise NotImplementedError(f'{var}')
 
@@ -475,10 +519,52 @@ def initialise_flowline(unknown_parameters, data_logger):
                     prefix = '_'
                 else:
                     prefix = ''
-                fl_vars_total[var] = torch.tensor(getattr(fl_init, prefix + var),
-                                                  dtype=torch_type,
-                                                  device=device,
-                                                  requires_grad=False)
+
+                if var == 'surface_h' and \
+                    data_logger.spinup_type in ['perfect_sfc_h',
+                                                'perfect_thickness',
+                                                'perfect_section',
+                                                'section']:
+                    if data_logger.spinup_type == 'perfect_sfc_h':
+                        perfect_sfc_h = torch.tensor(
+                            data_logger.perfect_spinup_value,
+                            dtype=torch_type,
+                            device=device,
+                            requires_grad=False)
+                        fl_vars_total['surface_h'] = perfect_sfc_h
+                    elif data_logger.spinup_type == 'perfect_thickness':
+                        perfect_thickness = torch.tensor(
+                            data_logger.perfect_spinup_value,
+                            dtype=torch_type,
+                            device=device,
+                            requires_grad=False)
+                    elif data_logger.spinup_type in ['perfect_section',
+                                                     'section']:
+                        fl_vars_total['surface_h'] = torch.tensor(
+                            getattr(fl_init, prefix + var) * 0,
+                            dtype=torch_type,
+                            device=device,
+                            requires_grad=False
+                        )
+                    else:
+                        raise NotImplementedError(f'{data_logger.spinup_options}')
+                elif var == 'bed_h' and \
+                        'perfect_bed_h' in data_logger.spinup_options.keys():
+                    fl_vars_total['bed_h'] = torch.tensor(
+                        data_logger.perfect_spinup_value,
+                        dtype=torch_type,
+                        device=device,
+                        requires_grad=False
+                    )
+                else:
+                    fl_vars_total[var] = torch.tensor(
+                        getattr(fl_init, prefix + var),
+                        dtype=torch_type,
+                        device=device,
+                        requires_grad=False)
+
+    if data_logger.spinup_type == 'perfect_thickness':
+        fl_vars_total['surface_h'] = perfect_thickness + fl_vars_total['bed_h']
 
     fl = MixedBedFlowline(line=fl_init.line,
                           dx=fl_init.dx,
@@ -495,6 +581,16 @@ def initialise_flowline(unknown_parameters, data_logger):
                           torch_type=torch_type,
                           device=device)
 
+    if data_logger.spinup_type == 'perfect_section':
+        fl.section = torch.tensor(
+            data_logger.perfect_spinup_value,
+            dtype=torch_type,
+            device=device,
+            requires_grad=False
+        )
+    elif data_logger.spinup_type == 'section':
+        fl.section = fl_vars_total['section']
+
     return fl
 
 
@@ -581,6 +677,8 @@ def do_height_shift_spinup(flowline, unknown_parameters, data_logger):
 
 def get_cost_terms(observations_mdl,
                    flowline,
+                   flux,
+                   unknown_parameters,
                    data_logger):
     """
     Returns the individual terms of the cost function. TODO
@@ -589,6 +687,7 @@ def get_cost_terms(observations_mdl,
     ----------
     observations_mdl
     flowline
+    flux
     data_logger
 
     Returns
@@ -669,6 +768,32 @@ def get_cost_terms(observations_mdl,
                           'bed_h_grad_scale']:
             # this regularisation terms are handled elsewhere
             pass
+        elif reg_term in ['smoothed_flux']:
+            reg_scale = torch.tensor(data_logger.reg_scaling_terms[reg_term],
+                                     dtype=data_logger.torch_type,
+                                     device=data_logger.device,
+                                     requires_grad=False)
+            reg_terms[i_costs] = reg_scale * \
+                (torch.sum(torch.abs(torch.diff(flux))) / torch.max(flux) - 2)
+
+            reg_term_np = reg_terms[i_costs].detach().to('cpu').numpy()
+            cost_terms_description[reg_term] = reg_term_np
+            individual_reg_terms.append(reg_term_np)
+        elif reg_term in ['distance_from_fg']:
+            reg_scale = torch.tensor(data_logger.reg_scaling_terms[reg_term],
+                                     dtype=data_logger.torch_type,
+                                     device=data_logger.device,
+                                     requires_grad=False)
+            first_guess = torch.tensor(data_logger.first_guess,
+                                       dtype=data_logger.torch_type,
+                                       device=data_logger.device,
+                                       requires_grad=False)
+            reg_terms[i_costs] = torch.sum(
+                reg_scale * (first_guess - unknown_parameters) /
+                torch.tensor([data_logger.len_unknown_parameter],
+                             dtype=data_logger.torch_type,
+                             device=data_logger.device,
+                             requires_grad=False))
         else:
             raise NotImplementedError(f'{reg_term}')
         i_costs += 1

combine1d/core/data_logging.py

@@ -60,11 +60,36 @@ def __init__(self, gdir, fls_init, inversion_input, climate_filename,
         elif 'surface_h' in list(self.spinup_options.keys()):
             spinup_type = 'surface_h'
             inversion_input['control_vars'].append(spinup_type)
+        elif 'section' in list(self.spinup_options.keys()):
+            spinup_type = 'section'
+            inversion_input['control_vars'].append(spinup_type)
         elif 'height_shift' in list(self.spinup_options.keys()):
             spinup_type = 'height_shift_spinup'
             inversion_input['control_vars'].append(spinup_type)
+        elif list(self.spinup_options.keys())[0] in ['perfect_sfc_h',
+                                                     'perfect_thickness',
+                                                     'perfect_section']:
+            spinup_type = list(self.spinup_options.keys())[0]
+            fls_true_init = gdir.read_pickle(
+                'model_flowlines',
+                filesuffix=self.spinup_options[spinup_type])[0]
+            if spinup_type == 'perfect_sfc_h':
+                self.perfect_spinup_value = fls_true_init.surface_h
+            elif spinup_type == 'perfect_thickness':
+                self.perfect_spinup_value = fls_true_init.thick
+            elif spinup_type == 'perfect_section':
+                self.perfect_spinup_value = fls_true_init.section
+            else:
+                raise NotImplementedError(f'{spinup_type}')
         else:
-            raise NotImplementedError
+            raise NotImplementedError(f'{self.spinup_options.keys()}')
+
+        if spinup_type is not None:
+            if 'perfect_bed_h' in list(self.spinup_options.keys()):
+                fls_true_init = gdir.read_pickle(
+                    'model_flowlines',
+                    filesuffix=self.spinup_options['perfect_bed_h'])[0]
+                self.perfect_spinup_value = fls_true_init.bed_h
         self.spinup_type = spinup_type
         self.control_vars = inversion_input['control_vars']
 
@@ -115,6 +140,7 @@ def __init__(self, gdir, fls_init, inversion_input, climate_filename,
         self.grads = None
         self.flowlines = None
         self.sfc_h_start = None
+        self.section_start = None
         self.observations_mdl = None
         self.end_time = None
         self.known_parameters = None
@@ -127,6 +153,7 @@ def __init__(self, gdir, fls_init, inversion_input, climate_filename,
         self.minimize_status = None
         self.memory_error = False
         self.extra_bed_h = None
+        self.first_guess = None
 
         self.filename = gdir.name + '_' + \
                         inversion_input['experiment_description']
@@ -175,7 +202,7 @@ def callback_fct(self, x0):
             for key in self.c_terms_description[-1][0].keys():
                 if key in ['J_obs', 'J_reg']:
                     continue
-                elif key in ['smoothed_bed']:
+                elif key in ['smoothed_bed', 'smoothed_flux']:
                     reg_terms_str.append(
                         f'{key}: {self.c_terms_description[-1][0][key]:.2e}')
                 else:
@@ -221,6 +248,7 @@ def filter_data_from_optimization(self):
         self.time_needed = self.squeeze_generic(self.time_needed[index])
         self.flowlines = self.squeeze_generic(self.flowlines[index])
         self.sfc_h_start = self.squeeze_generic(self.sfc_h_start[index])
+        self.section_start = self.squeeze_generic(self.section_start[index])
         self.observations_mdl = self.squeeze_generic(self.observations_mdl[index])
         self.unknown_parameters = self.unknown_parameters[index]
         self.fct_calls = self.squeeze_generic(self.fct_calls[index + 1])
@@ -241,6 +269,7 @@ def create_and_save_dataset(self):
         ds['grads'] = (['iteration', 'nr_unknown_parameters'], self.grads)
         ds['flowlines'] = (['iteration'], self.flowlines)
         ds['sfc_h_start'] = (['iteration', 'x'], self.sfc_h_start)
+        ds['section_start'] = (['iteration', 'x'], self.section_start)
         ds['observations_mdl'] = (['iteration'], self.observations_mdl)
         ds['c_terms'] = (['iteration', 'nr_cost_terms'], self.c_terms)
         ds['reg_terms'] = (['iteration', 'nr_reg_terms'], self.reg_terms)

combine1d/core/dynamics.py

@@ -93,6 +93,8 @@ def add_model_variable(year, variable):
                 unit_check(obs_name, obs_unit, ['m'])
             elif obs_name in ['us']:
                 unit_check(obs_name, obs_unit, ['myr-1'])
+            elif obs_name in ['volume']:
+                unit_check(obs_name, obs_unit, ['m3', 'km3'])
             else:
                 raise NotImplementedError(f'{obs_name} not implemented!')
 
@@ -148,6 +150,8 @@ def run_model_and_get_model_values(flowline, dynamic_model, mb_models,
                     actual_model_data[obs_yr][var] = dyn_model.area_km2
                 elif var == 'volume:m3':
                     actual_model_data[obs_yr][var] = dyn_model.volume_m3
+                elif var == 'volume:km3':
+                    actual_model_data[obs_yr][var] = dyn_model.volume_km3
                 elif var == 'fl_surface_h:m':
                     actual_model_data[obs_yr][var] = dyn_model.fls[0].surface_h
                 elif var == 'fl_widths:m':
@@ -198,7 +202,7 @@ def calculate_model_observations(observations, actual_model_data):
         # first all observations where nothing need to be calculated
         if var_key in ['area:m2', 'area:km2', 'fl_total_area:m2',
                        'fl_total_area:km2', 'fl_surface_h:m', 'fl_widths:m',
-                       'us:myr-1']:
+                       'us:myr-1', 'volume:m3', 'volume:km3']:
             for year in out[var_key].keys():
                 out[var_key][year] = actual_model_data[int(year)][var_key]
 

combine1d/core/first_guess.py

@@ -1,8 +1,9 @@
 import copy
 
 import numpy as np
-from oggm.core.flowline import FluxBasedModel
-from oggm.core.massbalance import MultipleFlowlineMassBalance, ConstantMassBalance
+from oggm.core.flowline import FluxBasedModel, SemiImplicitModel
+from oggm.core.massbalance import MultipleFlowlineMassBalance, \
+    ConstantMassBalance, MonthlyTIModel
 
 
 def get_first_guess(data_logger):
@@ -26,6 +27,15 @@ def get_first_guess(data_logger):
         elif ind in ['height_shift_spinup']:
             ind_first_guess = \
                 data_logger.spinup_options['height_shift']['mb_model']['fg_height_shift']
+        elif ind in ['section']:
+            extra_grid_points =\
+                data_logger.spinup_options['section']['extra_grid_points']
+            nx = sum(ice_mask)
+            dyn_model = SemiImplicitModel(
+                fl, mb_model=MonthlyTIModel(data_logger.gdir), y0=1980)
+            dyn_model.run_until(
+                1980 + data_logger.spinup_options['section']['fg_years'])
+            ind_first_guess = dyn_model.fls[0].section[:nx + extra_grid_points]
         else:
             raise NotImplementedError(f'{ind} is not implemented!')
 
@@ -38,6 +48,8 @@ def get_first_guess(data_logger):
     scale = data_logger.control_vars_characteristic_scale
     first_guess = (first_guess - min_bound) / (max_bound - min_bound) * scale
 
+    data_logger.first_guess = first_guess
+
     return first_guess