Merge pull request #239 from thetisproject/sed_slide

Sed slide

test/sediment/test_sed_slide.py

@@ -0,0 +1,89 @@
+"""
+Unphysical Slope Test case
+=======================
+
+Tests the sediment slide mechanism in the Exner equation, by starting with an
+unphysical slope and reducing the slope angle over time
+
+"""
+
+from thetis import *
+
+# define mesh
+mesh2d = RectangleMesh(20, 10, 4, 2)
+x, y = SpatialCoordinate(mesh2d)
+
+vectorP1_2d = VectorFunctionSpace(mesh2d, 'DG', 1)
+V = FunctionSpace(mesh2d, 'CG', 1)
+
+# define initial bathymetry
+bathymetry_2d = Function(V, name='Bathymetry')
+z_init = conditional(x < 2, 0, conditional(x <= 4, 0.5*x-1, 0))
+bathymetry_2d.interpolate(z_init)
+
+# define initial conditions
+uv_init = Function(vectorP1_2d).interpolate(as_vector((Constant(0.46), Constant(0.0))))
+elev_init = Constant(4)
+
+# set up solver
+solver_obj = solver2d.FlowSolver2d(mesh2d, bathymetry_2d)
+options = solver_obj.options
+options.simulation_export_time = 1
+options.simulation_end_time = 20
+options.no_exports = True
+
+options.horizontal_viscosity = Constant(1e-6)
+
+# for the test, only using bedload with sediment slide mechanism
+options.sediment_model_options.solve_suspended_sediment = False
+options.sediment_model_options.use_bedload = True
+options.sediment_model_options.use_slope_mag_correction = False
+options.sediment_model_options.use_angle_correction = False
+options.sediment_model_options.use_sediment_slide = True
+options.sediment_model_options.solve_exner = True
+options.sediment_model_options.average_sediment_size = Constant(2.6e-4)
+options.sediment_model_options.bed_reference_height = Constant(0.0002)
+# average meshgrid stepsize
+options.sediment_model_options.sed_slide_length_scale = Constant(0.2)
+# maximum angle of repose which the slope should have (this is the target angle)
+options.sediment_model_options.max_angle = Constant(22)
+options.sediment_model_options.morphological_acceleration_factor = Constant(20)
+options.sediment_model_options.use_advective_velocity_correction = False
+# using nikuradse friction
+options.nikuradse_bed_roughness = Constant(3*options.sediment_model_options.average_sediment_size)
+
+# crank-nicholson used to integrate in time system of ODEs resulting from application of galerkin FEM
+options.timestepper_type = 'CrankNicolson'
+options.timestepper_options.implicitness_theta = 1.0
+options.timestep = 0.1
+
+# set boundary conditions
+left_bnd_id = 1
+right_bnd_id = 2
+
+swe_bnd = {}
+uv_vector = as_vector((0.46, 0.0))
+swe_bnd[left_bnd_id] = {'uv': uv_vector}
+swe_bnd[right_bnd_id] = {'elev': Constant(4)}
+solver_obj.bnd_functions['shallow_water'] = swe_bnd
+
+# set initial conditions
+solver_obj.assign_initial_conditions(uv=uv_init, elev=elev_init)
+
+beta = Function(V)
+max_beta_list = []
+
+
+def update_forcing(t_new):
+    # record maximum slope angle and check it is decreasing
+    beta.interpolate(solver_obj.sediment_model.betaangle)
+    max_beta_list.append(max(beta.dat.data[:])*180/pi)
+
+    if len(max_beta_list) > 30:
+        assert max_beta_list[-1] < max_beta_list[-10], 'Sediment slide mechanism is not causing\
+                                                         the angle to decrease'
+
+
+solver_obj.iterate(update_forcings=update_forcing)
+
+assert np.round(max_beta_list[-1], 1) == 24.6, 'Sediment slide mechanism has changed'

thetis/exner_eq.py

@@ -117,9 +117,12 @@ def residual(self, solution, solution_old, fields, fields_old, bnd_conditions=No
             keys = [*bnd_conditions[bnd_marker].keys()]
             values = [*bnd_conditions[bnd_marker].values()]
             for i in range(len(keys)):
-                if keys[i] != 'elev' and float(values[i]) == 0.0:
-                    no_contr = True
-
+                if keys[i] not in ('elev', 'uv'):
+                    if float(values[i]) == 0.0:
+                        no_contr = True
+                elif keys[i] == 'uv':
+                    if all(j == 0.0 for j in [float(j) for j in values[i]]):
+                        no_contr = True
             if not no_contr:
                 f += -self.test*(fac*qbx*self.n[0] + fac*qby*self.n[1])*self.ds(bnd_marker)
 
@@ -128,6 +131,26 @@ def residual(self, solution, solution_old, fields, fields_old, bnd_conditions=No
         return -f
 
 
+class ExnerSedimentSlideTerm(ExnerTerm):
+    r"""
+    Term which adds component to bedload transport to ensure the slope angle does not exceed a certain value
+    """
+    def residual(self, solution, solution_old, fields, fields_old, bnd_conditions=None):
+        f = 0
+
+        diff_tensor = self.sediment_model.get_sediment_slide_term(solution)
+
+        diff_flux = dot(diff_tensor, grad(-solution))
+        f += inner(grad(self.test), diff_flux)*dx
+        f += -avg(self.sediment_model.sigma)*inner(jump(self.test, self.sediment_model.n),
+                                                   dot(avg(diff_tensor), jump(solution,
+                                                                              self.sediment_model.n)))*dS
+        f += -inner(avg(dot(diff_tensor, grad(self.test))), jump(solution, self.sediment_model.n))*dS
+        f += -inner(jump(self.test, self.sediment_model.n), avg(dot(diff_tensor, grad(solution))))*dS
+
+        return -f
+
+
 class ExnerEquation(Equation):
     """
     Exner equation
@@ -151,3 +174,5 @@ def __init__(self, function_space, depth, sediment_model, depth_integrated_sedim
             self.add_term(ExnerSourceTerm(*args), 'source')
         if sediment_model.use_bedload:
             self.add_term(ExnerBedloadTerm(*args), 'implicit')
+        if sediment_model.use_sediment_slide:
+            self.add_term(ExnerSedimentSlideTerm(*args), 'implicit')

thetis/options.py

@@ -521,10 +521,15 @@ class SedimentModelOptions(FrozenHasTraits):
     solve_suspended_sediment = Bool(False, help='Solve suspended sediment transport equation').tag(config=True)
     use_sediment_conservative_form = Bool(False, help='Solve 2D sediment transport in the conservative form').tag(config=True)
     use_bedload = Bool(False, help='Use bedload transport in sediment model').tag(config=True)
+    use_sediment_slide = Bool(False, help='Use sediment slide mechanism in sediment model').tag(config=True)
     use_angle_correction = Bool(True, help='Switch to use slope effect angle correction').tag(config=True)
     use_slope_mag_correction = Bool(True, help='Switch to use slope effect magnitude correction').tag(config=True)
     use_secondary_current = Bool(False, help='Switch to use secondary current for helical flow effect').tag(config=True)
     average_sediment_size = FiredrakeScalarExpression(None, allow_none=True, help='Average sediment size').tag(config=True)
+    slide_region = FiredrakeScalarExpression(None, allow_none=True, help="""Region where sediment slide occurs.
+
+                                             If None then sediment slide is applied over whole domain.
+                                             """).tag(config=True)
     bed_reference_height = FiredrakeScalarExpression(None, allow_none=True, help='Bottom bed reference height').tag(config=True)
     use_advective_velocity_correction = Bool(True, help="""
         Switch to apply correction to advective velocity used in sediment equation
@@ -534,6 +539,13 @@ class SedimentModelOptions(FrozenHasTraits):
         """).tag(config=True)
     porosity = FiredrakeCoefficient(
         Constant(0.4), help="Bed porosity for exner equation").tag(config=True)
+    max_angle = FiredrakeConstantTraitlet(
+        Constant(32), help="Angle of repose for sediment slide mechanism in degrees").tag(config=True)
+    sed_slide_length_scale = FiredrakeConstantTraitlet(
+        Constant(0), help="""Length scale for sediment slide mechanism.
+
+                           This should normally be the average meshgrid size.
+                           """).tag(config=True)
     morphological_acceleration_factor = FiredrakeConstantTraitlet(
         Constant(1), help="""Rate at which timestep in exner equation is accelerated compared to timestep for model
 

thetis/sediment_model.py

@@ -85,11 +85,14 @@ def __init__(self, options, mesh2d, uv, elev, depth):
         self.options = options
         self.solve_suspended_sediment = options.sediment_model_options.solve_suspended_sediment
         self.use_bedload = options.sediment_model_options.use_bedload
+        self.use_sediment_slide = options.sediment_model_options.use_sediment_slide
         self.use_angle_correction = options.sediment_model_options.use_angle_correction
         self.use_slope_mag_correction = options.sediment_model_options.use_slope_mag_correction
         self.use_advective_velocity_correction = options.sediment_model_options.use_advective_velocity_correction
         self.use_secondary_current = options.sediment_model_options.use_secondary_current
 
+        self.mesh2d = mesh2d
+
         if not self.use_bedload:
             if self.use_angle_correction:
                 warning('Slope effect angle correction only applies to bedload transport which is not used in this simulation')
@@ -105,8 +108,11 @@ def __init__(self, options, mesh2d, uv, elev, depth):
         # define function spaces
         self.P1DG_2d = get_functionspace(mesh2d, "DG", 1)
         self.P1_2d = get_functionspace(mesh2d, "CG", 1)
+        self.R_1d = get_functionspace(mesh2d, "R", 0)
         self.P1v_2d = VectorFunctionSpace(mesh2d, "CG", 1)
 
+        self.n = FacetNormal(mesh2d)
+
         # define parameters
         self.g = physical_constants['g_grav']
         self.rhow = physical_constants['rho0']
@@ -222,6 +228,10 @@ def get_bedload_term(self, bathymetry):
             the bedlevel (positive downwards).
         """
 
+        # define bed gradient
+        dzdx = self.old_bathymetry_2d.dx(0)
+        dzdy = self.old_bathymetry_2d.dx(1)
+
         if self.use_slope_mag_correction:
             # slope effect magnitude correction due to gravity where beta is a parameter normally set to 1.3
             # we use z_n1 and equals so that we can use an implicit method in Exner
@@ -234,10 +244,6 @@ def get_bedload_term(self, bathymetry):
             cparam = Function(self.P1_2d).interpolate((self.rhos-self.rhow)*self.g*self.average_size*(self.surbeta2**2))
             tt1 = conditional(self.stress > Constant(1e-10), sqrt(cparam/self.stress), sqrt(cparam/Constant(1e-10)))
 
-            # define bed gradient
-            dzdx = self.old_bathymetry_2d.dx(0)
-            dzdy = self.old_bathymetry_2d.dx(1)
-
             # add on a factor of the bed gradient to the normal
             aa = self.salfa + tt1*dzdy
             bb = self.calfa + tt1*dzdx
@@ -303,6 +309,50 @@ def get_bedload_term(self, bathymetry):
 
         return qbx, qby
 
+    def get_sediment_slide_term(self, bathymetry):
+        # add component to bedload transport to ensure the slope angle does not exceed a certain value
+
+        # maximum gradient allowed by sediment slide mechanism
+        self.tanphi = tan(self.options.sediment_model_options.max_angle*pi/180)
+        # approximate mesh step size for sediment slide mechanism
+        L = self.options.sediment_model_options.sed_slide_length_scale
+
+        degree_h = self.P1_2d.ufl_element().degree()
+
+        if degree_h == 0:
+            self.sigma = 1.5 / CellSize(self.mesh2d)
+        else:
+            self.sigma = 5.0*degree_h*(degree_h + 1)/CellSize(self.mesh2d)
+
+        # define bed gradient
+        x, y = SpatialCoordinate(self.mesh2d)
+
+        if self.options.sediment_model_options.slide_region is not None:
+            dzdx = self.options.sediment_model_options.slide_region*bathymetry.dx(0)
+            dzdy = self.options.sediment_model_options.slide_region*bathymetry.dx(1)
+        else:
+            dzdx = bathymetry.dx(0)
+            dzdy = bathymetry.dx(1)
+
+        # calculate normal to the bed
+        nz = 1/sqrt(1 + (dzdx**2 + dzdy**2))
+
+        self.betaangle = asin(sqrt(1 - (nz**2)))
+        self.tanbeta = sqrt(1 - (nz**2))/nz
+
+        morfac = self.options.sediment_model_options.morphological_acceleration_factor
+
+        # calculating magnitude of added component
+        qaval = conditional(self.tanbeta - self.tanphi > 0, (1-self.options.sediment_model_options.porosity)
+                            * 0.5*(L**2)*(self.tanbeta - self.tanphi)/(cos(self.betaangle*self.options.timestep
+                                                                           * morfac)), 0)
+        # multiplying by direction
+        alphaconst = conditional(sqrt(1 - (nz**2)) > 0, - qaval*(nz**2)/sqrt(1 - (nz**2)), 0)
+
+        diff_tensor = as_matrix([[alphaconst, 0, ], [0, alphaconst, ]])
+
+        return diff_tensor
+
     def get_deposition_coefficient(self):
         """Returns coefficient :math:`C` such that :math:`C/H*sediment` is deposition term in sediment equation