Skip to content
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.

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

Compressible vorticity #495

Open
wants to merge 22 commits into
base: main
Choose a base branch
from
Open

Compressible vorticity #495

wants to merge 22 commits into from

Conversation

Witt-D
Copy link
Collaborator

@Witt-D Witt-D commented Mar 27, 2024

Added compressible vorticity diagnostics so that a vector vorticity can be calculated for the compressible euler equations for both the absolute and relative versions of vorticity

# TODO: add linearisation and label for this
residual += subject(prognostic(
inner(w, cross(2*Omega, u))*dx, "u"), self.X)

Copy link
Contributor

Choose a reason for hiding this comment

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

we better not lose the comment reminding us to add the linearisation here

def __init__(self, parameters, space=None, method='solve'):
u"""
Args:
space (:class:`FunctionSpace`, optional): the function space to
Copy link
Contributor

Choose a reason for hiding this comment

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

please update docstring to include parameters argument

Copy link
Contributor

Choose a reason for hiding this comment

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

I think it's just CompressibleParameters - this should also work for the Boussinesq equations so you can remove any mention of Euler... in fact, it also works for incompressible...

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':
Copy link
Contributor

Choose a reason for hiding this comment

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

why not if vorticity_type == 'absolute' as above?

@@ -1687,6 +1687,110 @@ def setup(self, domain, state_fields):
super().setup(domain, state_fields, vorticity_type="relative")


class CompressibleVorticity(DiagnosticField):
u""" Base diagnostic Field for three dimensional Vorticity """
Copy link
Contributor

Choose a reason for hiding this comment

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

does this (should this?) also work in a vertical slice?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

This is a good question, which im not sure on so will have a test around and decide if it does / should it, i think it should work but may require a bit of thinking

Copy link
Contributor

Choose a reason for hiding this comment

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

Since you've put in effort to make this work in a vertical slice, you can change this comment!

Copy link
Contributor

@jshipton jshipton left a comment

Choose a reason for hiding this comment

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

Thanks @Witt-D, it'll be great to have this in! I've just left a couple of minor comments to address...

@tommbendall
Copy link
Contributor

I'm afraid that this now conflicts with changes that have just gone onto main -- could you fix these and then we can look at getting this on?

Copy link
Contributor

@tommbendall tommbendall left a comment

Choose a reason for hiding this comment

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

Sorry it's taken me a long time to look at this. I'm just leaving a few thoughts which hopefully shouldn't be too complicated to deal with!

I've looked at the failing tests, and I think this is just because there are a couple of places (examples/shallow_water/thermal_williamson2.py and integration_tests/equations/test_sw_triangle.py) with existing vorticity diagnostics that need the parameters argument adding

You also have some lint failures to fix!

self.expression = curl(u)
elif vorticity_type == 'absolute':
Omega = as_vector((0, 0, self.parameters.Omega))
self.expression = curl(u) + 2*Omega
Copy link
Contributor

Choose a reason for hiding this comment

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

Should this be self.expr?

u = state_fields('u')
if self.method != 'solve':
if vorticity_type == 'relative':
self.expression = curl(u)
Copy link
Contributor

Choose a reason for hiding this comment

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

Should this be self.expr?

be computed ('relative', 'absolute'). Defaults to
None.
"""
# TODO Do we eventually want a potential voriticy?
Copy link
Contributor

Choose a reason for hiding this comment

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

Yes I think we will -- you can remove this comment for now though

@@ -1687,6 +1687,110 @@ def setup(self, domain, state_fields):
super().setup(domain, state_fields, vorticity_type="relative")


class CompressibleVorticity(DiagnosticField):
u""" Base diagnostic Field for three dimensional Vorticity """
Copy link
Contributor

Choose a reason for hiding this comment

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

Since you've put in effort to make this work in a vertical slice, you can change this comment!

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
# TODO implement absolute version, unsure atm how to get corioilis in vertical slice smartly
Copy link
Contributor

Choose a reason for hiding this comment

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

This isn't obvious to me, and it isn't obvious that we'll ever want the absolute vorticity in a vertical slice. I suggest that you put an if check here so that we fail with a NotImplementedError in this situation

residual += coriolis(subject(prognostic(
if parameters.Omega is not None:
# TODO add linerisation and label for this
Omega = as_vector((0, 0, parameters.Omega))
Copy link
Contributor

Choose a reason for hiding this comment

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

Maybe we should call this variable Omega_vec to avoid confusion



class AbsoluteVorticity(Vorticity):
u""" Diagnostic field for compressible euler absolute vorticity """
Copy link
Contributor

Choose a reason for hiding this comment

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

As Jemma says, you could remove reference to compressible euler here?

@tommbendall tommbendall added enhancement Involves adding a new capability tidying Involves tidying up code labels Aug 7, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement Involves adding a new capability tidying Involves tidying up code
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants