Merge pull request #279 from thetisproject/recovery

HessianRecoverer2D

test/operations/test_diagnostics.py

@@ -0,0 +1,78 @@
+from thetis import *
+from thetis.diagnostics import *
+import pytest
+
+
+@pytest.fixture(params=[True, False])
+def interp(request):
+    return request.param
+
+
+def test_hessian_recovery2d(interp, tol=1.0e-08):
+    r"""
+    Apply Hessian recovery techniques to the quadratic
+    polynomial :math:`f(x, y) = \frac12(x^2 + y^2)` and
+    check that the result is the identity matrix.
+
+    We get the gradient approximation as part of the
+    recovery procedure, so check that this is also as
+    expected.
+
+    :arg interp: should the polynomial be interpolated as
+        as a :class:`Function`, or left as a UFL expression?
+    :kwarg tol: relative tolerance for value checking
+    """
+    mesh2d = UnitSquareMesh(4, 4)
+    x, y = SpatialCoordinate(mesh2d)
+    bowl = 0.5*(x**2 + y**2)
+    if interp:
+        P2_2d = get_functionspace(mesh2d, 'CG', 2)
+        bowl = interpolate(bowl, P2_2d)
+    P1v_2d = get_functionspace(mesh2d, 'CG', 1, vector=True)
+    P1t_2d = get_functionspace(mesh2d, 'CG', 1, tensor=True)
+
+    # Recover derivatives
+    g = Function(P1v_2d, name='Gradient')
+    H = Function(P1t_2d, name='Hessian')
+    hessian_recoverer = HessianRecoverer2D(bowl, H, g)
+    hessian_recoverer.solve()
+
+    # Check values
+    g_expect = Function(P1v_2d, name='Expected gradient')
+    g_expect.interpolate(as_vector([x, y]))
+    H_expect = Function(P1t_2d, name='Expected Hessian')
+    H_expect.interpolate(Identity(2))
+    err_g = errornorm(g, g_expect)/norm(g_expect)
+    assert err_g < tol, f'Gradient approximation error {err_g:.4e}'
+    err_H = errornorm(H, H_expect)/norm(H_expect)
+    assert err_H < tol, f'Hessian approximation error {err_H:.4e}'
+
+
+def test_vorticity_calculation2d(interp, tol=1.0e-08):
+    r"""
+    Calculate the vorticity of the velocity field
+    :math:`\mathbf u(x, y) = \frac12(y, -x)` and check
+    that the result is unity.
+
+    :arg interp: should the velocity be interpolated as
+        as a :class:`Function`, or left as a UFL expression?
+    :kwarg tol: relative tolerance for value checking
+    """
+    mesh2d = UnitSquareMesh(4, 4)
+    x, y = SpatialCoordinate(mesh2d)
+    uv = 0.5*as_vector([y, -x])
+    if interp:
+        P1v_2d = get_functionspace(mesh2d, 'CG', 1, vector=True)
+        uv = interpolate(uv, P1v_2d)
+    P1_2d = get_functionspace(mesh2d, 'CG', 1)
+
+    # Recover vorticity
+    omega = Function(P1_2d, name='Vorticity')
+    vorticity_calculator = VorticityCalculator2D(uv, omega)
+    vorticity_calculator.solve()
+
+    # Check values
+    omega_expect = Function(P1_2d, name='Expected vorticity')
+    omega_expect.assign(1.0)
+    err = errornorm(omega, omega_expect)/norm(omega_expect)
+    assert err < tol, f'Vorticity approximation error {err:.4e}'

test/swe2d/test_anisotropic.py

@@ -57,7 +57,7 @@ def run(**model_options):
         'pc_factor_mat_solver_type': 'mumps',
     }
     options.output_directory = 'outputs'
-    options.fields_to_export = ['uv_2d', 'elev_2d']
+    options.fields_to_export = ['uv_2d', 'elev_2d', 'hessian_2d']
     options.use_grad_div_viscosity_term = False
     options.horizontal_viscosity = viscosity
     options.quadratic_drag_coefficient = drag_coefficient
@@ -68,6 +68,16 @@ def run(**model_options):
     options.update(model_options)
     solver_obj.create_equations()
 
+    # recover Hessian
+    if 'hessian_2d' in field_metadata:
+        field_metadata.pop('hessian_2d')
+    P1t_2d = get_functionspace(mesh2d, 'CG', 1, tensor=True)
+    hessian_2d = Function(P1t_2d)
+    u_2d = solver_obj.fields.uv_2d[0]
+    hessian_calculator = thetis.diagnostics.HessianRecoverer2D(u_2d, hessian_2d)
+    solver_obj.add_new_field(hessian_2d, 'hessian_2d', 'Hessian of x-velocity', 'Hessian2d',
+                             preproc_func=hessian_calculator.solve)
+
     # apply boundary conditions
     solver_obj.bnd_functions['shallow_water'] = {
         1: {'uv': inflow_velocity},  # inflow condition used upstream
@@ -143,4 +153,4 @@ def test_sipg(family):
 
 
 if __name__ == '__main__':
-    print_output(run(no_exports=False, element_family='rt-dg'))
+    print_output(f"Converged in {run(no_exports=False, element_family='rt-dg')} iterations")

test/swe2d/test_rossby_wave.py

@@ -152,7 +152,7 @@ def run(refinement_level, **model_options):
     # Physics
     g = physical_constants['g_grav'].values()[0]
     physical_constants['g_grav'].assign(1.0)
-    P1_2d = FunctionSpace(mesh2d, "CG", 1)
+    P1_2d = get_functionspace(mesh2d, "CG", 1)
     bathymetry2d = Function(P1_2d).assign(1.0)
 
     # Create solver object
@@ -163,15 +163,26 @@ def run(refinement_level, **model_options):
         options.swe_timestepper_options.use_automatic_timestep = False
     options.timestep = 0.96/refinement_level if stepper == 'SSPRK33' else 9.6/refinement_level
     options.simulation_export_time = 5.0
-    options.simulation_end_time = model_options.get('simulation_end_time')
+    options.simulation_end_time = model_options.get('simulation_end_time', 120.0)
     options.use_grad_div_viscosity_term = False
     options.use_grad_depth_viscosity_term = False
     options.horizontal_viscosity = None
     solver_obj.create_function_spaces()
     options.coriolis_frequency = interpolate(y, solver_obj.function_spaces.P1_2d)
     options.swe_timestepper_options.solver_parameters['ksp_rtol'] = 1.0e-04
     options.no_exports = True
+    options.fields_to_export = ['uv_2d', 'elev_2d', 'vorticity_2d']
     options.update(model_options)
+    solver_obj.create_equations()
+
+    # Calculate vorticity
+    if 'vorticity_2d' in field_metadata:
+        field_metadata.pop('vorticity_2d')
+    vorticity_2d = Function(P1_2d)
+    uv_2d = solver_obj.fields.uv_2d
+    vorticity_calculator = thetis.diagnostics.VorticityCalculator2D(uv_2d, vorticity_2d)
+    solver_obj.add_new_field(vorticity_2d, 'vorticity_2d', 'Fluid vorticity', 'Vorticity2d',
+                             preproc_func=vorticity_calculator.solve)
 
     # Apply boundary conditions
     for tag in mesh2d.exterior_facets.unique_markers:
@@ -267,4 +278,5 @@ def test_convergence(stepper, family):
 # ---------------------------
 
 if __name__ == '__main__':
-    test_convergence('DIRK22', 'bdm-dg')
+    run(24, swe_timestepper_type='DIRK22', element_family='bdm-dg',
+        expansion_order=1, no_exports=False)

test/swe2d/test_standing_wave.py

@@ -114,4 +114,4 @@ def test_standing_wave_channel(timesteps, max_rel_err, timestepper, tmpdir, do_e
 
 
 if __name__ == '__main__':
-    test_standing_wave_channel(do_export=True)
+    test_standing_wave_channel(10, 1.6e-02, 'CrankNicolson', 'outputs', do_export=True)

thetis/__init__.py

@@ -14,6 +14,7 @@
 import thetis.timezone as timezone      # NOQA
 import thetis.turbines  # NOQA
 import thetis.optimisation  # NOQA
+import thetis.diagnostics  # NOQA
 from thetis._version import get_versions
 from thetis.assembledschur import AssembledSchurPC  # NOQA
 from thetis.options import TidalTurbineFarmOptions  # NOQA

thetis/diagnostics.py

@@ -0,0 +1,149 @@
+"""
+Classes for computing diagnostics.
+"""
+from .utility import *
+
+
+__all__ = ["VorticityCalculator2D", "HessianRecoverer2D"]
+
+
+class VorticityCalculator2D(object):
+    r"""
+    Linear solver for recovering fluid vorticity.
+
+    It is recommended that the vorticity is sought
+    in :math:`\mathbb P1` space.
+
+    :arg uv_2d: horizontal velocity :class:`Function`.
+    :arg vorticity_2d: :class:`Function` to hold calculated vorticity.
+    :kwargs: to be passed to the :class:`LinearVariationalSolver`.
+    """
+    @PETSc.Log.EventDecorator("thetis.VorticityCalculator2D.__init__")
+    def __init__(self, uv_2d, vorticity_2d, **kwargs):
+        fs = vorticity_2d.function_space()
+        dim = fs.mesh().topological_dimension()
+        if dim != 2:
+            raise ValueError(f'Dimension {dim} not supported')
+        if element_continuity(fs.ufl_element()).horizontal != 'cg':
+            raise ValueError('Vorticity must be calculated in a continuous space')
+        n = FacetNormal(fs.mesh())
+
+        # Weak formulation
+        test = TestFunction(fs)
+        a = TrialFunction(fs)*test*dx
+        L = -inner(perp(uv_2d), grad(test))*dx \
+            + dot(perp(uv_2d), test*n)*ds \
+            + dot(avg(perp(uv_2d)), jump(test, n))*dS
+
+        # Setup vorticity solver
+        prob = LinearVariationalProblem(a, L, vorticity_2d)
+        kwargs.setdefault('solver_parameters', {
+            "ksp_type": "cg",
+            "pc_type": "bjacobi",
+            "sub_pc_type": "ilu",
+        })
+        self.solver = LinearVariationalSolver(prob, **kwargs)
+
+    @PETSc.Log.EventDecorator("thetis.VorticityCalculator2D.solve")
+    def solve(self):
+        self.solver.solve()
+
+
+class HessianRecoverer2D(object):
+    r"""
+    Linear solver for recovering Hessians.
+
+    Hessians are recoved using double :math:`L^2`
+    projection, which is implemented using a
+    mixed finite element method.
+
+    It is recommended that gradients and Hessians
+    are sought in :math:`\mathbb P1` space of
+    appropriate dimension.
+
+    :arg field_2d: :class:`Function` to recover the Hessian of.
+    :arg hessian_2d: :class:`Function` to hold recovered Hessian.
+    :kwarg gradient_2d: :class:`Function` to hold recovered gradient.
+    :kwargs: to be passed to the :class:`LinearVariationalSolver`.
+    """
+    @PETSc.Log.EventDecorator("thetis.HessianRecoverer2D.__init__")
+    def __init__(self, field_2d, hessian_2d, gradient_2d=None, **kwargs):
+        self.hessian_2d = hessian_2d
+        self.gradient_2d = gradient_2d
+        Sigma = hessian_2d.function_space()
+        mesh = Sigma.mesh()
+        dim = mesh.topological_dimension()
+        if dim != 2:
+            raise ValueError(f'Dimension {dim} not supported')
+        n = FacetNormal(mesh)
+
+        # Extract function spaces
+        if element_continuity(Sigma.ufl_element()).horizontal != 'cg':
+            raise ValueError('Hessians must be calculated in a continuous space')
+        if Sigma.dof_dset.dim != (2, 2):
+            raise ValueError('Expecting a square tensor function')
+        if gradient_2d is None:
+            V = get_functionspace(mesh, 'CG', 1, vector=True)
+        else:
+            V = gradient_2d.function_space()
+            if element_continuity(V.ufl_element()).horizontal != 'cg':
+                raise ValueError('Gradients must be calculated in a continuous space')
+            if V.dof_dset.dim != (2,):
+                raise ValueError('Expecting a 2D vector function')
+
+        # Setup function spaces
+        W = V*Sigma
+        g, H = TrialFunctions(W)
+        phi, tau = TestFunctions(W)
+        sol = Function(W)
+        self._gradient, self._hessian = sol.split()
+
+        # The formulation is chosen such that f does not need to have any
+        # finite element derivatives
+        a = inner(tau, H)*dx \
+            + inner(div(tau), g)*dx \
+            + inner(phi, g)*dx \
+            - dot(g, dot(tau, n))*ds \
+            - dot(avg(g), jump(tau, n))*dS
+        L = field_2d*dot(phi, n)*ds \
+            + avg(field_2d)*jump(phi, n)*dS \
+            - field_2d*div(phi)*dx
+
+        # Apply stationary preconditioners in the Schur complement to get away
+        # with applying GMRES to the whole mixed system
+        sp = {
+            "mat_type": "aij",
+            "ksp_type": "gmres",
+            "ksp_max_it": 20,
+            "pc_type": "fieldsplit",
+            "pc_fieldsplit_type": "schur",
+            "pc_fieldsplit_0_fields": "1",
+            "pc_fieldsplit_1_fields": "0",
+            "pc_fieldsplit_schur_precondition": "selfp",
+            "fieldsplit_0_ksp_type": "preonly",
+            "fieldsplit_1_ksp_type": "preonly",
+            "fieldsplit_1_pc_type": "gamg",
+            "fieldsplit_1_mg_levels_ksp_max_it": 5,
+        }
+        if COMM_WORLD.size == 1:
+            sp["fieldsplit_0_pc_type"] = "ilu"
+            sp["fieldsplit_1_mg_levels_pc_type"] = "ilu"
+        else:
+            sp["fieldsplit_0_pc_type"] = "bjacobi"
+            sp["fieldsplit_0_sub_ksp_type"] = "preonly"
+            sp["fieldsplit_0_sub_pc_type"] = "ilu"
+            sp["fieldsplit_1_mg_levels_pc_type"] = "bjacobi"
+            sp["fieldsplit_1_mg_levels_sub_ksp_type"] = "preonly"
+            sp["fieldsplit_1_mg_levels_sub_pc_type"] = "ilu"
+
+        # Setup solver
+        prob = LinearVariationalProblem(a, L, sol)
+        kwargs.setdefault('solver_parameters', sp)
+        self.solver = LinearVariationalSolver(prob, **kwargs)
+
+    @PETSc.Log.EventDecorator("thetis.HessianRecoverer2D.solve")
+    def solve(self):
+        self.solver.solve()
+        self.hessian_2d.assign(self._hessian)
+        if self.gradient_2d is not None:
+            self.gradient_2d.assign(self._gradient)

thetis/exporter.py

@@ -21,13 +21,15 @@ def get_visu_space(fs):
     :arg fs: function space
     :return: function space for VTK visualization
     """
-    is_vector = len(fs.ufl_element().value_shape()) == 1
     mesh = fs.mesh()
     family = 'Lagrange' if is_cg(fs) else 'Discontinuous Lagrange'
-    if is_vector:
+    if len(fs.ufl_element().value_shape()) == 1:
         dim = fs.ufl_element().value_shape()[0]
         visu_fs = get_functionspace(mesh, family, 1, family, 1,
                                     vector=True, dim=dim)
+    elif len(fs.ufl_element().value_shape()) == 2:
+        visu_fs = get_functionspace(mesh, family, 1, family, 1,
+                                    tensor=True)
     else:
         visu_fs = get_functionspace(mesh, family, 1, family, 1)
     # make sure that you always get the same temp work function