Discretisation changes for RT-DG pair:

test/swe2d/supported_elements.py

@@ -0,0 +1,9 @@
+import pytest
+
+cases = [('rt-dg', 0), ('rt-dg', 1), ('dg-dg', 1), ('dg-cg', 1), ('bdm-dg', 1)]
+case_ids = [f'{fam}{deg}' for fam, deg in cases]
+
+
+@pytest.fixture(params=cases, ids=case_ids)
+def element_family_and_degree(request):
+    return request.param

test/swe2d/test_anisotropic.py

@@ -16,10 +16,10 @@
     pp.1053--1063 (2020), DOI: 10.1007/s42452-020-2745-9, URL: https://rdcu.be/b35wZ.
 """
 from thetis import *
+from supported_elements import *
 import thetis.diagnostics as diagnostics
 from petsc4py import PETSc
 import numpy as np
-import pytest
 import os
 
 
@@ -215,25 +215,22 @@ def estimate_error(mesh, **model_options):
 # standard tests for pytest
 # ---------------------------
 
-@pytest.fixture(params=['dg-cg', 'dg-dg', 'rt-dg', 'bdm-dg'])
-def family(request):
-    return request.param
-
-
-def test_sipg(family):
-    snes_it = run(element_family=family)
+def test_sipg(element_family_and_degree):
+    family, degree = element_family_and_degree
+    snes_it = run(element_family=family, polynomial_degree=degree)
     expected = 3
     msg = f'snes iterations exceed expected: {snes_it} > {expected}'
     assert snes_it <= expected, msg
 
 
-def test_dwr(family):
+def test_dwr(element_family_and_degree):
     n = 5
     mesh = RectangleMesh(12 * n, 5 * n, 1200, 500)
+    family, degree = element_family_and_degree
     # NOTE: Building a MeshHierarchy on a mesh loaded
     #       from a DMPlex stored as HDF5 appears to
     #       be broken.
-    estimate_error(mesh, element_family=family, no_exports=True)
+    estimate_error(mesh, element_family=family, polynomial_degree=degree, no_exports=True)
 
 # ---------------------------
 # run individual setup for debugging

test/swe2d/test_atmospheric_pressure.py

@@ -12,14 +12,15 @@
 #
 from thetis import *
 import pytest
+import supported_elements
 
 
-@pytest.mark.parametrize("element_family", [
-    'dg-dg', 'rt-dg', 'dg-cg', 'bdm-dg'])
+@pytest.mark.parametrize("element_family_and_degree",
+                         supported_elements.cases)
 @pytest.mark.parametrize("timestepper", [
     'CrankNicolson', 'SSPRK33', ])
-def test_pressure_forcing(element_family, timestepper):
-    order = 1
+def test_pressure_forcing(element_family_and_degree, timestepper):
+    element_family, order = element_family_and_degree
 
     lx = 10000
     ly = 10000

test/swe2d/test_rossby_wave.py

@@ -15,6 +15,7 @@
 (2008), Journal of Geophysical Research: Oceans, 113(C7).
 """
 from thetis import *
+from supported_elements import *
 import pytest
 
 
@@ -267,14 +268,10 @@ def stepper(request):
     return request.param
 
 
-@pytest.fixture(params=['dg-dg', 'dg-cg', 'rt-dg', 'bdm-dg'])
-def family(request):
-    return request.param
-
-
-def test_convergence(stepper, family):
+def test_convergence(stepper, element_family_and_degree):
+    family, degree = element_family_and_degree
     run_convergence([24, 48], swe_timestepper_type=stepper,
-                    simulation_end_time=30.0, polynomial_degree=1, element_family=family,
+                    simulation_end_time=30.0, polynomial_degree=degree, element_family=family,
                     no_exports=True, expansion_order=1)
 
 

test/swe2d/test_steady_state_basin_mms.py

@@ -10,6 +10,7 @@
 from thetis import *
 from scipy import stats
 import pytest
+from supported_elements import *
 
 
 def setup7(x, lx, ly, h0, f0, nu0, g):
@@ -177,7 +178,7 @@ def run(setup, refinement, order, do_export=True, options=None,
         coriolis_func.project(sdict['cori_expr'])
         solver_obj.options.coriolis_frequency = coriolis_func
     if 'visc_expr' in sdict:
-        viscosity_space = FunctionSpace(solver_obj.mesh2d, "CG", order)
+        viscosity_space = FunctionSpace(solver_obj.mesh2d, "CG", max(order, 1))
         viscosity_func = Function(viscosity_space, name='viscosity')
         viscosity_func.project(sdict['visc_expr'])
         solver_obj.options.horizontal_viscosity = viscosity_func
@@ -325,22 +326,24 @@ def setup(request):
     return request.param
 
 
-@pytest.fixture(params=['rt-dg', 'dg-dg', 'dg-cg', 'bdm-dg'])
-def element_family(request):
-    return request.param
-
-
 @pytest.fixture(params=['CrankNicolson'])
 def timestepper_type(request):
     return request.param
 
 
-def test_steady_state_basin_convergence(setup, element_family, timestepper_type):
+def test_steady_state_basin_convergence(setup, element_family_and_degree, timestepper_type):
     sp = {'ksp_type': 'preonly', 'pc_type': 'lu', 'snes_monitor': None,
           'mat_type': 'aij'}
+    if setup == setup9 and element_family_and_degree == ('rt-dg', 0):
+        # SIPG is invalid for lowest order RT:
+        # because the basis functions do not support shear within the element,
+        # first two terms with either grad(trial) or grad(test) filter out
+        # the shear component
+        pytest.xfail("No valid viscosity scheme for RT-DG0 available")
+    element_family, degree = element_family_and_degree
     options = {
         'element_family': element_family,
         'swe_timestepper_type': timestepper_type
     }
-    run_convergence(setup, [1, 2, 4, 6], 1, options=options,
+    run_convergence(setup, [1, 2, 4, 6], degree, options=options,
                     solver_parameters=sp, save_plot=False)

test/swe2d/test_steady_state_channel_mms.py

@@ -1,18 +1,13 @@
 from thetis import *
+from supported_elements import *
 import math
-import pytest
 
 
-@pytest.fixture(params=['rt-dg', 'dg-dg', 'dg-cg', 'bdm-dg'])
-def element_family(request):
-    return request.param
-
-
-def test_steady_state_channel_mms(element_family, do_exports=False):
+def test_steady_state_channel_mms(element_family_and_degree, do_exports=False):
     lx = 5e3
     ly = 1e3
 
-    order = 1
+    element_family, order = element_family_and_degree
     # minimum resolution
     min_cells = 48
     n = 1  # number of timesteps
@@ -24,14 +19,17 @@ def test_steady_state_channel_mms(element_family, do_exports=False):
     k = 4.0*math.pi/lx
     Q = H0*1.0  # flux (depth-integrated velocity)
     eta0 = 1.0  # free surface amplitude
-    C_D = 0.0025  # quadratic drag coefficient
+    C_D = 0.01  # quadratic drag coefficient
 
     xhat = Identity(2)[0, :]
 
     eta_errs = []
     u_errs = []
     for i in range(4):
-        mesh2d = RectangleMesh(min_cells*2**i, 1, lx, ly)
+        # for rt/bdm accuracy depends on triangle quality;
+        # avoid elements getting flatter by making channel narrower
+        # (wouldn't be an issue if we'd refine in both x and y)
+        mesh2d = RectangleMesh(min_cells*2**i, 1, lx, ly/2**i)
         x = mesh2d.coordinates
         eta_expr = eta0*cos(k*x[0])
         H = H0+eta0*cos(k*x[0])

thetis/shallowwater_eq.py

@@ -359,11 +359,12 @@ def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions)
 
         if grad_eta_by_parts:
             f = -g_grav*head*nabla_div(self.u_test)*self.dx
-            if uv is not None:
-                head_star = avg(head) + sqrt(avg(total_h)/g_grav)*jump(uv, self.normal)
-            else:
-                head_star = avg(head)
-            f += g_grav*head_star*jump(self.u_test, self.normal)*self.dS
+            if self.u_continuity == 'dg':
+                if uv is not None:
+                    head_star = avg(head) + sqrt(avg(total_h)/g_grav)*jump(uv, self.normal)
+                else:
+                    head_star = avg(head)
+                f += g_grav*head_star*jump(self.u_test, self.normal)*self.dS
             for bnd_marker in self.boundary_markers:
                 funcs = bnd_conditions.get(bnd_marker)
                 ds_bnd = ds(int(bnd_marker), degree=self.quad_degree)
@@ -421,10 +422,14 @@ def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions)
         if hu_by_parts:
             f = -inner(grad(self.eta_test), total_h*uv)*self.dx
             if self.eta_is_dg:
-                h = avg(total_h)
-                uv_rie = avg(uv) + sqrt(g_grav/h)*jump(eta, self.normal)
-                hu_star = h*uv_rie
-                f += inner(jump(self.eta_test, self.normal), hu_star)*self.dS
+                if self.u_continuity == 'dg':
+                    h = avg(total_h)
+                    uv_rie = avg(uv) + sqrt(g_grav/h)*jump(eta, self.normal)
+                    hu_star = h*uv_rie
+                    f += inner(jump(self.eta_test, self.normal), hu_star)*self.dS
+                else:
+                    f += inner(jump(self.eta_test, self.normal), avg(total_h*uv))*self.dS
+
             for bnd_marker in self.boundary_markers:
                 funcs = bnd_conditions.get(bnd_marker)
                 ds_bnd = ds(int(bnd_marker), degree=self.quad_degree)