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.

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

Improve tracer density #492

Merged
merged 12 commits into from
Mar 28, 2024
Merged

Improve tracer density #492

Show file tree
Hide file tree
Changes from 7 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
f9cad21
start better tracer_density diagnostic
ta440 Mar 22, 2024
7a7a099
start with more accurate space for tracer density
ta440 Mar 22, 2024
a411171
create space for tracer density diagnostic with TensorProductElement
ta440 Mar 25, 2024
baa1b66
better space used for tracer_density
ta440 Mar 26, 2024
c320322
Merge branch 'main' into improve_tracer_density
ta440 Mar 26, 2024
d0a12bc
lint
ta440 Mar 26, 2024
e1a5028
Merge branch 'main' into improve_tracer_density
tommbendall Mar 27, 2024
3a66b4e
allow user defined space in tracer density diagnostic
ta440 Mar 27, 2024
1cb388f
Merge branch 'improve_tracer_density' of github.com:firedrakeproject/…
ta440 Mar 27, 2024
0aebae3
add more comments for the tracer density diagnostic
ta440 Mar 27, 2024
3a9c38d
lint
ta440 Mar 27, 2024
5dad4ca
Merge branch 'main' into improve_tracer_density
tommbendall Mar 28, 2024
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
28 changes: 23 additions & 5 deletions gusto/active_tracers.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,8 @@ class WaterVapour(ActiveTracer):
"""An object encoding the details of water vapour as a tracer."""
def __init__(self, name='water_vapour', space='theta',
variable_type=TracerVariableType.mixing_ratio,
transport_eqn=TransportEquationType.advective):
transport_eqn=TransportEquationType.advective,
density_name=None):
"""
Args:
name (str, optional): the variable's name. Defaults to
Expand All @@ -113,16 +114,22 @@ def __init__(self, name='water_vapour', space='theta',
transport_eqn (:class:`TransportEquationType`, optional): enumerator
indicating the type of transport equation to be solved (e.g.
advective). Defaults to `TransportEquationType.advective`.
density_name (str): the name of the associated density for a mixing
ratio when using the tracer_conservative transport. Defaults to None,
as this argument is not needed for other transport types.
"""
super().__init__(f'{name}', space, variable_type,
transport_eqn=transport_eqn, phase=Phases.gas, chemical='H2O')
transport_eqn=transport_eqn,
density_name=density_name,
phase=Phases.gas, chemical='H2O')


class CloudWater(ActiveTracer):
"""An object encoding the details of cloud water as a tracer."""
def __init__(self, name='cloud_water', space='theta',
variable_type=TracerVariableType.mixing_ratio,
transport_eqn=TransportEquationType.advective):
transport_eqn=TransportEquationType.advective,
density_name=None):
"""
Args:
name (str, optional): the variable name. Default is 'cloud_water'.
Expand All @@ -134,15 +141,21 @@ def __init__(self, name='cloud_water', space='theta',
transport_eqn (:class:`TransportEquationType`, optional): enumerator
indicating the type of transport equation to be solved (e.g.
advective). Defaults to `TransportEquationType.advective`.
density_name (str): the name of the associated density for a mixing
ratio when using the tracer_conservative transport. Defaults to None,
as this argument is not needed for other transport types.
"""
super().__init__(f'{name}', space, variable_type,
transport_eqn=transport_eqn, phase=Phases.liquid, chemical='H2O')
transport_eqn=transport_eqn,
density_name=density_name,
phase=Phases.liquid, chemical='H2O')


class Rain(ActiveTracer):
"""An object encoding the details of rain as a tracer."""
def __init__(self, name='rain', space='theta',
variable_type=TracerVariableType.mixing_ratio,
density_name=None,
Copy link
Contributor

Choose a reason for hiding this comment

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

should the density_name argument come after the transport_eqn argument for consistency?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Oops, good spot, thanks!

transport_eqn=TransportEquationType.advective):
"""
Args:
Expand All @@ -155,6 +168,11 @@ def __init__(self, name='rain', space='theta',
transport_eqn (:class:`TransportEquationType`, optional): enumerator
indicating the type of transport equation to be solved (e.g.
advective). Defaults to `TransportEquationType.advective`.
density_name (str): the name of the associated density for a mixing
ratio when using the tracer_conservative transport. Defaults to None,
as this argument is not needed for other transport types.
"""
super().__init__(f'{name}', space, variable_type,
transport_eqn=transport_eqn, phase=Phases.liquid, chemical='H2O')
transport_eqn=transport_eqn,
density_name=density_name,
phase=Phases.liquid, chemical='H2O')
41 changes: 38 additions & 3 deletions gusto/diagnostics.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,8 @@
LinearVariationalProblem, LinearVariationalSolver, FacetNormal, \
ds_b, ds_v, ds_t, dS_h, dS_v, ds, dS, div, avg, jump, pi, \
TensorFunctionSpace, SpatialCoordinate, as_vector, \
Projector, Interpolator, FunctionSpace
Projector, Interpolator, FunctionSpace, FiniteElement, \
TensorProductElement
from firedrake.assign import Assigner

from abc import ABCMeta, abstractmethod, abstractproperty
Expand Down Expand Up @@ -1691,7 +1692,10 @@ class TracerDensity(DiagnosticField):
"""Diagnostic for computing the density of a tracer. This is
computed as the product of a mixing ratio and dry density"""

name = "TracerDensity"
@property
def name(self):
"""Gives the name of this diagnostic field."""
return "TracerDensity_"+self.mixing_ratio_name+'_'+self.density_name

def __init__(self, mixing_ratio_name, density_name, space=None, method='interpolate'):
"""
Expand All @@ -1705,7 +1709,9 @@ def __init__(self, mixing_ratio_name, density_name, space=None, method='interpol
for this diagnostic. Valid options are 'interpolate', 'project' and
'assign'. Defaults to 'interpolate'.
"""

super().__init__(method=method, required_fields=(mixing_ratio_name, density_name))

self.mixing_ratio_name = mixing_ratio_name
self.density_name = density_name

Expand All @@ -1717,7 +1723,36 @@ def setup(self, domain, state_fields):
domain (:class:`Domain`): the model's domain object.
state_fields (:class:`StateFields`): the model's field container.
"""

m_X = state_fields(self.mixing_ratio_name)
rho_d = state_fields(self.density_name)

m_X_space = m_X.function_space()
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 we want all of this code to be protected behind an if self.space_defined is None -- currently this diagnostic allows the user to pass in a space (and we want the user to be able to do that) but it will be overridden by this code. I think we only want to override if space is None

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Cool, I've changed it to enable this, although I've used if self.space is None instead of self.space_defined

rho_d_space = rho_d.function_space()

if domain.spaces.extruded_mesh:
m_X_horiz = m_X_space.ufl_element().sub_elements[0]
m_X_vert = m_X_space.ufl_element().sub_elements[1]
rho_d_horiz = rho_d_space.ufl_element().sub_elements[0]
rho_d_vert = rho_d_space.ufl_element().sub_elements[1]

horiz_degree = m_X_horiz.degree() + rho_d_horiz.degree()
vert_degree = m_X_vert.degree() + rho_d_vert.degree()

cell = domain.mesh._base_mesh.ufl_cell().cellname()
horiz_elt = FiniteElement('DG', cell, horiz_degree)
vert_elt = FiniteElement('DG', cell, vert_degree)
elt = TensorProductElement(horiz_elt, vert_elt)

else:
m_X_degree = m_X_space.ufl_element().degree()
rho_d_degree = rho_d_space.ufl_element().degree()
degree = m_X_degree + rho_d_degree

cell = domain.mesh.ufl_cell().cellname()
elt = FiniteElement('DG', cell, degree)

tracer_density_space = FunctionSpace(domain.mesh, elt, name='tracer_density_space')

self.expr = m_X*rho_d
super().setup(domain, state_fields)
super().setup(domain, state_fields, space=tracer_density_space)
Loading