Merge branch 'master' into tsunami

demos/demo_2d_north_sea.py

@@ -39,34 +39,30 @@
 # correspond to open ocean (tag 100) or coasts (tag 200). This is because we impose
 # different boundary conditions in each case.
 #
-# We set up the (UTM) mesh and calculate its longitude-latitude coordinates as follows:
-
-mesh2d = Mesh("north_sea.msh")
-lonlat = coord_system.get_mesh_lonlat_function(mesh2d)
-lon, lat = lonlat
-
-# With the mesh, we can now move on to set up fields defined upon it. For the
-# bathymetry data, we use the
-# `ETOPO1 data set <https://www.ngdc.noaa.gov/mgg/global>`__
-# :cite:`ETOPO1:2009`, :cite:`ETOPO1tech:2009`.
-# A NetCDF file containing such data for the North Sea can be downloaded from the
-# webpage, stored as `etopo1.nc` and then interpolated onto the unstructured mesh
-# using SciPy. An `interpolate_bathymetry` script is provided in the
-# `example <https://github.com/thetisproject/thetis/tree/master/examples/north_sea>`__
-# in the Thetis source code, which follows the
-# `recommendations in the Firedrake documentation for interpolating data <https://firedrakeproject.org/interpolation.html#interpolation-from-external-data>`__.
+# Normally, we would read in this mesh using ``mesh2d = Mesh('north_sea.msh)``;
+# Here however, we use a pre-prepared bathymetry that is stored in a checkpoint
+# file and therefore also read in the mesh from the checkpoint file.
+#
+# For the bathymetry data, we used the `ETOPO1 data set
+# <https://www.ngdc.noaa.gov/mgg/global>`__ :cite:`ETOPO1:2009`,
+# :cite:`ETOPO1tech:2009`.  A NetCDF file containing such data for the North
+# Sea can be downloaded from the webpage, stored as `etopo1.nc` and then
+# interpolated onto the unstructured mesh using SciPy. An
+# `interpolate_bathymetry` function is provided in model_config module in the
+# `example
+# <https://github.com/thetisproject/thetis/tree/master/examples/north_sea>`__
+# in the Thetis source code, which follows the `recommendations in the
+# Firedrake documentation for interpolating data
+# <https://firedrakeproject.org/interpolation.html#interpolation-from-external-data>`__.
 # However, the NetCDF file cannot be included here for
 # copyright reasons, so we insteady provide a HDF5 file containing the data already
 # interpolated onto the mesh. Note that HDF5 files currently have to be saved and
 # loaded using the same number of processors. The bathymetry field was generated by
 # a serial run, so the following will not work in parallel.
 
-P1_2d = get_functionspace(mesh2d, "CG", 1)
-bathymetry_2d = Function(P1_2d, name="Bathymetry")
 with CheckpointFile("north_sea_bathymetry.h5", "r") as f:
-    m = f.load_mesh("firedrake_default")
-    g = f.load_function(m, "Bathymetry")
-    bathymetry_2d.assign(g)
+    mesh2d = f.load_mesh("firedrake_default")
+    bathymetry_2d = f.load_function(mesh2d, "Bathymetry")
 
 # The resulting bathymetry field is plotted below.
 #
@@ -94,11 +90,13 @@ def read_station_data():
 # We also require fields for the Manning friction coefficient and the Coriolis forcing.
 # These can be set up as follows:
 
+P1_2d = get_functionspace(mesh2d, "CG", 1)
 manning_2d = Function(P1_2d, name="Manning coefficient")
 manning_2d.assign(3.0e-02)
 
 omega = 7.292e-05
 coriolis_2d = Function(P1_2d, name="Coriolis forcing")
+lon, lat = coord_system.get_mesh_lonlat_function(mesh2d)
 coriolis_2d.interpolate(2 * omega * sin(lat * pi / 180.0))
 
 # We also need to choose a time window of interest and discretise it appropriately.
@@ -206,7 +204,7 @@ def read_station_data():
 #
 # For the purposes of this demo, we have included HDF5 files containing spun-up
 # elevation and velocity fields in the `outputs_spinup` directory. These can be
-# used to initialise the model as follows. Again, the spun-up HDF5 data were
+# used to initialise the model as follows. The spun-up HDF5 data were
 # generated by a serial run, so this demo will not work in parallel.
 
 solver_obj.load_state(14, outputdir="outputs_spinup", t=0, iteration=0)

docs/Makefile

@@ -26,7 +26,7 @@ copy_demos:
 	cp ../demos/*.py source/demos
 	cp ../demos/figures/*.png source/demos
 	cp ../demos/demo_references.bib source/demos
-	for file in source/demos/*.py; do ${VIRTUAL_ENV}/src/firedrake/pylit/pylit.py -c $$file; mv $$file.txt $$file.rst; done
+	for file in source/demos/*.py; do pylit -c $$file; mv $$file.txt $$file.rst; done
 	install -d $(BUILDDIR)/html/demos
 	cp source/demos/*.py $(BUILDDIR)/html/demos
 

docs/source/docutils.conf

@@ -0,0 +1,2 @@
+[html writers]
+table-style: colwidths-grid

docs/source/outputs_and_visu.rst

@@ -98,15 +98,17 @@ Restarting a simulation
 ~~~~~~~~~~~~~~~~~~~~~~~
 
 If you have stored the required HDF5 files, you can continue a simulation
-using :py:meth:`~.FlowSolver.load_state` method, provided that you use the same
-mesh and the same number of MPI processes. This call replaces the
+using :py:meth:`~.FlowSolver.load_state` method, provided that you use the
+mesh from that checkpoint file. This call replaces the
 :py:meth:`~.FlowSolver.assign_initial_conditions` call.
 If initial conditions are not set, add ``load_state`` call above
 the :py:meth:`~.FlowSolver.iterate` call.
 
-In the simplest form, one only defines the export index that is used as initial
-condition::
+In the simplest form, the mesh needs to be loaded from a HDF5 checkpointfile
+and then the index given to the :py:meth:`~.FlowSolver.load_state` method::
 
+    mesh2d = read_mesh_from_checkpoint(outputdir)
+    ...set-up the thetis solver object with mesh2d...
     solver_obj.load_state(155)
 
 This also loads simulation time from the stored state.

docs/source/publications.rst

@@ -8,9 +8,38 @@ Citing Thetis
 If you publish results obtained with Thetis, please cite the `2018 GMD paper <https://doi.org/10.5194/gmd-11-4359-2018>`_.
 Also see `instructions for citing Firedrake <https://firedrakeproject.org/citing.html>`_.
 
+2023
+----
+
+Wallwork, J. G., Angeloudis, A., Barral, N., Mackie, L., Kramer, S. C., Piggott, M. D.:
+Tidal turbine array modelling using goal-oriented mesh adaptation.
+*Journal of Ocean Engineering and Marine Energy*,
+doi: `10.1007/s40722-023-00307-9 <https://doi.org/10.1007/s40722-023-00307-9>`_, 2023.
+
+Kärnä, T., Wallwork, J. G., Kramer, S. C.:
+Adjoint-based optimization of a regional water elevation model.
+*Journal of Advances in Modeling Earth Systems*,
+doi: `10.1029/2022MS003169 <https://doi.org/10.1029/2022MS003169>`_, 2023.
+
+Woodroffe, S., Hill, J., Bustamante-Fernandez, E., Lloyd, JM., Luff, J., Richards, S., Shennan, I:
+On the varied impact of the Storegga tsunami in northwest Scotland.
+*Journal of Quaternary Science*, 
+doi: `10.1002/jqs.3539 <https://doi.org/10.1002/jqs.3539>`_, 2023.
+
+Hill, J., Rush, G., Peakall, J., Johnson, M., Hodson, L., Barlow, N.L.M.,
+Bowman, E.T., Gehrels, W.R., Hodgson, D.M., Kesserwani, G:
+Resolving tsunami wave dynamics: integrating sedimentology and numerical modelling.
+*The Depositional Record*
+doi: `10.1002/dep2.247 <https://doi.org/10.1002/dep2.247>`_, 2023.
+
 2022
 ----
 
+Lee, K.C., Webster, J.M., Salles, T., Mawson, E.E., Hill, J.:
+Tidal dynamics drive ooid formation in the Capricorn Channel since the Last Glacial Maximum.
+*Marine Geology*, 454:106944,
+doi:`10.1016/j.margeo.2022.106944 <https://doi.org/10.1016/j.margeo.2022.106944>`_, 2022.
+
 Mawson, E.E., Lee, K.C., Hill, J.:
 Sea level rise and the Great Barrier Reef: the future implications on reef tidal dynamics.
 *Journal of Geophysical Research: Oceans*, 127:e2021JC017823, 
@@ -24,13 +53,38 @@ Clare, M. C. A., Wallwork, J. G., Kramer, S. C., Weller, H., Cotter, C. J., Pigg
 Multi-scale hydro-morphodynamic modelling using mesh movement methods. *GEM - International Journal on Geomathematics*, 13:1,
 doi:`10.1007/s13137-021-00191-1 <https://doi.org/10.1007/s13137-021-00191-1>`_, 2022.
 
+Clare, M. C. A., Kramer, S. C., Cotter, C. J., Piggott, M. D.:
+Calibration, inversion and sensitivity analysis for hydro-morphodynamic models through the application of adjoint methods. 
+*Computers & Geosciences* 163, p. 105104,
+doi: `10.1016/j.cageo.2022.105104 <https://doi.org/10.1016/j.cageo.2022.105104>`_, 2022.
+
+Pan, W., Kramer, S. C., Piggott, M. D., Xiping, Y.:
+Modeling landslide generated waves using the discontinuous finite element method.
+*International Journal for Numerical Methods in Fluids*, 94 (8), pp. 1298-1330,
+doi: `10.1002/fld.5090 <https://doi.org/10.1002/fld.5090>`_, 2022.
+
+Zhang, C., Kramer, S. C.,  Angeloudis, A., Zhang, J., Lin, X., Piggott, M. D.:
+Improving tidal turbine array performance through the optimisation of layout and yaw angles.
+*International Marine Energy Journal* 5.3, pp. 273–280,
+doi: `10.36688/imej.5.273-280 <https://doi.org/10.36688/imej.5.273-280>`_, 2022
+
+Zhang, J., Zhang, C., Angeloudis, A., Kramer, S. C., He, R., Piggott, M. D.:
+Interactions between tidal stream turbine arrays and their hydrodynamic impact around Zhoushan Island, China
+*Ocean Engineering* 246, p. 110431,
+doi: `10.1016/j.oceaneng.2021.110431 <https://doi.org/10.1016/j.oceaneng.2021.110431>`_, 2022
+
 2021
 ----
 
 Fofonova, V., Kärnä, T., Klingbeil, K., Androsov, A., Kuznetsov, I., Sidorenko, D., Danilov, S., Burchard, H., Wiltshire, K. H.: Plume spreading test case for coastal ocean models. *Geosci. Model Dev.*, 14:6945–6975, doi:`10.5194/gmd-14-6945-2021 <https://doi.org/10.5194/gmd-14-6945-2021>`_, 2021.
 
 Bateman, M. D., Kinnaird, T. C., Hill, J., Ashurst, R. A., Mohan, J., Bateman, R. B. I., Robinson, R.: Detailing the impact of the Storegga Tsunami at Montrose, Scotland. *Boreas*, bor.12532, doi:`10.1111/bor.12532 <https://doi.org/10.1111/bor.12532>`_, 2021.
 
+Mackie, L., Kramer, S. C., Piggott, M. D., Angeloudis, A.:
+Assessing impacts of tidal power lagoons of a consistent design.
+*Ocean Engineering* 240, p. 109879. 
+doi: `10.1016/j.oceaneng.2021.109879 <https://doi.org/10.1016/j.oceaneng.2021.109879>`
+
 Pan, W., Kramer, S. C., and Piggott, M. D.: A sigma-coordinate non-hydrostatic discontinuous finite element coastal ocean model. *Ocean Modelling*, 157:101732, doi:`10.1016/j.ocemod.2020.101732 <https://doi.org/10.1016/j.ocemod.2020.101732>`_, 2021.
 
 Rasheed, S., Warder, S. C., Plancherel, Y., and Piggott, M. D.:
@@ -46,6 +100,11 @@ doi:`10.1016/j.cageo.2020.104658 <https://doi.org/10.1016/j.cageo.2020.104658>`_
 Warder, S. C., Horsburgh, K. J. and Piggott, M. D.:
 Adjoint-based sensitivity analysis for a numerical storm surge model. *Ocean Modelling*, 160:101766, doi:`10.1016/j.ocemod.2021.101766 <https://doi.org/10.1016/j.ocemod.2021.101766>`_, 2021.
 
+Goss, Z. L., Coles, D. S., Kramer, S. C., Piggott, M. D.:
+Efficient economic optimisation of large-scale tidal stream arrays.
+*Applied Energy* 295, p. 116975.
+doi: `10.1016/j.apenergy.2021.116975 <https://doi.org/10.1016/j.apenergy.2021.116975>`_, 2021.
+
 
 2020
 ----

examples/bottomFriction/ekman_bottom.py

@@ -21,8 +21,11 @@ def bottom_ekman_test(layers=50, verify=True, iterate=True,
     lx = nx*dx
     ny = 3
     ly = ny*dx
-    mesh2d = PeriodicRectangleMesh(nx, ny, lx, ly, direction='both',
-                                   reorder=True)
+    if load_export_ix:
+        mesh2d = read_mesh_from_checkpoint(outputdir)
+    else:
+        mesh2d = PeriodicRectangleMesh(nx, ny, lx, ly, direction='both',
+                                       reorder=True)
 
     dt = 90.0
     t_end = 5 * 3600.0  # sufficient to reach ~steady state

examples/bottomFriction/ekman_surface.py

@@ -20,8 +20,11 @@ def surface_ekman_test(layers=50, verify=True, iterate=True,
     lx = nx*dx
     ny = 3
     ly = ny*dx
-    mesh2d = PeriodicRectangleMesh(nx, ny, lx, ly, direction='both',
-                                   reorder=True)
+    if load_export_ix:
+        mesh2d = read_mesh_from_checkpoint(outputdir)
+    else:
+        mesh2d = PeriodicRectangleMesh(nx, ny, lx, ly, direction='both',
+                                       reorder=True)
 
     dt = 90.0
     t_end = 6 * 3600.0

examples/bottomFriction/steadyChannel.py

@@ -38,7 +38,10 @@ def bottom_friction_test(layers=25, gls_closure='k-omega',
     lx = nx*dx
     ny = 3  # nb elements in cross direction
     ly = ny*dx
-    mesh2d = PeriodicRectangleMesh(nx, ny, lx, ly, direction='x', reorder=True)
+    if load_export_ix:
+        mesh2d = read_mesh_from_checkpoint(outputdir)
+    else:
+        mesh2d = PeriodicRectangleMesh(nx, ny, lx, ly, direction='x', reorder=True)
 
     dt = 25.0
     t_end = 12 * 3600.0  # sufficient to reach ~steady state

examples/channel_inversion/inverse_problem.py

@@ -1,5 +1,5 @@
 from thetis import *
-from firedrake_adjoint import *
+from firedrake.adjoint import *
 import numpy
 import thetis.inversion_tools as inversion_tools
 from model_config import construct_solver
@@ -31,6 +31,9 @@
 pwd = os.path.abspath(os.path.dirname(__file__))
 output_dir = f'{pwd}/outputs_new_' + '-'.join(controls) + '-opt'
 
+# annotate all Firedrake operations of the forward run
+continue_annotation()
+
 # Create the solver object
 # NOTE We disable all exports here as they would only be active during
 # the initial forward run
@@ -90,7 +93,7 @@
 sta_manager.set_model_field(solver_obj.fields.elev_2d)
 
 # Define the scaling for the cost function so that J ~ O(1)
-J_scalar = Constant(solver_obj.dt / options.simulation_end_time)
+J_scalar = Constant(solver_obj.dt / options.simulation_end_time, domain=mesh2d)
 
 # Create inversion manager and add controls
 inv_manager = inversion_tools.InversionManager(

examples/columbia_plume/bathymetry.py

@@ -150,6 +150,6 @@ def get_boundary_relaxation_field(mask_func, bnd_id, dist_scale,
     # remove negative values
     mask_func.dat.data[mask_func.dat.data < 0.0] = 0.0
     if scalar is not None:
-        mask_func.assign(mask_func*scalar)
+        mask_func.interpolate(mask_func*scalar)
 
     return mask_func

examples/columbia_plume/cre-plume.py

@@ -221,8 +221,8 @@
 
 def interp_ocean_bnd(time):
     oce_bnd_interp.set_fields(time)
-    uvel_bnd_3d.assign(uvel_bnd_3d*ncom_vel_mask_3d)
-    vvel_bnd_3d.assign(vvel_bnd_3d*ncom_vel_mask_3d)
+    uvel_bnd_3d.interpolate(uvel_bnd_3d*ncom_vel_mask_3d)
+    vvel_bnd_3d.interpolate(vvel_bnd_3d*ncom_vel_mask_3d)
 
 
 # tides