Skip to content

New method for calculating vorticity diagnostic #734

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jshipton merged 4 commits into from
Jun 3, 2026
Merged

New method for calculating vorticity diagnostic #734

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
e119325
an attempt at new vorticity diagnostic
tommbendall Apr 30, 2026
596f582
implement new method for calculating vorticity diagnostic
tommbendall May 1, 2026
72d1559
fix vorticity space logic
tommbendall May 19, 2026
b5328e9
Merge branch 'main' into TBendall/NewVorticity
tommbendall May 21, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
77 changes: 62 additions & 15 deletions gusto/diagnostics/shallow_water_diagnostics.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,11 @@

from firedrake import (
dx, TestFunction, TrialFunction, grad, inner, curl, Function, assemble,
LinearVariationalProblem, LinearVariationalSolver, conditional, interpolate
LinearVariationalProblem, LinearVariationalSolver, conditional, interpolate,
Projector
)
from gusto.diagnostics.diagnostics import DiagnosticField, Energy
from gusto.core.function_spaces import is_cg

__all__ = ["ShallowWaterKineticEnergy", "ShallowWaterPotentialEnergy",
"ShallowWaterPotentialEnstrophy", "PotentialVorticity",
Expand Down Expand Up @@ -157,28 +159,49 @@ def setup(self, domain, state_fields, vorticity_type=None):
vorticity_types = ["relative", "absolute", "potential"]
if vorticity_type not in vorticity_types:
raise ValueError(f"vorticity type must be one of {vorticity_types}, not {vorticity_type}")
space = domain.spaces("H1")

u = state_fields("u")
# Set default space
if self.space is None:
self.space = domain.spaces("H1")

# Work out continuity requirements of space
if not is_cg(self.space):
# Using DG, calculate u in HCurl space as intermediary
self.hcurl_intermediary = True
assert self.method != 'solve', (

@jshipton jshipton May 5, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This assertion is causing test failures... I'm not quite sure what you intend here...

@tommbendall tommbendall May 5, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Interesting... I wonder if that actually might explain a real problem with the vorticity diagnosics!

'vorticity diagnostics can only be output in L2 space when '
+ 'not using "solve" method'
)
u_hdiv = state_fields("u")
self.u_hcurl = Function(domain.spaces("HCurl"))
u = self.u_hcurl
curl_operator = domain.divperp

else:
# Using CG, calculate vorticity directly from HDiv u
self.hcurl_intermediary = False
u = state_fields("u")
curl_operator = curl

if vorticity_type in ["absolute", "potential"]:
f = state_fields("coriolis")
if vorticity_type == "potential":
D = state_fields("D")

if self.method != 'solve':
if vorticity_type == "potential":
self.expr = (curl(u) + f) / D
self.expr = (curl_operator(u) + f) / D
elif vorticity_type == "absolute":
self.expr = curl(u) + f
self.expr = curl_operator(u) + f
elif vorticity_type == "relative":
self.expr = curl(u)
self.expr = curl_operator(u)

super().setup(domain, state_fields, space=space)
super().setup(domain, state_fields, space=self.space)

# Set up problem now that self.field has been set up
if self.method == 'solve':
gamma = TestFunction(space)
q = TrialFunction(space)
gamma = TestFunction(self.space)
q = TrialFunction(self.space)

if vorticity_type == "potential":
a = q*gamma*D*dx
Expand All @@ -196,6 +219,15 @@ def setup(self, domain, state_fields, vorticity_type=None):
constant_jacobian=constant_jacobian)
self.evaluator = LinearVariationalSolver(problem, solver_parameters={"ksp_type": "cg"})

elif self.hcurl_intermediary:
self.hcurl_projection = Projector(u_hdiv, self.u_hcurl)

def compute(self):
"""Compute the vorticity, possibly via velocity in HCurl space."""
if self.hcurl_intermediary:
self.hcurl_projection.project()
super().compute()


class PotentialVorticity(Vorticity):
u"""Diagnostic field for shallow-water potential vorticity, q=(∇×(u+f))/D"""
Expand All @@ -206,10 +238,15 @@ def __init__(self, space=None, method='solve'):
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.
case a default space will be chosen for this diagnostic, which
is the H1 space. Note that if the diagnostic is in the L2 space,
it should be evaluated with the 'interpolate' or 'project'
method, and will project the velocity into the HCurl space as
an intermediary step.
method (str, optional): a string specifying the method of evaluation
for this diagnostic. Valid options are 'interpolate', 'project',
'assign' and 'solve'. Defaults to 'solve'.
'assign' and 'solve'. Defaults to 'solve'. Note that the 'solve'
method is only available when the diagnostic is in the H1 space.
"""
self.solve_implemented = True
super().__init__(space=space, method=method,
Expand All @@ -235,10 +272,15 @@ def __init__(self, space=None, method='solve'):
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.
case a default space will be chosen for this diagnostic, which
is the H1 space. Note that if the diagnostic is in the L2 space,
it should be evaluated with the 'interpolate' or 'project'
method, and will project the velocity into the HCurl space as
an intermediary step.
method (str, optional): a string specifying the method of evaluation
for this diagnostic. Valid options are 'interpolate', 'project',
'assign' and 'solve'. Defaults to 'solve'.
'assign' and 'solve'. Defaults to 'solve'. Note that the 'solve'
method is only available when the diagnostic is in the H1 space.
"""
self.solve_implemented = True
super().__init__(space=space, method=method, required_fields=('u', 'coriolis'))
Expand All @@ -263,10 +305,15 @@ def __init__(self, space=None, method='solve'):
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.
case a default space will be chosen for this diagnostic, which
is the H1 space. Note that if the diagnostic is in the L2 space,
it should be evaluated with the 'interpolate' or 'project'
method, and will project the velocity into the HCurl space as
an intermediary step.
method (str, optional): a string specifying the method of evaluation
for this diagnostic. Valid options are 'interpolate', 'project',
'assign' and 'solve'. Defaults to 'solve'.
'assign' and 'solve'. Defaults to 'solve'. Note that the 'solve'
method is only available when the diagnostic is in the H1 space.
"""
self.solve_implemented = True
super().__init__(space=space, method=method, required_fields=('u',))
Expand Down
Loading