diff --git a/test/sediment/test_sed_slide.py b/test/sediment/test_sed_slide.py new file mode 100755 index 000000000..24bf708ce --- /dev/null +++ b/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' diff --git a/thetis/exner_eq.py b/thetis/exner_eq.py index b8b28983b..308bb1f8e 100644 --- a/thetis/exner_eq.py +++ b/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') diff --git a/thetis/options.py b/thetis/options.py index eed55c352..5b9880467 100644 --- a/thetis/options.py +++ b/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 diff --git a/thetis/sediment_model.py b/thetis/sediment_model.py index 7f76ca3e8..182bb3b85 100644 --- a/thetis/sediment_model.py +++ b/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