From 8e2889ad10a26b12423962ac42ca0c31dd444249 Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Fri, 26 Nov 2021 15:09:09 +0000 Subject: [PATCH 01/10] Extend get_functionspace to tensor case --- thetis/utility.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/thetis/utility.py b/thetis/utility.py index 59b3741b3..090629ccc 100755 --- a/thetis/utility.py +++ b/thetis/utility.py @@ -114,7 +114,7 @@ def __setattr__(self, key, value): def get_functionspace(mesh, h_family, h_degree, v_family=None, v_degree=None, - vector=False, hdiv=False, variant=None, v_variant=None, + vector=False, tensor=False, hdiv=False, variant=None, v_variant=None, **kwargs): cell_dim = mesh.cell_dimension() assert cell_dim in [2, (2, 1)], 'Unsupported cell dimension' @@ -146,7 +146,8 @@ def get_functionspace(mesh, h_family, h_degree, v_family=None, v_degree=None, else: elt = FiniteElement(h_family, mesh.ufl_cell(), h_degree, variant=variant) - constructor = VectorFunctionSpace if vector else FunctionSpace + assert not (vector and tensor) + constructor = TensorFunctionSpace if tensor else VectorFunctionSpace if vector else FunctionSpace return constructor(mesh, elt, **kwargs) From 615af9611e0952bc7d7228d3656ab4dd21a8f484 Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Fri, 26 Nov 2021 15:09:27 +0000 Subject: [PATCH 02/10] Extend get_visu_space to tensor case --- thetis/exporter.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/thetis/exporter.py b/thetis/exporter.py index 6831575e6..2eb4f079f 100644 --- a/thetis/exporter.py +++ b/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 From 5a659e96c0e559579f16d52b81d15d5c803dd273 Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Fri, 26 Nov 2021 14:28:26 +0000 Subject: [PATCH 03/10] Rework VorticityCalculator2D in a new file --- thetis/diagnostics.py | 46 +++++++++++++++++++++++++++++++++++++++++++ thetis/utility.py | 42 --------------------------------------- 2 files changed, 46 insertions(+), 42 deletions(-) create mode 100644 thetis/diagnostics.py diff --git a/thetis/diagnostics.py b/thetis/diagnostics.py new file mode 100644 index 000000000..ba81a47bc --- /dev/null +++ b/thetis/diagnostics.py @@ -0,0 +1,46 @@ +""" +Classes for computing diagnostics. +""" +from .utility import * + + +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() diff --git a/thetis/utility.py b/thetis/utility.py index 090629ccc..df56f9d5f 100755 --- a/thetis/utility.py +++ b/thetis/utility.py @@ -893,48 +893,6 @@ def get_total_depth(self, eta): return total_h -class VorticityCalculator2D(object): - r""" - Linear solver for recovering fluid vorticity. - - It is recommended that the vorticity is sought - in :math:`\mathbb P1` space. - - :arg vorticity_2d: :class:`Function` to hold calculated vorticity. - :arg solver_obj: :class:`FlowSolver2d` object. - :kwargs: to be passed to the :class:`LinearVariationalSolver`. - """ - @PETSc.Log.EventDecorator("thetis.VorticityCalculator2D.__init__") - def __init__(self, vorticity_2d, solver_obj, **kwargs): - fs = vorticity_2d.function_space() - if element_continuity(fs.ufl_element()).horizontal != 'cg': - raise NotImplementedError - n = FacetNormal(fs.mesh()) - - # Weak formulation - test = TestFunction(fs) - uv, elev = split(solver_obj.fields.solution_2d) - a = TrialFunction(fs)*test*dx - L = -inner(perp(uv), grad(test))*dx \ - + dot(perp(uv), n)*test*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) - - if 'vorticity_2d' not in solver_obj.options.fields_to_export: - print_output("NOTE: 'vorticity_2d' will be calculated but not exported.") - - @PETSc.Log.EventDecorator("thetis.VorticityCalculator2D.solve") - def solve(self): - self.solver.solve() - - class DepthIntegratedPoissonSolver(object): r""" Construct solvers for Poisson equation and updating velocities From 46a0dfe60faad29b99823e25d3cfccbc57e4cb04 Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Fri, 26 Nov 2021 14:43:40 +0000 Subject: [PATCH 04/10] HessianRecoverer2D --- thetis/diagnostics.py | 100 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 100 insertions(+) diff --git a/thetis/diagnostics.py b/thetis/diagnostics.py index ba81a47bc..591eac69e 100644 --- a/thetis/diagnostics.py +++ b/thetis/diagnostics.py @@ -44,3 +44,103 @@ def __init__(self, uv_2d, vorticity_2d, **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) From aa966cc1bb9621232cc46d4b4f64dbaffe064b4d Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Wed, 1 Dec 2021 09:20:06 +0000 Subject: [PATCH 05/10] Import diagnostics --- thetis/__init__.py | 1 + thetis/diagnostics.py | 3 +++ 2 files changed, 4 insertions(+) diff --git a/thetis/__init__.py b/thetis/__init__.py index 2b24a8af9..fd8da6d7d 100644 --- a/thetis/__init__.py +++ b/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 diff --git a/thetis/diagnostics.py b/thetis/diagnostics.py index 591eac69e..6653d43d7 100644 --- a/thetis/diagnostics.py +++ b/thetis/diagnostics.py @@ -4,6 +4,9 @@ from .utility import * +__all__ = ["VorticityCalculator2D", "HessianRecoverer2D"] + + class VorticityCalculator2D(object): r""" Linear solver for recovering fluid vorticity. From 9d8d65b9983cdedce29fbb6b7f16e057a97e000d Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Wed, 1 Dec 2021 09:20:39 +0000 Subject: [PATCH 06/10] Valid debug for test_standing_wave --- test/swe2d/test_standing_wave.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/swe2d/test_standing_wave.py b/test/swe2d/test_standing_wave.py index d4e12adb6..79b05ca68 100644 --- a/test/swe2d/test_standing_wave.py +++ b/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) From d9a15077ecdd447885de17504ed289b9937e4878 Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Wed, 1 Dec 2021 09:32:26 +0000 Subject: [PATCH 07/10] Test vorticity calculation --- test/swe2d/test_rossby_wave.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/test/swe2d/test_rossby_wave.py b/test/swe2d/test_rossby_wave.py index 3910c5499..42b279c32 100644 --- a/test/swe2d/test_rossby_wave.py +++ b/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,7 +163,7 @@ 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 @@ -171,7 +171,18 @@ def run(refinement_level, **model_options): 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) From bd82e5fadc1fa388a9089fa730eee868bd039668 Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Wed, 1 Dec 2021 09:41:24 +0000 Subject: [PATCH 08/10] Test Hessian recovery --- test/swe2d/test_anisotropic.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/test/swe2d/test_anisotropic.py b/test/swe2d/test_anisotropic.py index fe3c104e6..6f5bcfd0e 100644 --- a/test/swe2d/test_anisotropic.py +++ b/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") From 3303f8509e4786dd90aff55ea78fdba613290b61 Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Thu, 2 Dec 2021 08:26:02 +0000 Subject: [PATCH 09/10] Test Hessian recovery --- test/operations/test_diagnostics.py | 48 +++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 test/operations/test_diagnostics.py diff --git a/test/operations/test_diagnostics.py b/test/operations/test_diagnostics.py new file mode 100644 index 000000000..111197e3c --- /dev/null +++ b/test/operations/test_diagnostics.py @@ -0,0 +1,48 @@ +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}' From 9bc0815c5fc55cf83201bc410070b2c0acb832b2 Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Thu, 2 Dec 2021 08:36:31 +0000 Subject: [PATCH 10/10] Test vorticity calculation --- test/operations/test_diagnostics.py | 30 +++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/test/operations/test_diagnostics.py b/test/operations/test_diagnostics.py index 111197e3c..cbe4701cd 100644 --- a/test/operations/test_diagnostics.py +++ b/test/operations/test_diagnostics.py @@ -46,3 +46,33 @@ def test_hessian_recovery2d(interp, tol=1.0e-08): 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}'