Merge pull request #133 from dngoldberg/branch_vaf_mask

VAF masking

fenics_ice/config.py

@@ -265,6 +265,7 @@ class ErrorPropCfg(ConfigPrinter):
     Configuration related to error propagation
     """
     qoi: str = 'vaf'
+    qoi_apply_vaf_mask: bool = False
     phase_name: str = 'error_prop'
     phase_suffix: str = ''
 
@@ -487,6 +488,7 @@ class IOCfg(ConfigPrinter):
     alpha_data_file: str = None
     melt_depth_therm_data_file: str = None
     melt_max_data_file: str = None
+    vaf_mask_data_file: str = None
 
     thick_field_name: str = "thick"
     bed_field_name: str = "bed"
@@ -498,6 +500,7 @@ class IOCfg(ConfigPrinter):
     alpha_field_name: str = "alpha"
     melt_depth_therm_field_name: str = "melt_depth"
     melt_max_field_name: str = "melt_max"
+    vaf_mask_field_name: str = "vaf_mask"
 
     inversion_file: str = None
     qoi_file: str = None  # "Qval_ts.p"

fenics_ice/inout.py

@@ -497,7 +497,7 @@ def __init__(self, params):
 
         # List of fields to search for
         field_list = ["thick", "bed", "bmelt", "smb", "Bglen", "Bglenmask", "alpha", \
-                      "melt_depth_therm", "melt_max"]
+                      "melt_depth_therm", "melt_max", "vaf_mask"]
 
         # Dictionary of filenames & field names (i.e. field to get from HDF5 file)
         # Possibly equal to None for variables which have sensible defaults

fenics_ice/model.py

@@ -161,9 +161,6 @@ def init_fields_from_data(self):
          self.H_np = self.field_from_data("thick", self.M, min_val=min_thick)
 
         if self.params.melt.use_melt_parameterisation:
-
-         melt_depth_therm_const: float = -999.0
-         melt_max_const: float = -999.0
 
          if (self.params.melt.melt_depth_therm_const == -999.0 or \
           self.params.melt.melt_max_const == -999.0):
@@ -177,6 +174,9 @@ def init_fields_from_data(self):
           method='nearest')
          self.melt_max = self.field_from_data("melt_max", self.M2, default=melt_max, method='nearest')
 
+        if self.params.error_prop.qoi_apply_vaf_mask:
+         self.vaf_mask = self.field_from_data("vaf_mask", self.M2, default=1.0, method='nearest')
+
         self.H = self.H_np.copy(deepcopy=True)
         self.H.rename("thick_H", "")
         self.H_DG = Function(self.M2, name="H_DG")

fenics_ice/solver.py

@@ -130,6 +130,8 @@ def __init__(self, model, mixed_space=False, obs_sensitivity=False):
         if self.params.melt.use_melt_parameterisation:
          self.melt_depth_therm = model.melt_depth_therm
          self.melt_max = model.melt_max
+        if self.params.error_prop.qoi_apply_vaf_mask:
+         self.vaf_mask = model.vaf_mask
 
         self.set_inv_params()
 
@@ -1454,8 +1456,12 @@ def comp_Q_vaf(self, verbose=False):
 
         b_ex = conditional(bed < 0.0, 1.0, 0.0)
         HAF = ufl.Max(b_ex * (H + (rhow/rhoi)*bed) + (1-b_ex)*(H), 0.0)
-
-        Q_vaf = HAF * dx
+
+        if self.params.error_prop.qoi_apply_vaf_mask:
+         msk_ex = conditional(self.vaf_mask > 0.0, 1.0, 0.0)
+         Q_vaf = msk_ex * HAF * dx
+        else:
+         Q_vaf = HAF * dx
 
         if verbose:
             info(f"Q_vaf: {assemble(Q_vaf)}")