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

Generalizing WMT (both 3D and 2D) with restructured classes; no redundant methods/functions #40

Open
wants to merge 91 commits into
base: dev_wmt3d
Choose a base branch
from

Conversation

hdrake
Copy link
Collaborator

@hdrake hdrake commented May 3, 2023

This PR picks up from #35 and aims to simultaneously integrate swmt.py into wmt.py while also addressing various structural limitations discussed in the following issues:

The commits containing the most important major changes are:

Remaining to do list:

As of xgcm/xgcm#470, xgcm.Grid.interp no longer supports
the `boundary='extrapolate'` option. It is not clear that this was justified
anyways, so I replaced it with the safer `boundary='extend'` option. This works
with xgcm version 0.8.1 (up from 0.5.2 in https://github.com/NOAA-GFDL/xwmt/blob/797cd8827e91128e1c842c1d149bc0dc1cd84e88/notebooks/xwmt_full_3D_WMT_tutorial.ipynb).
- Also reduced the amount of hard-coded logic, and clarified the terminology, for `_group_process` -> `_group_processes`, which now includes advection (vertical + horizontal) and diabatic (external forcing + diffusion), and `_group_tend` -> `_sum_components`.

- Changed the convention that deleted the components that go into these groupings -- it is useful to have all of the different components accessible by default.
- WaterMass and WaterMassTransformation now require information about the
xgcm.Grid coordinates and metrics
- I realized that plotting low-resolution water mass transformation rates as line plots,
and comparing them to high-resolution analytical solution curves, was misleading. The
numerical water mass transformation terms represent averaged transformation rates within
a finite volume, while the analytical solutions are point-wise estimates. To clarify this
distinction, I am not plotting the watermass transformation rates uses "matplotlib.pyplot.step".
This demonstrates that even at low resolution, the water mass transformation rates are not actually erroneous!

- I also added some more non-monotic stratification examples to probe the behavior of the singularities
when the magnitude of the gradient of lambda tends towards zero. As expected, watermass transformations
blow up to ±infinity in this limit.
- Main takeaway is that there are problems when there are lambda bins
that do not correspond to any cells in the dataset.
Regardless of whether trying to compute a single term (specified by a string)
or all of the terms in the budget, return a xr.Dataset. Partially addresses the
last point in NOAA-GFDL#8
…onverge

- Considers four different idealized experiments, with various vertical profiles
of lambda and extensive lambda tendencies. Includes two cases with vanishing
gradients (water mass transformations that diverge to ±ifinity) and one case
in which z(lambda) has multiple values.
- A first test checks that, with very high vertical resolution, the numerical results
converge are close (rtol < 1.e-5) to the layer-integrated analytical solution.
- A second test checks that, as vertical resolution (in z or in lambda bins) is increased,
the numerical results converge on the point-wise analytical solution in lambda-space.

All tests passed locally.
This completes the separation of the WaterMass class (consisting of the ds, grid,
and variables/constants/methods related to thermodynamics and density) from the
WaterMassTransformation class. The names of the temperature and salinity variables,
required for density calculations, are now keyword arguments to WaterMass.__init__,
rather than inferred from the watermass transformation variable dictionaries required
for WaterMassTransformations.__init__.

Also changed the structure such that all derived thermodynamic variables (e.g. alpha,
beta, p, sa, ct) are now just added to WaterMass.ds, rather than as top-level class
attributes.

- Begins to address NOAA-GFDL#36
- Adresses all of NOAA-GFDL#8, except the part
about expanding 2D fields to 3D. I have not yet brought swmt.py up to date with
the changes herein.
- Removed `swmt.py` and the `SurfaceWaterMassTransformation` class
within it, with the idea that all of its functionality is already
included within the core of `wmt.py` (and now `wm.py`).

- Moved the core (parent) `WaterMass` class to `wm.py`. This was
actually done in the last two commits, but here we include the
associated updates to the exporting in `xwmt.__init__.py`. The
keyword argument `t_name='temperature'` is now passed to `WaterMass`
in the analytical convergence tests, which had broken.

- Generalized core `WaterMassTransformations.datadict` method to
now identify surface fluxes as tendencies that do not have any
vertical coordinates (taken from `self.grid.axes['Z'].coords`)
in their dimensions. Addresses

- The vertical coordinate used in both calculations for pressure
(in the `wm.WaterMass.get_density` method) and for the 2D->3D
expansion that is currently used to compute WMT from surface fluxes
(see `compute.expand_surface_to_3d`) has been renamed to `z` (with
levels `z_l` and interfaces `z_i`) to emphasize that these only
support depth coordinates as currently implemented. For WMT calculations
in arbitrary vertical coordinates, these will need to be fundamentally
changed. (See related discussion at NOAA-GFDL#3).

- Added `/conventions/MOM6.yaml` file that can be read in as a
python dictionary that is formatted as required by the
`budgets_dict` argument of `WaterMassTransformations`, and
includes the default names of the core heat, salt, and mass
tendencies required for full WMT calculations in temperature,
salinity, and density space.

Issues addressed by this commit (or previous ones building up to)
- Makes many of the re-labellings requested in NOAA-GFDL#28

- Partially addresses NOAA-GFDL#13, although the information
is still stored as a dictionary. It should not be fairly straight forward to encode this information
in attributes of the data arrays in `self.ds`.

- Completely addresses NOAA-GFDL#7 and NOAA-GFDL#8,
I think, except for the expansion of 2D surface fluxes (see above).
There are now no redundant function or method definitions, as we
have removed `swmt.py` completely and made use of class inheritance.
@hdrake hdrake marked this pull request as draft May 3, 2023 15:13
@hdrake
Copy link
Collaborator Author

hdrake commented May 3, 2023

@gmacgilchrist, do you have any suggestions for how to replace the xwmt.compute.expand_surface_to_3d function with a more robust approach (that would work for any vertical coordinate)?

hdrake added 3 commits May 3, 2023 15:17
In the previous commit, I was not carefully treating the case of
2D surface flux terms, which should be fed through `hlamdot_from_Jlam`.
I now expand them to 3D, closely following the relevant lines in `swmt.py`:
- https://github.com/NOAA-GFDL/xwmt/blob/07b2e6fc95fd5bb715b461656509a8497c21b87a/xwmt/swmt.py#L122-L125
- https://github.com/NOAA-GFDL/xwmt/blob/07b2e6fc95fd5bb715b461656509a8497c21b87a/xwmt/swmt.py#L194-L214
- https://github.com/NOAA-GFDL/xwmt/blob/07b2e6fc95fd5bb715b461656509a8497c21b87a/xwmt/compute.py#L163-L178

I confirmed the 2D surface flux terms now give the same total WMT
as the 3D extensive "boundary_forcing" term.
Integrating 2D boundary fluxes into `xwmt`'s Water Mass Transformation
framework requires placing the correspoinding flux convergence at
the correct vertical level of the model (or, at least, into the correct
lambda bin). Finding the correct lambda bin is an equivalent problem to
finding the outcrop level (for surface fluxes) or incrop level (for
seafloor fluxes, e.g. geothermal or benthic tracer fluxes).

While finding outcrops is trivial in `zstr` coordinates (top level),
it is more challenging in generalized vertical coordinates such as
a `rho2` diagnostic coordinate or the native `hybrid` coordinate.
Flipping the kwarg incrop to True returns incrops instead, which are
non-trivial in all coordinates.

Note that these methods assume the coordinate are monotonically
increasing, with larger values being "down" in physical space
(which is currently true of all supported diagnostic vertical
coordinates, but may not always be true).
@hdrake
Copy link
Collaborator Author

hdrake commented May 4, 2023

Commit 32c6b44 addresses #3 (comment)

Major changes:
- Manually set vertical thickness metrics within a new
"Z_metrics" attribute in `WaterMass.grid` because `xgcm`
methods do not currently correctly handle the case where
the metrics have larger dimensions than the dimension they
describe the metric for, e.g. vertical cell thicknesses
(`thkcello`) which describe thickness in the `Z` dimension
but also vary over ('X','Y') and "time". In particular,
the `grid.get_metrics` method was having a hard time picking
out the correct thickness metrics for vertical derivatives or
integrals.

- `WaterMass.get_density` now always uses cumulative integrals of
thickness (interpolated onto interface levels) to define the center
depth of each layer. Previously, this function only worked for `zstr`
vertical coordinate because it assumed `z_l` was in the dataset.

- `expand_surface_to_3d` has been renamed to `expand_surface_array_vertically`
and refactored from being a function in `compute.py` to a method on `wm.WaterMass`.
It now correctly identifies the outcrop level in all cases (assuming
a monotonic diagnostic coordinate that increases from the surface down)
and masks all over levels. Previously, it only correctly worked for `zstr` because it
masked everything except for the top level.

- Added docstrings to key classes and methods (but have not yet updated
all of the functions in `compute.py` or the internal methods in
`wmt.py`).

Minor changes:
- Made `get_outcrop_lev` and `sel_outcrop_lev` functions more robust.

- Moved `zonal_mean` and `bin_percentile` from functions in
`compute.py` to methods on `WaterMass` in `wm.py`. I've also
reduced the number of hard-coded grid variables that they relied
on, such as `z_l`, `areacello`, `wet`.

- Incremented (significantly) the pinned `xgcm` version in
`ci/environment.yml` and removed `xhistogram` dependency.
First, if the dataset already contains lambda as coordinate, then
only rebin if `rebin=True` kwarg is passed. If it does rebin, then
it uses xgcm to do conservative one-dimensional rebinning to target bins.
Otherwise, either use 3D lamda field for 3D rebinning, with the option
for calculating density from temperature and salinity offline if needed.

- Rework organization of under-the-hood transformation functions.

- Changed default of `bin_percentile` to [0., 1.] (i.e. min and max)
- Work directly with `self.grid._ds` instead of `self.ds`
This prevents awkward errors when trying to feed the same dataset or grid object
into multiple different WaterMass or WaterMassBudget objects. For example, lazy
variables for convergent transport sections will have different sizes of the
`sect` dimension for differently shaped regions. These will conflict instead of being
overridden by new variables.
hdrake and others added 23 commits December 19, 2023 18:33
Previously, we were using in-situ thermal expansion and haline contraction
coefficients to convert temperature and salinity tendencies into potential density
tendencies. Upon further thought, we believe this to be incorrect, as a more consistent
version of the evolution equation for potential density (referenced to p_ref) would use
EOS coefficients that also assume constant pressure (i.e. are also referenced to p_ref).
Fixed inconsistency in installation instructions noticed by @cspencerjones
Previously, we assumed that budget tendencies provided to `xwmt` were per-unit-area.
This makes sense for MOM6 because this is how budget terms are output, but is a poor
choice for finite-volume models in general, which deal in discrete integrals.

Requires `xbudget >= v0.0.3`
Assume the `t_name` and `s_name` variables provided as arguments
to the `xwmt.WaterMass` instance constructor refer to Conservative
Temperature (CT) and Absolute Salinity (SA) by default. Optional
key word argument can be switched to "potential" to instead take in
potential temperature and practical salinity.

Also, estimate layer depths by integrating thickness downwards by default,
instead of only doing this when calling `xwmt.WaterMass.get_density`.
We have changed the sign of WMTs such that they are positive when they
increase the mass of water with tracer values less than the threshold.
This is a breaking change, but earlier results can be recovered by
simply flipping the sign of the results. The package now agrees
with the sign convention in the original Walin 1982 paper.

Some other minor changes:
-Bumped the minor version number because of quasi-breaking change in
sign of transformation rates
-Add condition so that method still works for 2D (surface flux) diagnostics
-Rename LHS and RHS terms from "material_derivative" to "transformation"

Update all notebooks accordingly
Reverse sign of water mass transformations
@hdrake hdrake mentioned this pull request Sep 23, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
5 participants