Merge pull request #216 from thetisproject/tracer_example

Sediment and morphodynamics model in 2D

examples/migrating_trench/bed.csv

@@ -0,0 +1,81 @@
+,x,Bathymetry
+0,0.0,6.61524310225572e-05
+1,0.2,0.00011704290464905355
+2,0.4,0.0001679223371808428
+3,0.6000000000000001,0.00021879467737886466
+4,0.8,0.0002696608201286439
+5,1.0,0.0003205208638387773
+6,1.2000000000000002,0.0003713735133016911
+7,1.4000000000000001,0.000422217820224626
+8,1.6,0.0004730518052333996
+9,1.8,0.0005238732433755796
+10,2.0,0.0005746794367063345
+11,2.2,0.000625466875872126
+12,2.4000000000000004,0.0006762313424690991
+13,2.6,0.0007269680029720688
+14,2.8000000000000003,0.0007776713405283933
+15,3.0,0.0008283353239552515
+16,3.2,0.0008789536272198742
+17,3.4000000000000004,0.0009295199322537436
+18,3.6,0.0009800283163689941
+19,3.8000000000000003,0.0010304737007652543
+20,4.0,0.00108085234898932
+21,4.2,0.0011311624109656322
+22,4.4,0.001181404542399653
+23,4.6000000000000005,0.0012315826953069777
+24,4.800000000000001,0.0012817052857363666
+25,5.0,0.0013317871253028702
+26,5.2,0.0013818527668980824
+27,5.4,0.0014319421682434325
+28,5.6000000000000005,0.0014821206766726849
+29,5.800000000000001,0.001532495080815121
+30,6.0,0.0015832390762258767
+31,6.2,0.0016346347654113763
+32,6.4,0.0016871365340177236
+33,6.6000000000000005,0.001741468013980367
+34,6.800000000000001,0.0017987733108787473
+35,7.0,0.001860844072378903
+36,7.2,0.0019304644516949204
+37,7.4,0.0020119329622704167
+38,7.6000000000000005,0.0021118515334252034
+39,7.800000000000001,0.0022403185238300603
+40,8.0,0.0024127341218373836
+41,8.200000000000001,0.002652536543595405
+42,8.4,0.0029953543568204608
+43,8.6,0.0034953005223820584
+44,8.8,0.004234424901386591
+45,9.0,0.0053364892047708164
+46,9.200000000000001,0.0069854693406525585
+47,9.4,0.009445294264367514
+48,9.600000000000001,0.0130655653699628
+49,9.8,0.01823400018411456
+50,10.0,0.02521526107363761
+51,10.200000000000001,0.03386997566933762
+52,10.4,0.04345565474922068
+53,10.600000000000001,0.05284834049036116
+54,10.8,0.06110241718986721
+55,11.0,0.06776598767639053
+56,11.200000000000001,0.07277505387935086
+57,11.4,0.07625845471203488
+58,11.600000000000001,0.07842709674020365
+59,11.8,0.07951686480449995
+60,12.0,0.07975477604943
+61,12.200000000000001,0.07934115580894457
+62,12.4,0.07844305461712654
+63,12.600000000000001,0.07719430838610628
+64,12.8,0.07569931544362198
+65,13.0,0.07404176248232183
+66,13.200000000000001,0.07227488814782924
+67,13.4,0.07044922641649917
+68,13.600000000000001,0.06859565852041681
+69,13.8,0.06673868879411395
+70,14.0,0.06489560068127793
+71,14.200000000000001,0.06307868527607642
+72,14.4,0.06129648069555598
+73,14.600000000000001,0.05955473307317007
+74,14.8,0.057857116049952605
+75,15.0,0.056205786000608514
+76,15.200000000000001,0.05460174515382923
+77,15.4,0.05304515550019863
+78,15.600000000000001,0.0515358199359191
+79,15.8,0.050074327868421774

examples/migrating_trench/experimental_data.csv

@@ -0,0 +1,31 @@
+5.078050288,0.000479042
+6.601538937,0.000838323
+8.564198385,-0.006347305
+8.807985994,-0.011017964
+9.052067854,-0.01497006
+9.31998411,-0.020718563
+9.541261457,-0.020359281
+9.562888964,-0.02754491
+9.732377996,-0.033652695
+9.999264371,-0.041916168
+10.29027939,-0.051257485
+10.4602098,-0.056287425
+10.75093057,-0.066347305
+10.94572526,-0.070658683
+11.01737557,-0.075688623
+11.16494284,-0.075329341
+11.48244052,-0.08
+11.75285792,-0.079640719
+11.95000662,-0.078203593
+12.19673675,-0.075688623
+12.51673557,-0.074251497
+12.76243582,-0.074251497
+12.95958452,-0.072814371
+13.20572614,-0.071736527
+13.50130206,-0.06994012
+13.69815651,-0.069221557
+14.01859671,-0.066706587
+14.21589254,-0.06491018
+14.46247554,-0.062754491
+14.75775721,-0.061676647
+14.95417028,-0.062035928

examples/migrating_trench/trench_example.py

@@ -0,0 +1,168 @@
+"""
+Migrating Trench Test case
+=======================
+
+Solves the test case of a migrating trench.
+
+We use this test case to validate the implementation of the mathematical and
+numerical methods used in Thetis to model sediment transport and morphological changes.
+In the figure produced, we compare our results with experimental data from a lab study
+
+For more details, see
+[1] Clare et al. 2020. “Hydro-morphodynamics 2D Modelling Using a Discontinuous
+    Galerkin Discretisation.” EarthArXiv. January 9. doi:10.31223/osf.io/tpqvy.
+"""
+
+from thetis import *
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+conservative = False
+
+# Note it is necessary to run trench_hydro first to get the hydrodynamics simulation
+
+# define mesh
+lx = 16
+ly = 1.1
+nx = lx*5
+ny = 5
+mesh2d = RectangleMesh(nx, ny, lx, ly)
+
+x, y = SpatialCoordinate(mesh2d)
+
+# define function spaces
+V = FunctionSpace(mesh2d, "CG", 1)
+DG_2d = FunctionSpace(mesh2d, "DG", 1)
+vector_dg = VectorFunctionSpace(mesh2d, "DG", 1)
+
+# define underlying bathymetry
+bathymetry_2d = Function(V, name='bathymetry_2d')
+initialdepth = Constant(0.397)
+depth_riv = Constant(initialdepth - 0.397)
+depth_trench = Constant(depth_riv - 0.15)
+depth_diff = depth_trench - depth_riv
+
+trench = conditional(le(x, 5), depth_riv, conditional(le(x, 6.5), (1/1.5)*depth_diff*(x-6.5) + depth_trench,
+                     conditional(le(x, 9.5), depth_trench, conditional(le(x, 11), -(1/1.5)*depth_diff*(x-11) + depth_riv, depth_riv))))
+bathymetry_2d.interpolate(-trench)
+
+# choose directory to output results
+outputdir = 'outputs'
+
+print_output('Exporting to '+outputdir)
+
+morfac = 100
+dt = 0.3
+end_time = 15*3600
+
+diffusivity = 0.15
+viscosity_hydro = Constant(1e-6)
+
+if os.getenv('THETIS_REGRESSION_TEST') is not None:
+    # the example is being run as a test
+    # run the spin-up by importing it
+    import trench_hydro  # NOQA
+    end_time = 3600.
+
+# initialise velocity and elevation
+chk = DumbCheckpoint("hydrodynamics_trench/elevation", mode=FILE_READ)
+elev = Function(DG_2d, name="elevation")
+chk.load(elev)
+chk.close()
+
+chk = DumbCheckpoint('hydrodynamics_trench/velocity', mode=FILE_READ)
+uv = Function(vector_dg, name="velocity")
+chk.load(uv)
+chk.close()
+
+# set up solver
+solver_obj = solver2d.FlowSolver2d(mesh2d, bathymetry_2d)
+options = solver_obj.options
+
+options.sediment_model_options.solve_suspended_sediment = True
+options.sediment_model_options.use_bedload = True
+options.sediment_model_options.solve_exner = True
+
+options.sediment_model_options.use_sediment_conservative_form = conservative
+options.sediment_model_options.average_sediment_size = 160*(10**(-6))
+options.sediment_model_options.bed_reference_height = 0.025
+options.sediment_model_options.morphological_acceleration_factor = Constant(morfac)
+
+options.simulation_end_time = end_time/morfac
+options.simulation_export_time = options.simulation_end_time/45
+
+options.output_directory = outputdir
+options.check_volume_conservation_2d = True
+
+if options.sediment_model_options.solve_suspended_sediment:
+    options.fields_to_export = ['sediment_2d', 'uv_2d', 'elev_2d', 'bathymetry_2d']  # note exporting bathymetry must be done through export func
+    options.sediment_model_options.check_sediment_conservation = True
+else:
+    options.fields_to_export = ['uv_2d', 'elev_2d', 'bathymetry_2d']  # note exporting bathymetry must be done through export func
+
+# using nikuradse friction
+options.nikuradse_bed_roughness = Constant(3*options.sediment_model_options.average_sediment_size)
+
+# set horizontal diffusivity parameter
+options.horizontal_diffusivity = Constant(diffusivity)
+options.horizontal_viscosity = Constant(viscosity_hydro)
+
+# 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.norm_smoother = Constant(0.1)
+
+if not hasattr(options.timestepper_options, 'use_automatic_timestep'):
+    options.timestep = dt
+
+# set boundary conditions
+
+left_bnd_id = 1
+right_bnd_id = 2
+
+swe_bnd = {}
+
+swe_bnd[left_bnd_id] = {'flux': Constant(-0.22)}
+swe_bnd[right_bnd_id] = {'elev': Constant(0.397)}
+
+solver_obj.bnd_functions['shallow_water'] = swe_bnd
+
+if options.sediment_model_options.solve_suspended_sediment:
+    # setting an equilibrium boundary conditions results in the sediment value at the boundary
+    # being chosen so that erosion and deposition are equal here (ie. in equilibrium) and the bed is immobile at this boundary
+    solver_obj.bnd_functions['sediment'] = {
+        left_bnd_id: {'flux': Constant(-0.22), 'equilibrium': None},
+        right_bnd_id: {'elev': Constant(0.397)}}
+
+    # set initial conditions
+    solver_obj.assign_initial_conditions(uv=uv, elev=elev)
+
+else:
+    # set initial conditions
+    solver_obj.assign_initial_conditions(uv=uv, elev=elev)
+
+# run model
+solver_obj.iterate()
+
+# record final bathymetry for plotting
+xaxisthetis1 = []
+baththetis1 = []
+
+for i in np.linspace(0, 15.8, 80):
+    xaxisthetis1.append(i)
+    if conservative:
+        baththetis1.append(-solver_obj.fields.bathymetry_2d.at([i, 0.55]))
+    else:
+        baththetis1.append(-solver_obj.fields.bathymetry_2d.at([i, 0.55]))
+
+if os.getenv('THETIS_REGRESSION_TEST') is None:
+    # Compare model and experimental results
+    # (this part is skipped when run as a test)
+    data = np.genfromtxt('experimental_data.csv', delimiter=',')
+
+    plt.scatter([i[0] for i in data], [i[1] for i in data], label='Experimental Data')
+
+    plt.plot(xaxisthetis1, baththetis1, label='Thetis')
+    plt.legend()
+    plt.show()