Merge pull request #105 from dngoldberg/cornford_sliding_law

Cornford sliding law

Author: dngoldberg

fenics_ice/config.py

@@ -389,7 +389,7 @@ class IceDynamicsCfg(ConfigPrinter):
 
     def __post_init__(self):
         """Check options valid"""
-        assert self.sliding_law in ['linear', 'budd']
+        assert self.sliding_law in ['linear', 'budd', 'corn']
         if self.min_thickness is not None:
             assert self.min_thickness >= 0.0
 
@@ -400,6 +400,7 @@ class MomsolveCfg(ConfigPrinter):
     Configuration of MomentumSolver with sensible defaults for picard & newton params
     """
 
+    quadrature_degree: int = -1
     picard_params: dict = field(default_factory=lambda: {
         'nonlinear_solver': 'newton',
         'newton_solver': {'linear_solver': 'cg',

fenics_ice/model.py

@@ -92,6 +92,7 @@ def __init__(self, mesh_in, input_data, param_in, init_fields=True,
             self.H_DG = None
             self.H_np = None
             self.surf = None
+            self.Umag_DG = None
 
         if init_vel_obs:
             # Load the velocity observations
@@ -145,7 +146,7 @@ def init_fields_from_data(self):
         self.H_DG = Function(self.M2, name="H_DG")
         self.bed_DG = Function(self.M2, name="bed_DG")
         self.H_DG.assign(project(self.H,self.M2))        
-        self.bed_DG.assign(project(self.bed,self.M2))        
+        self.bed_DG.assign(project(self.bed,self.M2))
 
         self.gen_surf()  # surf = bed + thick
 
@@ -455,26 +456,40 @@ def gen_surf(self):
     def bdrag_to_alpha(self, B2):
         """Convert basal drag to alpha"""
         sl = self.params.ice_dynamics.sliding_law
+
+
         if sl == 'linear':
             alpha = sqrt(B2)
 
-        elif sl == 'budd':
+        elif sl  in ['budd','corn']:
             bed = self.bed
             H = self.H
             g = self.params.constants.g
             rhoi = self.params.constants.rhoi
             rhow = self.params.constants.rhow
+            H_flt = -rhow/rhoi * bed
             u_obs = self.u_obs_M
             v_obs = self.v_obs_M
             vel_rp = self.params.constants.vel_rp
 
             # Flotation Criterion
-            H_flt = -rhow/rhoi * bed
             fl_ex = conditional(H <= H_flt, 1.0, 0.0)
-
             N = (1-fl_ex)*(H*rhoi*g + ufl.Min(bed, 0.0)*rhow*g)
             U_mag = sqrt(u_obs**2 + v_obs**2 + vel_rp**2)
-            alpha = (1-fl_ex)*sqrt(B2 * ufl.Max(N, 0.01)**(-1.0/3.0) * U_mag**(2.0/3.0))
+            
+            if sl == 'budd':
+                alpha = (1-fl_ex)*sqrt(B2 * ufl.Max(N, 0.01)**(-1.0/3.0) * U_mag**(2.0/3.0))
+
+            elif sl == 'corn':
+
+                # the relationship between alpha and B2 is too nontrivial to "invert", and an exact
+                # solution is not sought. Rather, since the sliding law is expected to deviate from 
+                # the weertman law (B2 = alpha^2 U^(-2/3)) only within a few km of the grounding line,
+                # we initialise based on the weertman sliding law
+                alpha = (1-fl_ex)*sqrt(B2 * U_mag**(2.0/3.0))
+
+        else:
+            raise RuntimeError(f"Invalid sliding law: '{sl:s}'")
 
         return alpha
 

fenics_ice/solver.py

@@ -29,6 +29,7 @@
 import time
 import ufl
 import weakref
+from IPython import embed
 
 log = logging.getLogger("fenics_ice")
 
@@ -492,15 +493,31 @@ def sliding_law(self, alpha, U):
         if sl == 'linear':
             B2 = C
 
-        elif sl == 'budd':
-            N = (1-fl_ex)*(H*rhoi*g + ufl.Min(bed, 0.0)*rhow*g)
+        elif sl in ['budd','corn']:
             U_mag = sqrt(U[0]**2 + U[1]**2 + vel_rp**2)
-            # Need to catch N <= 0.0 here, as it's raised to
-            # 1/3 (forward) and -2/3 (adjoint)
-            N_term = ufl.conditional(N > 0.0, N ** (1.0/3.0), 0)
-            B2 = (1-fl_ex)*(C * N_term * U_mag**(-2.0/3.0))
+            N = (1-fl_ex)*(H*rhoi*g + ufl.Min(bed, 0.0)*rhow*g)
+
+            if sl == 'budd':
+                # Need to catch N <= 0.0 here, as it's raised to
+                # 1/3 (forward) and -2/3 (adjoint)
+                N_term = ufl.conditional(N > 0.0, N ** (1.0/3.0), 0)
+                B2 = (1-fl_ex)*(C * N_term * U_mag**(-2.0/3.0))
+
+            elif sl == 'corn':
+                # see eq 11 of Asay-Davis et al, Experimental design for 
+                # three interrelated marine ice sheet and ocean model intercomparison projects ... 
+                # Geosci. Model Dev., 9, 2471–2497
+
+                # Need to catch N <= 0.0 here, as it's raised to
+                # 1/3 (forward) and -2/3 (adjoint)
+                N_term = ufl.conditional(N > 0.0, N, 0)
+                denom_term = (C**3 * U_mag + (0.5 * N_term)**3)**(1.0/3.0)
+                B2 = (1-fl_ex)*(C * 0.5 * N_term * U_mag**(-2.0/3.0) / denom_term)
+            
+            else:
+                raise RuntimeError(f"Invalid sliding law: '{sl:s}'")
 
-        return B2
+            return B2
 
     def solve_mom_eq(self, annotate_flag=None):
         """Solve the momentum equation defined in def_mom_eq"""
@@ -509,14 +526,16 @@ def solve_mom_eq(self, annotate_flag=None):
 
         newton_params = self.params.momsolve.newton_params
         picard_params = self.params.momsolve.picard_params
+        quad_degree = self.params.momsolve.quadrature_degree
         J_p = self.mom_Jac_p
 
         momsolver = MomentumSolver(self.mom_F == 0,
                                    self.U,
                                    bcs=self.flow_bcs,
                                    J_p=J_p,
                                    picard_params=picard_params,
-                                   solver_parameters=newton_params)
+                                   solver_parameters=newton_params,
+                                   form_compiler_parameters=None if quad_degree == -1 else {"quadrature_degree": quad_degree})
 
         momsolver.solve(annotate=annotate_flag)