From f3b77c07dae845694c9cb87aca8081a54e1bb13b Mon Sep 17 00:00:00 2001 From: Witt-D Date: Wed, 27 Mar 2024 12:28:39 +0000 Subject: [PATCH 01/17] addressed docstring and other minor issues --- gusto/diagnostics.py | 4 +++- gusto/equations.py | 1 + 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/gusto/diagnostics.py b/gusto/diagnostics.py index e7ffd62c2..efa353336 100644 --- a/gusto/diagnostics.py +++ b/gusto/diagnostics.py @@ -1726,7 +1726,7 @@ def setup(self, domain, state_fields, vorticity_type=None): n = FacetNormal(domain.mesh) a = inner(vort, w) * dx L = inner(u, curl(w)) * dx - jump(cross(w, u), n) * dS_h - if vorticity_type != 'relative': + if vorticity_type == 'absolute': Omega = as_vector((0, 0, self.parameters.Omega)) L += inner(2*Omega, w) * dx @@ -1769,6 +1769,8 @@ class CompressibleAbsoluteVorticity(CompressibleVorticity): def __init__(self, parameters, space=None, method='solve'): u""" Args: + parameters (:class:`CompressibleEulerParameters`): the configuration + object containing the physical parameters for this equation. space (:class:`FunctionSpace`, optional): the function space to evaluate the diagnostic field in. Defaults to None, in which case a default space will be chosen for this diagnostic. diff --git a/gusto/equations.py b/gusto/equations.py index 5d6197914..d24b19eb7 100644 --- a/gusto/equations.py +++ b/gusto/equations.py @@ -1074,6 +1074,7 @@ def __init__(self, domain, parameters, sponge=None, # -------------------------------------------------------------------- # if parameters.Omega is not None: + # TODO add linerisation and label for this Omega = as_vector((0, 0, parameters.Omega)) coriolis_form = coriolis(subject(prognostic( inner(w, cross(2*Omega, u))*dx, "u"), self.X)) From 2a9da25689e1e86bf6a23978c6c789489f0aed10 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Wed, 27 Mar 2024 15:21:51 +0000 Subject: [PATCH 02/17] renamed vorticities --- gusto/diagnostics.py | 49 +++++++++++++++++++++++--------------------- 1 file changed, 26 insertions(+), 23 deletions(-) diff --git a/gusto/diagnostics.py b/gusto/diagnostics.py index efa353336..4319b6420 100644 --- a/gusto/diagnostics.py +++ b/gusto/diagnostics.py @@ -26,7 +26,8 @@ "Perturbation", "Theta_e", "InternalEnergy", "PotentialEnergy", "ThermodynamicKineticEnergy", "Dewpoint", "Temperature", "Theta_d", "RelativeHumidity", "Pressure", "Exner_Vt", "HydrostaticImbalance", "Precipitation", - "PotentialVorticity", "RelativeVorticity", "AbsoluteVorticity", "Divergence", + "ShallowWaterPotentialVorticity", "ShallowWaterRelativeVorticity", "ShallowWaterAbsoluteVorticity", + "Vorticity", "RelativeVorticity", "AbsoluteVorticity", "Divergence", "BruntVaisalaFrequencySquared", "TracerDensity"] @@ -796,7 +797,7 @@ def setup(self, domain, state_fields): class ShallowWaterPotentialEnstrophy(DiagnosticField): """Diagnostic (dry) compressible kinetic energy density.""" - def __init__(self, base_field_name="PotentialVorticity", space=None, + def __init__(self, base_field_name="ShallowWaterPotentialVorticity", space=None, method='interpolate'): """ Args: @@ -809,15 +810,17 @@ def __init__(self, base_field_name="PotentialVorticity", space=None, for this diagnostic. Valid options are 'interpolate', 'project', 'assign' and 'solve'. Defaults to 'interpolate'. """ - base_enstrophy_names = ["PotentialVorticity", "RelativeVorticity", "AbsoluteVorticity"] + base_enstrophy_names = ["ShallowWaterPotentialVorticity", + "ShallowWaterRelativeVorticity", + "ShallowWaterAbsoluteVorticity"] if base_field_name not in base_enstrophy_names: raise ValueError( f"Don't know how to compute enstrophy with base_field_name={base_field_name};" + f"base_field_name should be one of {base_enstrophy_names}") # Work out required fields - if base_field_name in ["PotentialVorticity", "AbsoluteVorticity"]: + if base_field_name in ["ShallowWaterPotentialVorticity", "ShallowWaterAbsoluteVorticity"]: required_fields = (base_field_name, "D") - elif base_field_name == "RelativeVorticity": + elif base_field_name == "ShallowWaterRelativeVorticity": required_fields = (base_field_name, "D", "coriolis") else: raise NotImplementedError(f'Enstrophy with vorticity {base_field_name} not implemented') @@ -839,17 +842,17 @@ def setup(self, domain, state_fields): domain (:class:`Domain`): the model's domain object. state_fields (:class:`StateFields`): the model's field container. """ - if self.base_field_name == "PotentialVorticity": - pv = state_fields("PotentialVorticity") + if self.base_field_name == "ShallowWaterPotentialVorticity": + pv = state_fields("ShallowWaterPotentialVorticity") D = state_fields("D") self.expr = 0.5*pv**2*D - elif self.base_field_name == "RelativeVorticity": - zeta = state_fields("RelativeVorticity") + elif self.base_field_name == "ShallowWaterRelativeVorticity": + zeta = state_fields("ShallowWaterRelativeVorticity") D = state_fields("D") f = state_fields("coriolis") self.expr = 0.5*(zeta + f)**2/D - elif self.base_field_name == "AbsoluteVorticity": - zeta_abs = state_fields("AbsoluteVorticity") + elif self.base_field_name == "ShallowWaterAbsoluteVorticity": + zeta_abs = state_fields("ShallowWaterAbsoluteVorticity") D = state_fields("D") self.expr = 0.5*(zeta_abs)**2/D else: @@ -1547,7 +1550,7 @@ def compute(self): self.field.assign(self.field + self.flux) -class Vorticity(DiagnosticField): +class ShallowWaterVorticity(DiagnosticField): """Base diagnostic field class for shallow-water vorticity variables.""" def setup(self, domain, state_fields, vorticity_type=None): @@ -1602,9 +1605,9 @@ def setup(self, domain, state_fields, vorticity_type=None): self.evaluator = LinearVariationalSolver(problem, solver_parameters={"ksp_type": "cg"}) -class PotentialVorticity(Vorticity): +class ShallowWaterPotentialVorticity(ShallowWaterVorticity): u"""Diagnostic field for shallow-water potential vorticity, q=(∇×(u+f))/D""" - name = "PotentialVorticity" + name = "ShallowWaterPotentialVorticity" def __init__(self, space=None, method='solve'): """ @@ -1631,9 +1634,9 @@ def setup(self, domain, state_fields): super().setup(domain, state_fields, vorticity_type="potential") -class AbsoluteVorticity(Vorticity): +class ShallowWaterAbsoluteVorticity(ShallowWaterVorticity): u"""Diagnostic field for absolute vorticity, ζ=∇×(u+f)""" - name = "AbsoluteVorticity" + name = "ShallowWaterAbsoluteVorticity" def __init__(self, space=None, method='solve'): """ @@ -1659,9 +1662,9 @@ def setup(self, domain, state_fields): super().setup(domain, state_fields, vorticity_type="absolute") -class RelativeVorticity(Vorticity): +class ShallowWaterRelativeVorticity(ShallowWaterVorticity): u"""Diagnostic field for relative vorticity, ζ=∇×u""" - name = "RelativeVorticity" + name = "ShallowWaterRelativeVorticity" def __init__(self, space=None, method='solve'): """ @@ -1687,7 +1690,7 @@ def setup(self, domain, state_fields): super().setup(domain, state_fields, vorticity_type="relative") -class CompressibleVorticity(DiagnosticField): +class Vorticity(DiagnosticField): u""" Base diagnostic Field for three dimensional Vorticity """ def setup(self, domain, state_fields, vorticity_type=None): @@ -1734,9 +1737,9 @@ def setup(self, domain, state_fields, vorticity_type=None): self.evaluator = LinearVariationalSolver(problem, solver_parameters={'ksp_type': 'cg'}) -class CompressibleRelativeVorticity(CompressibleVorticity): +class RelativeVorticity(Vorticity): u""" Diagnostic field for compressible euler relative vorticity """ - name = 'CompressibleRelativeVorticity' + name = 'RelativeVorticity' def __init__(self, space=None, method='solve'): u""" @@ -1762,9 +1765,9 @@ def setup(self, domain, state_fields): super().setup(domain, state_fields, vorticity_type='relative') -class CompressibleAbsoluteVorticity(CompressibleVorticity): +class AbsoluteVorticity(Vorticity): u""" Diagnostic field for compressible euler absolute vorticity """ - name = 'CompressibleAbsoluteVorticity' + name = 'AbsoluteVorticity' def __init__(self, parameters, space=None, method='solve'): u""" From e4babe562f3e775c6b76d68f844dc3aa8032da29 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Wed, 1 May 2024 10:52:49 +0100 Subject: [PATCH 03/17] fixed merge problems and flake8 --- gusto/diagnostics.py | 8 ++++---- gusto/equations.py | 12 +++++++----- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/gusto/diagnostics.py b/gusto/diagnostics.py index e4423411d..00e377c4e 100644 --- a/gusto/diagnostics.py +++ b/gusto/diagnostics.py @@ -28,7 +28,7 @@ "Perturbation", "Theta_e", "InternalEnergy", "PotentialEnergy", "ThermodynamicKineticEnergy", "Dewpoint", "Temperature", "Theta_d", "RelativeHumidity", "Pressure", "Exner_Vt", "HydrostaticImbalance", "Precipitation", - "ShallowWaterPotentialVorticity", "ShallowWaterRelativeVorticity", "ShallowWaterAbsoluteVorticity", + "ShallowWaterPotentialVorticity", "ShallowWaterRelativeVorticity", "ShallowWaterAbsoluteVorticity", "Vorticity", "RelativeVorticity", "AbsoluteVorticity", "Divergence", "BruntVaisalaFrequencySquared", "TracerDensity"] @@ -812,8 +812,8 @@ def __init__(self, base_field_name="ShallowWaterPotentialVorticity", space=None, for this diagnostic. Valid options are 'interpolate', 'project', 'assign' and 'solve'. Defaults to 'interpolate'. """ - base_enstrophy_names = ["ShallowWaterPotentialVorticity", - "ShallowWaterRelativeVorticity", + base_enstrophy_names = ["ShallowWaterPotentialVorticity", + "ShallowWaterRelativeVorticity", "ShallowWaterAbsoluteVorticity"] if base_field_name not in base_enstrophy_names: raise ValueError( @@ -1740,7 +1740,7 @@ def setup(self, domain, state_fields, vorticity_type=None): class RelativeVorticity(Vorticity): - u""" Diagnostic field for compressible euler relative vorticity """ + u""" Diagnostic field for compressible relative vorticity """ name = 'RelativeVorticity' def __init__(self, space=None, method='solve'): diff --git a/gusto/equations.py b/gusto/equations.py index b7c9b82c2..23b01292d 100644 --- a/gusto/equations.py +++ b/gusto/equations.py @@ -1074,7 +1074,7 @@ class CompressibleEulerEquations(PrognosticEquationSet): pressure. """ - def __init__(self, domain, parameters, Omega=None, sponge_options=None, + def __init__(self, domain, parameters, sponge_options=None, extra_terms=None, space_names=None, linearisation_map='default', u_transport_option="vector_invariant_form", @@ -1256,10 +1256,12 @@ def __init__(self, domain, parameters, Omega=None, sponge_options=None, # -------------------------------------------------------------------- # # Extra Terms (Coriolis, Sponge, Diffusion and others) # -------------------------------------------------------------------- # - if Omega is not None: - # TODO: add linearisation and label for this - residual += subject(prognostic( - inner(w, cross(2*Omega, u))*dx, "u"), self.X) + if parameters.Omega is not None: + # TODO add linerisation and label for this + Omega = as_vector((0, 0, parameters.Omega)) + coriolis_form = coriolis(subject(prognostic( + inner(w, cross(2*Omega, u))*dx, "u"), self.X)) + residual += coriolis_form if sponge_options is not None: W_DG = FunctionSpace(domain.mesh, "DG", 2) From 4afb9f201627c0cce1bf2afd0486b3d5abf04287 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Wed, 1 May 2024 14:21:42 +0100 Subject: [PATCH 04/17] added rudimentary vertical slice version --- gusto/diagnostics.py | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/gusto/diagnostics.py b/gusto/diagnostics.py index 00e377c4e..edd2b50de 100644 --- a/gusto/diagnostics.py +++ b/gusto/diagnostics.py @@ -1726,14 +1726,22 @@ def setup(self, domain, state_fields, vorticity_type=None): super().setup(domain, state_fields, space=space) if self.method == 'solve': - vort = TrialFunction(space) - w = TestFunction(space) - n = FacetNormal(domain.mesh) - a = inner(vort, w) * dx - L = inner(u, curl(w)) * dx - jump(cross(w, u), n) * dS_h - if vorticity_type == 'absolute': - Omega = as_vector((0, 0, self.parameters.Omega)) - L += inner(2*Omega, w) * dx + if domain.mesh.topological_dimension() == 3: + vort = TrialFunction(space) + w = TestFunction(space) + n = FacetNormal(domain.mesh) + a = inner(vort, w) * dx + L = inner(u, curl(w)) * dx - jump(cross(w, u), n) * dS_h + if vorticity_type == 'absolute': + Omega = as_vector((0, 0, self.parameters.Omega)) + L += inner(2*Omega, w) * dx + else: + space = domain.spaces("H1") + vort = TrialFunction(space) + gamma = TestFunction(space) + a = inner(vort, gamma) * dx + L = (- inner(domain.perp(grad(gamma)), u))*dx + # TODO implement absolute version, unsure atm how to get corioilis in vertical slice smartly problem = LinearVariationalProblem(a, L, self.field) self.evaluator = LinearVariationalSolver(problem, solver_parameters={'ksp_type': 'cg'}) @@ -1754,7 +1762,7 @@ def __init__(self, space=None, method='solve'): 'assign' and 'solve'. Defaults to 'solve'. """ self.solve_implemented = True - super().__init__(space=space, method=method, required_fields=('u',)) + super().__init__(space=space, method=method, required_fields=('u')) def setup(self, domain, state_fields): u""" From d3c9381f8ca21d0b109cc98da3ff1e8aabf67086 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Tue, 7 May 2024 09:11:18 +0100 Subject: [PATCH 05/17] added source term for vertical slice --- gusto/diagnostics.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/gusto/diagnostics.py b/gusto/diagnostics.py index edd2b50de..a6fd26009 100644 --- a/gusto/diagnostics.py +++ b/gusto/diagnostics.py @@ -1713,8 +1713,11 @@ def setup(self, domain, state_fields, vorticity_type=None): raise ValueError('Potential vorticity has not yet been implemented') else: raise ValueError(f'vorticity type must be one of {vorticity_types}, not {vorticity_type}') + if domain.mesh.topological_dimension() == 3: + space = domain.spaces('HCurl') + else: + space = domain.spaces('H1') - space = domain.spaces('HCurl') u = state_fields('u') if self.method != 'solve': if vorticity_type == 'relative': @@ -1736,13 +1739,12 @@ def setup(self, domain, state_fields, vorticity_type=None): Omega = as_vector((0, 0, self.parameters.Omega)) L += inner(2*Omega, w) * dx else: - space = domain.spaces("H1") vort = TrialFunction(space) gamma = TestFunction(space) + n = FacetNormal(domain.mesh) a = inner(vort, gamma) * dx - L = (- inner(domain.perp(grad(gamma)), u))*dx + L = ( inner(domain.perp(grad(gamma)), u))*dx #- jump(inner(n, u)*gamma)*ds_v # TODO implement absolute version, unsure atm how to get corioilis in vertical slice smartly - problem = LinearVariationalProblem(a, L, self.field) self.evaluator = LinearVariationalSolver(problem, solver_parameters={'ksp_type': 'cg'}) From d76254478f3f6ab58d6730400223cc70e44932ee Mon Sep 17 00:00:00 2001 From: Witt-D Date: Tue, 7 May 2024 10:09:41 +0100 Subject: [PATCH 06/17] added perp --- gusto/diagnostics.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/gusto/diagnostics.py b/gusto/diagnostics.py index a6fd26009..e8d87d66c 100644 --- a/gusto/diagnostics.py +++ b/gusto/diagnostics.py @@ -1743,12 +1743,11 @@ def setup(self, domain, state_fields, vorticity_type=None): gamma = TestFunction(space) n = FacetNormal(domain.mesh) a = inner(vort, gamma) * dx - L = ( inner(domain.perp(grad(gamma)), u))*dx #- jump(inner(n, u)*gamma)*ds_v + L = ( inner(domain.perp(grad(gamma)), u))*dx - jump(inner(domain.perp(n), u)*gamma)*ds_h # TODO implement absolute version, unsure atm how to get corioilis in vertical slice smartly problem = LinearVariationalProblem(a, L, self.field) self.evaluator = LinearVariationalSolver(problem, solver_parameters={'ksp_type': 'cg'}) - class RelativeVorticity(Vorticity): u""" Diagnostic field for compressible relative vorticity """ name = 'RelativeVorticity' From 6cb1116e6adb41cec3fc9af95eea03db91390461 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Wed, 29 May 2024 09:30:42 +0100 Subject: [PATCH 07/17] saving current state --- gusto/diagnostics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gusto/diagnostics.py b/gusto/diagnostics.py index e8d87d66c..1a5e26b33 100644 --- a/gusto/diagnostics.py +++ b/gusto/diagnostics.py @@ -1743,7 +1743,7 @@ def setup(self, domain, state_fields, vorticity_type=None): gamma = TestFunction(space) n = FacetNormal(domain.mesh) a = inner(vort, gamma) * dx - L = ( inner(domain.perp(grad(gamma)), u))*dx - jump(inner(domain.perp(n), u)*gamma)*ds_h + L = -( inner(domain.perp(grad(gamma)), u)) * dx #- jump(inner(domain.perp(n), u) * gamma) * (dS_h + dS_v) # TODO implement absolute version, unsure atm how to get corioilis in vertical slice smartly problem = LinearVariationalProblem(a, L, self.field) self.evaluator = LinearVariationalSolver(problem, solver_parameters={'ksp_type': 'cg'}) From e5eabd51aa913531fc7677c077558dcca78c3af8 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Wed, 29 May 2024 11:05:22 +0100 Subject: [PATCH 08/17] edited diagnostic --- gusto/diagnostics.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/gusto/diagnostics.py b/gusto/diagnostics.py index 1a5e26b33..8917d613e 100644 --- a/gusto/diagnostics.py +++ b/gusto/diagnostics.py @@ -1734,16 +1734,15 @@ def setup(self, domain, state_fields, vorticity_type=None): w = TestFunction(space) n = FacetNormal(domain.mesh) a = inner(vort, w) * dx - L = inner(u, curl(w)) * dx - jump(cross(w, u), n) * dS_h + L = inner(u, curl(w)) * dx - jump(cross(w, u)) * (dS_h + dS_v) if vorticity_type == 'absolute': Omega = as_vector((0, 0, self.parameters.Omega)) L += inner(2*Omega, w) * dx else: vort = TrialFunction(space) gamma = TestFunction(space) - n = FacetNormal(domain.mesh) a = inner(vort, gamma) * dx - L = -( inner(domain.perp(grad(gamma)), u)) * dx #- jump(inner(domain.perp(n), u) * gamma) * (dS_h + dS_v) + L = -( inner(domain.perp(grad(gamma)), u)) * dx # TODO implement absolute version, unsure atm how to get corioilis in vertical slice smartly problem = LinearVariationalProblem(a, L, self.field) self.evaluator = LinearVariationalSolver(problem, solver_parameters={'ksp_type': 'cg'}) From 9037ea4eadf6537020acf7a151b66e44cb22f0f4 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Mon, 3 Jun 2024 15:57:27 +0100 Subject: [PATCH 09/17] added unit test --- unit-tests/test_vorticity_diagnostic.py | 45 +++++++++++++++++++++++++ 1 file changed, 45 insertions(+) create mode 100644 unit-tests/test_vorticity_diagnostic.py diff --git a/unit-tests/test_vorticity_diagnostic.py b/unit-tests/test_vorticity_diagnostic.py new file mode 100644 index 000000000..ce90a5324 --- /dev/null +++ b/unit-tests/test_vorticity_diagnostic.py @@ -0,0 +1,45 @@ + +from gusto.diagnostics import RelativeVorticity +from gusto.fields import StateFields, PrescribedFields, TimeLevelFields +from gusto import Domain, CompressibleParameters, CompressibleEulerEquations +from firedrake import (PeriodicIntervalMesh, ExtrudedMesh, Function, sin, cos, + SpatialCoordinate, pi, as_vector, errornorm, norm) + + +def test_vorticity(): + L = 10 + H = 10 + ncol = 100 + nlayers = 100 + + m = PeriodicIntervalMesh(ncol, L) + mesh = ExtrudedMesh(m, layers=nlayers, layer_height=H/nlayers) + _, z = SpatialCoordinate(mesh) + + domain = Domain(mesh, 0.1, 'CG', 1) + params = CompressibleParameters() + eqn = CompressibleEulerEquations(domain, params) + prog_field = TimeLevelFields(eqn) + + H1 = domain.spaces('H1') + HDiv = domain.spaces('HDiv') + + u_expr = 3 * sin(2*pi*z/H) + vort_exact_expr = -6*pi/H * cos(2*pi*z/H) + vorticity_analytic = Function(H1, name='analytic_vort').interpolate(vort_exact_expr) + + # Setting up test field for the diagnostic to use + prescribed_fields = PrescribedFields() + prescribed_fields('u', HDiv) + state = StateFields(prog_field, prescribed_fields) + state.u.project(as_vector([u_expr, 0])) + + Vorticity = RelativeVorticity() + Vorticity.setup(domain, state) + Vorticity.compute() + # Compare analytic vorticity expression to diagnostic + error = errornorm(vorticity_analytic, state.RelativeVorticity) / norm(vorticity_analytic) + print(error) + # We dont expect it to be zero as the discrete vorticity is not equal to analytic and dependent on resolution + assert error < 1e-6, \ + 'Relative Vorticity not in error tolerence' From e616591258f7768040b9d580e1cca8f11a98a172 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Tue, 11 Jun 2024 11:07:34 +0100 Subject: [PATCH 10/17] added current state of vorticity diagnostic --- gusto/diagnostics.py | 13 +-- unit-tests/test_vorticity_diagnostic.py | 100 ++++++++++++++++++++++-- 2 files changed, 102 insertions(+), 11 deletions(-) diff --git a/gusto/diagnostics.py b/gusto/diagnostics.py index 8917d613e..c132eb741 100644 --- a/gusto/diagnostics.py +++ b/gusto/diagnostics.py @@ -1731,22 +1731,23 @@ def setup(self, domain, state_fields, vorticity_type=None): if self.method == 'solve': if domain.mesh.topological_dimension() == 3: vort = TrialFunction(space) - w = TestFunction(space) - n = FacetNormal(domain.mesh) - a = inner(vort, w) * dx - L = inner(u, curl(w)) * dx - jump(cross(w, u)) * (dS_h + dS_v) + gamma = TestFunction(space) + a = inner(vort, gamma) * dx + L = -inner(curl(gamma), u) * dx if vorticity_type == 'absolute': Omega = as_vector((0, 0, self.parameters.Omega)) - L += inner(2*Omega, w) * dx + L += inner(2*Omega, gamma) * dx else: vort = TrialFunction(space) gamma = TestFunction(space) a = inner(vort, gamma) * dx - L = -( inner(domain.perp(grad(gamma)), u)) * dx + L = -(inner(domain.perp(grad(gamma)), u)) * dx # TODO implement absolute version, unsure atm how to get corioilis in vertical slice smartly + problem = LinearVariationalProblem(a, L, self.field) self.evaluator = LinearVariationalSolver(problem, solver_parameters={'ksp_type': 'cg'}) + class RelativeVorticity(Vorticity): u""" Diagnostic field for compressible relative vorticity """ name = 'RelativeVorticity' diff --git a/unit-tests/test_vorticity_diagnostic.py b/unit-tests/test_vorticity_diagnostic.py index ce90a5324..f5f80e1bd 100644 --- a/unit-tests/test_vorticity_diagnostic.py +++ b/unit-tests/test_vorticity_diagnostic.py @@ -1,12 +1,101 @@ -from gusto.diagnostics import RelativeVorticity +from gusto.diagnostics import RelativeVorticity, ZonalComponent, MeridionalComponent, RadialComponent from gusto.fields import StateFields, PrescribedFields, TimeLevelFields -from gusto import Domain, CompressibleParameters, CompressibleEulerEquations -from firedrake import (PeriodicIntervalMesh, ExtrudedMesh, Function, sin, cos, +from gusto import (Domain, CompressibleParameters, CompressibleEulerEquations, + GeneralCubedSphereMesh, lonlatr_from_xyz, xyz_vector_from_lonlatr) +from firedrake import (PeriodicIntervalMesh, ExtrudedMesh, Function, sin, cos, File, SpatialCoordinate, pi, as_vector, errornorm, norm) +import pytest -def test_vorticity(): +@pytest.mark.parametrize("topology", ["slice", "sphere"]) +def test_relative_vorticity(topology): + if topology=='slice': + vertical_slice_test() + + if topology=='sphere': + sphere_test() + + +def sphere_test(): + R = 1 # radius of ball + H = 5 # height of model top + nlayers=5 + c=8 + # Building model and state object + m = GeneralCubedSphereMesh(R, num_cells_per_edge_of_panel=c, degree=2) + mesh = ExtrudedMesh(m, layers=nlayers, layer_height=H/nlayers, extrusion_type='radial') + domain=Domain(mesh, 0.1, 'RTCF', degree=1) + params = CompressibleParameters(g=0, cp=0) + eqn = CompressibleEulerEquations(domain, params) + prog_field = TimeLevelFields(eqn) + HDiv = domain.spaces('HDiv') + HCurl = domain.spaces('HCurl') + pres_field = PrescribedFields() + pres_field('u', HDiv) + state = StateFields(prog_field, pres_field) + + # Getting spherical co-ordinates + xyz = SpatialCoordinate(mesh) + _, lat, r = lonlatr_from_xyz(xyz[0], xyz[1], xyz[2]) + e_zonal = xyz_vector_from_lonlatr(1, 0, 0, xyz) + e_rad = xyz_vector_from_lonlatr(0, 0, 1, xyz) + + + # Initlising vorticity field + zonal_u = sin(lat) / r + merid_u = 0.0 + radial_u = 0.0 + #u_expr = xyz_vector_from_lonlatr(zonal_u, merid_u, radial_u, xyz) + print('Projecting u') + state.u.project(e_zonal*zonal_u) + + # Analytic relative vorticity + radial_vort = 2*cos(lat) / r**2 + analytical_vort_expr = xyz_vector_from_lonlatr(0, 0, radial_vort, xyz) + print('Projecting analytical vorticity') + vorticity_analytic = Function(HCurl, name='exact_vorticity').project(e_rad*radial_vort) + + # Setup and compute diagnostic vorticity and zonal components + Vorticity = RelativeVorticity() + Vorticity.setup(domain, state) + print('diagnosing vorticity') + Vorticity.compute() + + Zonal_u = ZonalComponent('u') + Zonal_u.setup(domain, state) + Zonal_u.compute() + + Meridional_u = MeridionalComponent('u') + Meridional_u.setup(domain, state) + Meridional_u.compute() + + Radial_u = RadialComponent('u') + Radial_u.setup(domain, state) + Radial_u.compute() + + + # zonal_diff = Function(HDiv, name='zonal_diff').project(state.u - state.u_zonal) + + diff = Function(HCurl, name='difference').project(vorticity_analytic - state.RelativeVorticity) + + + radial_diagnostic_vort = RadialComponent('RelativeVorticity') + radial_diagnostic_vort.setup(domain, state) + radial_diagnostic_vort.compute() + + # Compare analytic vorticity expression to diagnostic + print('calculating error') + error = errornorm(vorticity_analytic, state.RelativeVorticity) / norm(vorticity_analytic) + print(error) + + outfile = File('spherical_vorticity.pvd') + outfile.write(state.u, vorticity_analytic, state.RelativeVorticity, diff, state.u_zonal, state.u_meridional, state.u_radial, state.RelativeVorticity_radial) + # We dont expect it to be zero as the discrete vorticity is not equal to analytic and dependent on resolution + assert error < 1e-6, \ + 'Relative Vorticity not in error tolerence' + +def vertical_slice_test(): L = 10 H = 10 ncol = 100 @@ -39,7 +128,8 @@ def test_vorticity(): Vorticity.compute() # Compare analytic vorticity expression to diagnostic error = errornorm(vorticity_analytic, state.RelativeVorticity) / norm(vorticity_analytic) - print(error) # We dont expect it to be zero as the discrete vorticity is not equal to analytic and dependent on resolution assert error < 1e-6, \ 'Relative Vorticity not in error tolerence' + +sphere_test() \ No newline at end of file From 7dc964b5f2211d8798a98c9f3a580a32def900f3 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Tue, 23 Jul 2024 11:03:24 +0100 Subject: [PATCH 11/17] re-instated vorticity changes after code refactor --- gusto/core/configuration.py | 1 + .../compressible_euler_diagnostics.py | 114 ++++++++++++++++++ .../equations/compressible_euler_equations.py | 13 +- 3 files changed, 123 insertions(+), 5 deletions(-) diff --git a/gusto/core/configuration.py b/gusto/core/configuration.py index 5aebf5a33..9abb74647 100644 --- a/gusto/core/configuration.py +++ b/gusto/core/configuration.py @@ -136,6 +136,7 @@ class CompressibleParameters(Configuration): w_sat2 = -17.27 # second const. in Teten's formula (no units) w_sat3 = 35.86 # third const. in Teten's formula (K) w_sat4 = 610.9 # fourth const. in Teten's formula (Pa) + Omega = None # Rotation rate class ShallowWaterParameters(Configuration): diff --git a/gusto/diagnostics/compressible_euler_diagnostics.py b/gusto/diagnostics/compressible_euler_diagnostics.py index ce03468d6..f66ae4d76 100644 --- a/gusto/diagnostics/compressible_euler_diagnostics.py +++ b/gusto/diagnostics/compressible_euler_diagnostics.py @@ -651,3 +651,117 @@ def compute(self): """Increment the precipitation diagnostic.""" self.solver.solve() self.field.assign(self.field + self.flux) + +class Vorticity(DiagnosticField): + u""" Base diagnostic Field for three dimensional Vorticity """ + + def setup(self, domain, state_fields, vorticity_type=None): + """ + Sets up the :class:`Function` for the diagnostic field. + + Args: + domain (:class:`Domain`): the model's domain object. + state_fields (:class:`StateFields`): the model's field container. + vorticity_type (str, optional): denotes which type of vorticity to + be computed ('relative', 'absolute'). Defaults to + None. + """ + # TODO Do we eventually want a potential voriticy? + vorticity_types = ['relative', 'absolute'] + if vorticity_type not in vorticity_types: + if vorticity_type == 'potential': + raise ValueError('Potential vorticity has not yet been implemented') + else: + raise ValueError(f'vorticity type must be one of {vorticity_types}, not {vorticity_type}') + if domain.mesh.topological_dimension() == 3: + space = domain.spaces('HCurl') + else: + space = domain.spaces('H1') + + u = state_fields('u') + if self.method != 'solve': + if vorticity_type == 'relative': + self.expression = curl(u) + elif vorticity_type == 'absolute': + Omega = as_vector((0, 0, self.parameters.Omega)) + self.expression = curl(u) + 2*Omega + + super().setup(domain, state_fields, space=space) + + if self.method == 'solve': + if domain.mesh.topological_dimension() == 3: + vort = TrialFunction(space) + gamma = TestFunction(space) + a = inner(vort, gamma) * dx + L = -inner(curl(gamma), u) * dx + if vorticity_type == 'absolute': + Omega = as_vector((0, 0, self.parameters.Omega)) + L += inner(2*Omega, gamma) * dx + else: + vort = TrialFunction(space) + gamma = TestFunction(space) + a = inner(vort, gamma) * dx + L = -(inner(domain.perp(grad(gamma)), u)) * dx + # TODO implement absolute version, unsure atm how to get corioilis in vertical slice smartly + + problem = LinearVariationalProblem(a, L, self.field) + self.evaluator = LinearVariationalSolver(problem, solver_parameters={'ksp_type': 'cg'}) + + +class RelativeVorticity(Vorticity): + u""" Diagnostic field for compressible relative vorticity """ + name = 'RelativeVorticity' + + def __init__(self, space=None, method='solve'): + u""" + Args: + space (:class:`FunctionSpace`, optional): the function space to + evaluate the diagnostic field in. Defaults to None, in which + case a default space will be chosen for this diagnostic. + method (str, optional): a string specifying the method of evaluation + for this diagnostic. Valid options are 'interpolate', 'project', + 'assign' and 'solve'. Defaults to 'solve'. + """ + self.solve_implemented = True + super().__init__(space=space, method=method, required_fields=('u')) + + def setup(self, domain, state_fields): + u""" + Sets up the :class:`Function` for the diagnostic field. + + Args: + domain (:class:`Domain`): the model's domain object. + state_fields (:class:`StateFields`): the model's field container. + """ + super().setup(domain, state_fields, vorticity_type='relative') + + +class AbsoluteVorticity(Vorticity): + u""" Diagnostic field for compressible euler absolute vorticity """ + name = 'AbsoluteVorticity' + + def __init__(self, parameters, space=None, method='solve'): + u""" + Args: + parameters (:class:`CompressibleEulerParameters`): the configuration + object containing the physical parameters for this equation. + space (:class:`FunctionSpace`, optional): the function space to + evaluate the diagnostic field in. Defaults to None, in which + case a default space will be chosen for this diagnostic. + method (str, optional): a string specifying the method of evaluation + for this diagnostic. Valid options are 'interpolate', 'project', + 'assign' and 'solve'. Defaults to 'solve'. + """ + self.solve_implemented = True + self.parameters = parameters + super().__init__(space=space, method=method, required_fields=('u')) + + def setup(self, domain, state_fields): + u""" + Sets up the :class:`Function` for the diagnostic field. + + Args: + domain (:class:`Domain`): the model's domain object. + state_fields (:class:`StateFields`): the model's field container. + """ + super().setup(domain, state_fields, vorticity_type='absolute') diff --git a/gusto/equations/compressible_euler_equations.py b/gusto/equations/compressible_euler_equations.py index b6b89cd8c..80dcf48d8 100644 --- a/gusto/equations/compressible_euler_equations.py +++ b/gusto/equations/compressible_euler_equations.py @@ -2,7 +2,7 @@ from firedrake import ( sin, pi, inner, dx, div, cross, FunctionSpace, FacetNormal, jump, avg, dS_v, - conditional, SpatialCoordinate, split, Constant + conditional, SpatialCoordinate, split, Constant, as_vector ) from firedrake.fml import subject, replace_subject from gusto.core.labels import ( @@ -34,7 +34,7 @@ class CompressibleEulerEquations(PrognosticEquationSet): pressure. """ - def __init__(self, domain, parameters, Omega=None, sponge_options=None, + def __init__(self, domain, parameters, sponge_options=None, extra_terms=None, space_names=None, linearisation_map='default', u_transport_option="vector_invariant_form", @@ -216,10 +216,13 @@ def __init__(self, domain, parameters, Omega=None, sponge_options=None, # -------------------------------------------------------------------- # # Extra Terms (Coriolis, Sponge, Diffusion and others) # -------------------------------------------------------------------- # - if Omega is not None: - # TODO: add linearisation - residual += coriolis(subject(prognostic( + if parameters.Omega is not None: + # TODO add linerisation and label for this + Omega = as_vector((0, 0, parameters.Omega)) + coriolis_form = coriolis(subject(prognostic( inner(w, cross(2*Omega, u))*dx, "u"), self.X)) + residual += coriolis_form + if sponge_options is not None: W_DG = FunctionSpace(domain.mesh, "DG", 2) From 0d52733e0520cd809ce0e5a9bda68c175f3bbe7f Mon Sep 17 00:00:00 2001 From: Witt-D Date: Mon, 19 Aug 2024 09:33:43 +0100 Subject: [PATCH 12/17] testing break --- unit-tests/test_vorticity_diagnostic.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/unit-tests/test_vorticity_diagnostic.py b/unit-tests/test_vorticity_diagnostic.py index f5f80e1bd..4e4c8f3f9 100644 --- a/unit-tests/test_vorticity_diagnostic.py +++ b/unit-tests/test_vorticity_diagnostic.py @@ -1,6 +1,6 @@ from gusto.diagnostics import RelativeVorticity, ZonalComponent, MeridionalComponent, RadialComponent -from gusto.fields import StateFields, PrescribedFields, TimeLevelFields +from gusto.core.fields import StateFields, PrescribedFields, TimeLevelFields from gusto import (Domain, CompressibleParameters, CompressibleEulerEquations, GeneralCubedSphereMesh, lonlatr_from_xyz, xyz_vector_from_lonlatr) from firedrake import (PeriodicIntervalMesh, ExtrudedMesh, Function, sin, cos, File, @@ -17,7 +17,7 @@ def test_relative_vorticity(topology): sphere_test() -def sphere_test(): +def sphere_test(): # try at a higher resolution R = 1 # radius of ball H = 5 # height of model top nlayers=5 @@ -76,7 +76,6 @@ def sphere_test(): # zonal_diff = Function(HDiv, name='zonal_diff').project(state.u - state.u_zonal) - diff = Function(HCurl, name='difference').project(vorticity_analytic - state.RelativeVorticity) From 6aacb34a1290c9371f7da66b28bdeb90cf6de184 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Mon, 19 Aug 2024 10:34:18 +0100 Subject: [PATCH 13/17] curremt status --- unit-tests/test_vorticity_diagnostic.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/unit-tests/test_vorticity_diagnostic.py b/unit-tests/test_vorticity_diagnostic.py index 4e4c8f3f9..9855701f7 100644 --- a/unit-tests/test_vorticity_diagnostic.py +++ b/unit-tests/test_vorticity_diagnostic.py @@ -21,7 +21,7 @@ def sphere_test(): # try at a higher resolution R = 1 # radius of ball H = 5 # height of model top nlayers=5 - c=8 + c=16 # Building model and state object m = GeneralCubedSphereMesh(R, num_cells_per_edge_of_panel=c, degree=2) mesh = ExtrudedMesh(m, layers=nlayers, layer_height=H/nlayers, extrusion_type='radial') @@ -52,7 +52,7 @@ def sphere_test(): # try at a higher resolution # Analytic relative vorticity radial_vort = 2*cos(lat) / r**2 - analytical_vort_expr = xyz_vector_from_lonlatr(0, 0, radial_vort, xyz) + #analytical_vort_expr = xyz_vector_from_lonlatr(0, 0, radial_vort, xyz) print('Projecting analytical vorticity') vorticity_analytic = Function(HCurl, name='exact_vorticity').project(e_rad*radial_vort) @@ -88,7 +88,7 @@ def sphere_test(): # try at a higher resolution error = errornorm(vorticity_analytic, state.RelativeVorticity) / norm(vorticity_analytic) print(error) - outfile = File('spherical_vorticity.pvd') + outfile = File('spherical_vorticity_high.pvd') outfile.write(state.u, vorticity_analytic, state.RelativeVorticity, diff, state.u_zonal, state.u_meridional, state.u_radial, state.RelativeVorticity_radial) # We dont expect it to be zero as the discrete vorticity is not equal to analytic and dependent on resolution assert error < 1e-6, \ From 3a66723315fc1dc699610f1a44b84e82bcd2fbef Mon Sep 17 00:00:00 2001 From: Witt-D Date: Mon, 19 Aug 2024 10:35:11 +0100 Subject: [PATCH 14/17] some adjustments to the diagnostics --- gusto/diagnostics/balance.py | 180 ++++++++++++++++++ .../compressible_euler_diagnostics.py | 6 +- .../diagnostics/shallow_water_diagnostics.py | 18 +- 3 files changed, 192 insertions(+), 12 deletions(-) create mode 100644 gusto/diagnostics/balance.py diff --git a/gusto/diagnostics/balance.py b/gusto/diagnostics/balance.py new file mode 100644 index 000000000..6062ac53d --- /dev/null +++ b/gusto/diagnostics/balance.py @@ -0,0 +1,180 @@ +from firedrake import (Constant, Function, dx, TestFunction, TrialFunction, as_vector, + inner, dot, cross, div, jump, LinearVariationalProblem, + LinearVariationalSolver, SpatialCoordinate, tan, FacetNormal, + sqrt, avg, dS_v) + +from gusto.core.coord_transforms import lonlatr_from_xyz +from gusto.diagnostics.diagnostics import DiagnosticField +from gusto.equations.compressible_euler_equations import CompressibleEulerEquations +from gusto.equations.thermodynamics import exner_pressure + +class GeostrophicImbalance(DiagnosticField): + """Geostrophic imbalance diagnostic field.""" + name = "GeostrophicImbalance" + + def __init__(self, equations, space=None, method='interpolate'): + """ + Args: + equations (:class:`PrognosticEquationSet`): the equation set being + solved by the model. + space (:class:`FunctionSpace`, optional): the function space to + evaluate the diagnostic field in. Defaults to None, in which + case a default space will be chosen for this diagnostic. + method (str, optional): a string specifying the method of evaluation + for this diagnostic. Valid options are 'interpolate', 'project', + 'assign' and 'solve'. Defaults to 'interpolate'. + """ + # Work out required fields + if isinstance(equations, CompressibleEulerEquations): + required_fields = ['rho', 'theta', 'u'] + self.equations = equations + self.parameters = equations.parameters + else: + raise NotImplementedError(f'Geostrophic Imbalance not implemented for {type(equations)}') + super().__init__(space=space, method=method, required_fields=required_fields) + + def setup(self, domain, state_fields): + """ + Sets up the :class:`Function` for the diagnostic field. + + Args: + domain (:class:`Domain`): the model's domain object. + state_fields (:class:`StateFields`): the model's field container. + """ + Vcurl = domain.spaces("HCurl") + u = state_fields('u') + rho = state_fields('rho') + theta = state_fields('theta') + + exner = exner_pressure(self.parameters, rho, theta) + cp = Constant(self.parameters.cp) + n = FacetNormal(domain.mesh) + # TODO: Generilise this for cases that aren't solid body rotation case + omega = Constant(7.292e-5) + Omega = as_vector((0., 0., omega)) + k = domain.k + + F = TrialFunction(Vcurl) + w = TestFunction(Vcurl) + + imbalance = Function(Vcurl) + a = inner(w, F)*dx + + L = (- cp*div(theta*w)*exner*dx + cp*jump(theta*w, n)*avg(exner)*dS_v # exner pressure grad discretisation + - inner(w, cross(2*Omega, u))*dx # coriolis + + inner(w, dot(k, cross(2*Omega, u) )*k)*dx #vertical part of coriolis + + cp*div((theta*k)*dot(k,w))*exner*dx # removes vertical part of the pressure divergence + - cp*jump((theta*k)*dot(k,w), n)*avg(exner)*dS_v) # removes vertical part of pressure jump condition + + + bcs = self.equations.bcs['u'] + + imbalanceproblem = LinearVariationalProblem(a, L, imbalance, bcs=bcs) + self.imbalance_solver = LinearVariationalSolver(imbalanceproblem) + self.expr = dot(imbalance, domain.k) + super().setup(domain, state_fields) + + def compute(self): + """Compute and return the diagnostic field from the current state. + """ + self.imbalance_solver.solve() + super().compute() + + +class SolidBodyImbalance(DiagnosticField): + """Solid Body imbalance diagnostic field.""" + name = "SolidBodyImbalance" + + def __init__(self, equations, space=None, method='interpolate'): + """ + Args: + equations (:class:`PrognosticEquationSet`): the equation set being + solved by the model. + space (:class:`FunctionSpace`, optional): the function space to + evaluate the diagnostic field in. Defaults to None, in which + case a default space will be chosen for this diagnostic. + method (str, optional): a string specifying the method of evaluation + for this diagnostic. Valid options are 'interpolate', 'project', + 'assign' and 'solve'. Defaults to 'interpolate'. + """ + # Work out required fields + if isinstance(equations, CompressibleEulerEquations): + required_fields = ['rho', 'theta', 'u'] + self.equations = equations + self.parameters = equations.parameters + else: + raise NotImplementedError(f'Geostrophic Imbalance not implemented for {type(equations)}') + super().__init__(space=space, method=method, required_fields=required_fields) + + def setup(self, domain, state_fields): + """ + Sets up the :class:`Function` for the diagnostic field. + + Args: + domain (:class:`Domain`): the model's domain object. + state_fields (:class:`StateFields`): the model's field container. + """ + Vu = domain.spaces("HDiv") + u = state_fields('u') + rho = state_fields('rho') + theta = state_fields('theta') + + exner = tde.exner_pressure(self.parameters, rho, theta) + cp = Constant(self.parameters.cp) + n = FacetNormal(domain.mesh) + # TODO: Generilise this for cases that aren't solid body rotation case + omega = Constant(7.292e-5) + Omega = as_vector((0., 0., omega)) + + # generating the spherical co-ords, and spherical components of velocity + x, y, z = SpatialCoordinate(domain.mesh) + x_hat = Constant(as_vector([1.0, 0.0, 0.0])) + y_hat = Constant(as_vector([0.0, 1.0, 0.0])) + z_hat = Constant(as_vector([0.0, 0.0, 1.0])) + R = sqrt(x**2 + y**2) # distance from z axis + r = sqrt(x**2 + y**2 + z**2) # distance from origin + lambda_hat = (x * y_hat - y * x_hat) / R + lon_dot = inner(u, lambda_hat) + phi_hat = (-x*z/R * x_hat - y*z/R * y_hat + R * z_hat) / r + lat_dot = inner(u, phi_hat) + r_hat = (x * x_hat + y * y_hat + z * z_hat) / r + r_dot = inner(u, r_hat) + mesh = domain.mesh + k=domain.k + #HACK loweing the quadrature manually + dx_low_quadrature = dx(degree=3) + lat = latlon_coords(mesh)[0] + F = TrialFunction(Vu) + w = TestFunction(Vu) + + imbalance = Function(Vu) + a = inner(w, F)*dx + + #TODO Segmentaion error somehwere here + + L = (- cp*div((theta)*w)*exner*dx + cp*jump((theta)*w, n)*avg(exner)*dS_v # exner pressure grad discretisation + - inner(w, cross(2*Omega, u))*dx # coriolis + + inner(w, dot(k, cross(2*Omega, u) )*k)*dx # vertical part of coriolis + + cp*div((theta*k)*dot(k,w))*exner*dx # vertical part of the pressure divergence + - cp*jump((theta*k)*dot(k,w), n)*avg(exner)*dS_v # vertical part of pressure jump condition + # BUG it is these terms arising from the non-linear that are the problem + - (lat_dot * lon_dot * tan(lat) / r)*inner(w, lambda_hat)*dx_low_quadrature + + (lon_dot * r_dot / r)*dot(w, lambda_hat)*dx_low_quadrature # lambda component of non linear term + + (lat_dot**2 * tan(lat) / r)*inner(w, phi_hat)*dx_low_quadrature # phi component 1 + + (lat_dot * r_dot / r)*inner(w, phi_hat)*dx_low_quadrature # phi component 1 + ) + + + bcs = self.equations.bcs['u'] + + imbalanceproblem = LinearVariationalProblem(a, L, imbalance, bcs=bcs) + self.imbalance_solver = LinearVariationalSolver(imbalanceproblem) + self.expr = dot(imbalance, domain.k) + super().setup(domain, state_fields) + + + def compute(self): + """Compute and return the diagnostic field from the current state. + """ + self.imbalance_solver.solve() + super().compute() \ No newline at end of file diff --git a/gusto/diagnostics/compressible_euler_diagnostics.py b/gusto/diagnostics/compressible_euler_diagnostics.py index f86db8240..62194672f 100644 --- a/gusto/diagnostics/compressible_euler_diagnostics.py +++ b/gusto/diagnostics/compressible_euler_diagnostics.py @@ -3,7 +3,7 @@ from firedrake import (dot, dx, Function, ln, TestFunction, TrialFunction, Constant, grad, inner, LinearVariationalProblem, LinearVariationalSolver, FacetNormal, ds_b, dS_v, div, - avg, jump, SpatialCoordinate) + avg, jump, SpatialCoordinate, curl, as_vector) from gusto.diagnostics.diagnostics import ( DiagnosticField, Energy, IterativeDiagnosticField @@ -18,7 +18,7 @@ "PotentialEnergy", "ThermodynamicKineticEnergy", "Temperature", "Theta_d", "RelativeHumidity", "Pressure", "Exner_Vt", "HydrostaticImbalance", "Precipitation", "BruntVaisalaFrequencySquared", - "WetBulbTemperature", "DewpointTemperature"] + "WetBulbTemperature", "DewpointTemperature", "RelativeVorticity", "AbsoluteVorticity"] class RichardsonNumber(DiagnosticField): @@ -852,7 +852,7 @@ def setup(self, domain, state_fields, vorticity_type=None): gamma = TestFunction(space) a = inner(vort, gamma) * dx L = -(inner(domain.perp(grad(gamma)), u)) * dx - # TODO implement absolute version, unsure atm how to get corioilis in vertical slice smartly + # TODO implement absolute version, unsure how to get corioilis in vertical slice smartly problem = LinearVariationalProblem(a, L, self.field) self.evaluator = LinearVariationalSolver(problem, solver_parameters={'ksp_type': 'cg'}) diff --git a/gusto/diagnostics/shallow_water_diagnostics.py b/gusto/diagnostics/shallow_water_diagnostics.py index dce6fac69..1ce4d58fc 100644 --- a/gusto/diagnostics/shallow_water_diagnostics.py +++ b/gusto/diagnostics/shallow_water_diagnostics.py @@ -6,8 +6,8 @@ from gusto.diagnostics.diagnostics import DiagnosticField, Energy __all__ = ["ShallowWaterKineticEnergy", "ShallowWaterPotentialEnergy", - "ShallowWaterPotentialEnstrophy", "PotentialVorticity", - "RelativeVorticity", "AbsoluteVorticity"] + "ShallowWaterPotentialEnstrophy", "ShallowWaterPotentialVorticity", + "ShallowWaterRelativeVorticity", "ShallowWaterAbsoluteVorticity"] class ShallowWaterKineticEnergy(Energy): @@ -136,7 +136,7 @@ def setup(self, domain, state_fields): super().setup(domain, state_fields) -class Vorticity(DiagnosticField): +class ShallowWaterVorticity(DiagnosticField): """Base diagnostic field class for shallow-water vorticity variables.""" def setup(self, domain, state_fields, vorticity_type=None): @@ -191,9 +191,9 @@ def setup(self, domain, state_fields, vorticity_type=None): self.evaluator = LinearVariationalSolver(problem, solver_parameters={"ksp_type": "cg"}) -class PotentialVorticity(Vorticity): +class ShallowWaterPotentialVorticity(ShallowWaterVorticity): u"""Diagnostic field for shallow-water potential vorticity, q=(∇×(u+f))/D""" - name = "PotentialVorticity" + name = "ShallowWaterPotentialVorticity" def __init__(self, space=None, method='solve'): """ @@ -220,9 +220,9 @@ def setup(self, domain, state_fields): super().setup(domain, state_fields, vorticity_type="potential") -class AbsoluteVorticity(Vorticity): +class ShallowWaterAbsoluteVorticity(ShallowWaterVorticity): u"""Diagnostic field for absolute vorticity, ζ=∇×(u+f)""" - name = "AbsoluteVorticity" + name = "ShallowWaterAbsoluteVorticity" def __init__(self, space=None, method='solve'): """ @@ -248,9 +248,9 @@ def setup(self, domain, state_fields): super().setup(domain, state_fields, vorticity_type="absolute") -class RelativeVorticity(Vorticity): +class ShallowWaterRelativeVorticity(ShallowWaterVorticity): u"""Diagnostic field for relative vorticity, ζ=∇×u""" - name = "RelativeVorticity" + name = "ShallowWaterRelativeVorticity" def __init__(self, space=None, method='solve'): """ From 998c2ce397cbc88c65ab074572ce1f1d0670482c Mon Sep 17 00:00:00 2001 From: Witt-D Date: Mon, 19 Aug 2024 10:37:53 +0100 Subject: [PATCH 15/17] current state --- unit-tests/test_vorticity_diagnostic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unit-tests/test_vorticity_diagnostic.py b/unit-tests/test_vorticity_diagnostic.py index 9855701f7..158799c88 100644 --- a/unit-tests/test_vorticity_diagnostic.py +++ b/unit-tests/test_vorticity_diagnostic.py @@ -20,7 +20,7 @@ def test_relative_vorticity(topology): def sphere_test(): # try at a higher resolution R = 1 # radius of ball H = 5 # height of model top - nlayers=5 + nlayers=8 c=16 # Building model and state object m = GeneralCubedSphereMesh(R, num_cells_per_edge_of_panel=c, degree=2) From 12667141c40a7d9c776928051a118e2b42da57b8 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Tue, 17 Sep 2024 15:22:05 +0100 Subject: [PATCH 16/17] changed sign --- gusto/diagnostics/compressible_euler_diagnostics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gusto/diagnostics/compressible_euler_diagnostics.py b/gusto/diagnostics/compressible_euler_diagnostics.py index 62194672f..f77f395f3 100644 --- a/gusto/diagnostics/compressible_euler_diagnostics.py +++ b/gusto/diagnostics/compressible_euler_diagnostics.py @@ -851,7 +851,7 @@ def setup(self, domain, state_fields, vorticity_type=None): vort = TrialFunction(space) gamma = TestFunction(space) a = inner(vort, gamma) * dx - L = -(inner(domain.perp(grad(gamma)), u)) * dx + L = (inner(domain.perp(grad(gamma)), u)) * dx # TODO implement absolute version, unsure how to get corioilis in vertical slice smartly problem = LinearVariationalProblem(a, L, self.field) From b024780dc3ef5ebfc5c04a55faf0575d06fcadd7 Mon Sep 17 00:00:00 2001 From: Witt-D Date: Thu, 19 Sep 2024 17:03:50 +0100 Subject: [PATCH 17/17] added jump term --- .../diagnostics/compressible_euler_diagnostics.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/gusto/diagnostics/compressible_euler_diagnostics.py b/gusto/diagnostics/compressible_euler_diagnostics.py index f77f395f3..9f80ed4c2 100644 --- a/gusto/diagnostics/compressible_euler_diagnostics.py +++ b/gusto/diagnostics/compressible_euler_diagnostics.py @@ -1,9 +1,9 @@ """Common diagnostic fields for the compressible Euler equations.""" -from firedrake import (dot, dx, Function, ln, TestFunction, TrialFunction, - Constant, grad, inner, LinearVariationalProblem, +from firedrake import (dot, dx, Function, ln, TestFunction, TrialFunction, + Constant, grad, inner, LinearVariationalProblem, dS, LinearVariationalSolver, FacetNormal, ds_b, dS_v, div, - avg, jump, SpatialCoordinate, curl, as_vector) + avg, jump, SpatialCoordinate, curl, as_vector, cross) from gusto.diagnostics.diagnostics import ( DiagnosticField, Energy, IterativeDiagnosticField @@ -842,16 +842,19 @@ def setup(self, domain, state_fields, vorticity_type=None): if domain.mesh.topological_dimension() == 3: vort = TrialFunction(space) gamma = TestFunction(space) + f = Coefficient(space) + breakpoint() + n = FacetNormal(domain.mesh) a = inner(vort, gamma) * dx - L = -inner(curl(gamma), u) * dx + L = inner(curl(gamma), u) * dx if vorticity_type == 'absolute': Omega = as_vector((0, 0, self.parameters.Omega)) - L += inner(2*Omega, gamma) * dx + L += inner(2*Omega, gamma) * dx - jump(inner(cross(u, n), gamma), n)*dS else: vort = TrialFunction(space) gamma = TestFunction(space) a = inner(vort, gamma) * dx - L = (inner(domain.perp(grad(gamma)), u)) * dx + L = (inner(domain.perp(grad(gamma)), u)) * dx - jump(inner(cross(u, n), gamma, n))*dS # TODO implement absolute version, unsure how to get corioilis in vertical slice smartly problem = LinearVariationalProblem(a, L, self.field)