Robustify extrapolate.py and its usage in the friction handling

Author: matthewhoffman

compass/landice/extrapolate.py

@@ -4,7 +4,9 @@
 from netCDF4 import Dataset
 
 
-def extrapolate_variable(nc_file, var_name, extrap_method, set_value=None):
+def extrapolate_variable(nc_file, var_name, extrap_method='auto',  # noqa: C901
+                         valid_region_method='auto',
+                         set_value=None):
     """
     Function to extrapolate variable values into undefined regions
 
@@ -19,6 +21,9 @@ def extrapolate_variable(nc_file, var_name, extrap_method, set_value=None):
     extrap_method : str
         idw, min, or value method of extrapolation
 
+    valid_region_method : str
+        choice of how to define region of valid data
+
     set_value : float
         value to set variable to outside keepCellMask
         when using -v value
@@ -36,13 +41,33 @@ def extrapolate_variable(nc_file, var_name, extrap_method, set_value=None):
     xCell = dataset.variables["yCell"][:]
     yCell = dataset.variables["xCell"][:]
 
+    # Define extrap method
+    if extrap_method == 'auto':
+        if var_name in ["effectivePressure", "beta", "muFriction"]:
+            extrap_method = 'min'
+        elif var_name in ["floatingBasalMassBal"]:
+            extrap_method = 'idw'
+        else:
+            extrap_method = 'idw'
+
     # Define region of good data to extrapolate from.
-    # Different methods for different variables
-    if var_name in ["effectivePressure", "beta", "muFriction"]:
+    if valid_region_method == 'auto':
+        if var_name in ["effectivePressure", "beta", "muFriction"]:
+            valid_region_method = 'grounded_ice'
+        elif var_name in ["floatingBasalMassBal"]:
+            valid_region_method = 'floating_ice'
+        else:
+            valid_region_method = 'ice'
+
+    print(f"Start {var_name} extrapolation using {extrap_method} method, "
+          f"defining good data using {valid_region_method} method")
+
+    # Calculate seed mask baed on method for defining valid data
+    if valid_region_method == 'positive_val':
+        keepCellMask = (varValue[:] > 0.0)
+    elif valid_region_method == 'grounded_ice':
         groundedMask = (thickness > (-1028.0 / 910.0 * bed))
         keepCellMask = np.copy(groundedMask)
-        extrap_method = "min"
-
         # grow mask by one cell oceanward of GL
         for iCell in range(nCells):
             for n in range(nEdgesOnCell[iCell]):
@@ -52,12 +77,14 @@ def extrapolate_variable(nc_file, var_name, extrap_method, set_value=None):
                     continue
         # ensure zero muFriction does not get extrapolated
         keepCellMask *= (varValue > 0.0)
-    elif var_name in ["floatingBasalMassBal"]:
+    elif valid_region_method == 'floating_ice':
         floatingMask = (thickness <= (-1028.0 / 910.0 * bed))
         keepCellMask = floatingMask * (varValue != 0.0)
-        extrap_method == "idw"
-    else:
+    elif valid_region_method == 'ice':
         keepCellMask = (thickness > 0.0)
+    else:
+        sys.exit("Unexpected value of valid_region_method encountered "
+                 "in landice/extrapolate.py")
 
     if keepCellMask.sum() == 0:
         sys.exit("In landice/extrapolate.py, initial keepCellMask has "
@@ -76,8 +103,6 @@ def extrapolate_variable(nc_file, var_name, extrap_method, set_value=None):
     # 5) Update mask
     # 6) go to step 1)
 
-    print("Start {} extrapolation using {} method".format(var_name,
-                                                          extrap_method))
     if extrap_method == 'value':
         varValue[np.where(np.logical_not(keepCellMask))] = float(set_value)
     else:

compass/landice/tests/ensemble_generator/ensemble_member.py

@@ -326,14 +326,18 @@ def _adjust_friction_exponent(orig_fric_exp, new_fric_exp, filename,
     uX = f.variables['uReconstructX'][0, :, -1]
     uY = f.variables['uReconstructY'][0, :, -1]
     spd = (uX**2 + uY**2)**0.5 * (60. * 60. * 24. * 365.)
+    if np.max(spd) == 0.0:
+        raise ValueError("uReconstructX/Y yielding max speed of 0")
     mu = mu * spd**(orig_fric_exp - new_fric_exp)
     # The previous step leads to infs or nans in ice-free areas.
     # Set them all to 0.0 for the extrapolation step
     mu[spd == 0.0] = 0.0
     f.variables['muFriction'][0, :] = mu[:]
     f.close()
 
-    extrapolate_variable(filename, 'muFriction', 'min')
+    extrapolate_variable(filename, 'muFriction',
+                         extrap_method='min',
+                         valid_region_method='positive_val')
 
     # now set exp in albany yaml file
     # using text replacement rather than yaml parser because yaml
@@ -384,8 +388,12 @@ def _apply_fric_sample(sample_num, sample_file, mu_opt_file,
     n_cells = len(f.dimensions['nCells'])
 
     # Set baseline mu outside of the Albany inverse region.
-    mu = 0.2 * np.ones(n_cells)
+    mu = -1.0 * np.ones(n_cells)
     mu[mpas_map[:] - 1] = np.exp(scaling * sample + log_mu_opt)
 
     f.variables['muFriction'][0, :] = mu[:]
     f.close()
+
+    extrapolate_variable(ic_file, 'muFriction',
+                         extrap_method='min',
+                         valid_region_method='positive_val')