Vector Reconstruction
MPAS models carry horizontal vector fields as a single scalar component
normal to each edge (normalVelocity, for example). Analysis, plotting and
coupling all need the full vector at cell centers or vertices instead.
polaris.mesh.reconstruct builds the least-squares weights that turn
edge-normal values into a Cartesian vector at those points, and applies them.
The method is the linear least-squares reconstruction of Peixoto and Barros (2014) on a two-ring edge stencil (their Figure 5). See the design document on vector reconstruction for the derivation and for why this method was chosen over Perot’s; this page covers only what a caller of the framework needs to know.
Weights are a property of the mesh
The weights depend on mesh geometry alone, not on the field being
reconstructed, so they are computed once per mesh and stored as mesh fields.
polaris.mesh.reconstruct.compute_reconstruction_weights() returns a
small dataset with just those fields, and
polaris.mesh.reconstruct.add_reconstruction_weights_to_dataset()
merges them into a mesh dataset:
from polaris.mesh.reconstruct import add_reconstruction_weights_to_dataset
ds_mesh = add_reconstruction_weights_to_dataset(ds_mesh, location='cell')
Both take a location of either 'cell' or 'vertex' and write three
fields, named after that location:
location |
stencil |
edge count |
weights |
|---|---|---|---|
|
|
|
|
|
|
|
|
The stencil holds 1-based edge indices, padded with zeros, and the edge count
says how many of them are valid. The weights have an extra R3 dimension
because the reconstructed vector is Cartesian.
compute_reconstruction_weights() selects only the coordinate and
connectivity variables it needs before loading anything, so it is safe to hand
it a production mesh with large state and tracer fields attached; those are
never read. It works with both MPAS-Ocean and Omega naming conventions and
writes output in whichever convention the input used. The standalone driver
add_reconstruction_weights wraps it for adding weights to an
existing mesh file from the command line.
The stencil is two rings of edges
For a cell, the stencil is the union of edgesOnVertex over the
verticesOnCell of that cell, with duplicates removed — 12 unique edges on a
hexagon and 10 on a pentagon, so it fits the existing maxEdges2 dimension.
For a vertex, it is the union of edgesOnVertex over the one-ring of
neighboring vertices, which is 9 edges on a triangular dual mesh and uses a
NINE dimension.
This is a wider stencil than the one-ring edgesOnCell that the legacy
coeffs_reconstruct field assumes, which is why the stencil is stored
explicitly rather than left implicit in the mesh connectivity.
Some meshes (QU240km among them) mark invalid connectivity entries with
nIndices + 1 rather than zero. compute_reconstruction_weights() normalizes
those to zero before building the stencil, so callers do not need to.
Applying the weights
polaris.mesh.reconstruct.tangential_reconstruction() contracts the
weights against an edge-normal field and returns the three Cartesian
components. Pass the precomputed stencil and weights when the mesh already
carries them; otherwise pass location and they are built on the fly, which is
much slower and worth avoiding when reconstructing many time slices:
from polaris.mesh.reconstruct import tangential_reconstruction
u_x, u_y, u_z = tangential_reconstruction(
ds_mesh,
ds.normalVelocity,
stencil=ds_mesh.reconstructStencilCell,
weights=ds_mesh.reconstructWeightsCell,
)
The edge-normal field may carry any extra dimensions besides nEdges —
nVertLevels, Time or both. They are broadcast through unchanged, so a
whole time series of a 3D field is reconstructed in one call.
On a sphere the Cartesian result is usually not what you want to plot.
polaris.mesh.reconstruct.cartesian_to_local_geographic() rotates it
into zonal, meridional and radial components at the reconstruction point:
from polaris.mesh.reconstruct import cartesian_to_local_geographic
u_zonal, u_merid, u_radial = cartesian_to_local_geographic(
ds_mesh, u_x, u_y, u_z
)
It infers whether the values are cell- or vertex-centered from the dimensions
of u_x, so there is no location argument. Because the reconstruction is
tangential, u_radial is zero up to truncation error.
On a planar mesh, do not apply this rotation. The tangent plane is the mesh
plane itself, so u_x and u_y are already the in-plane components and u_z
is zero.
Planar meshes
The same functions work on planar meshes, but two things differ internally, and both matter if you are reading the code:
the rotation to the local tangent plane is the identity, since every reconstruction point already lies in the mesh plane; and
because the identity rotation leaves the edge coordinates absolute rather than relative to the reconstruction point, they are translated explicitly. Without that, the least-squares fit would be anchored at the mesh origin and the constant term would be the field extrapolated there instead of the value at the reconstruction point.
Whether a mesh is planar is decided with Mesh Information, from the
on_a_sphere attribute, so a dataset missing that attribute raises a
ValueError rather than being reconstructed as the wrong kind of mesh.
Validating against a baseline
polaris.mesh.reconstruct.get_reconstruction_validate_vars() returns
the names of the three fields for a location, in the form
polaris.Step.add_output_file() expects:
self.add_output_file(
filename='reconstruction_weights.nc',
validate_vars=get_reconstruction_validate_vars(location='cell'),
)
See Baseline Validation Variables for the equivalent lists for mesh and cell-width files.
Internals
The remaining functions in the module — stencil construction, the rotation
matrix, the tangent-plane projection, the Renka (1984) least-squares weights
and the pseudo-inverse solve — are steps of the single pipeline that
compute_reconstruction_weights() drives. They are not part of the API this
page describes, and their signatures may change; read
polaris/mesh/reconstruct.py if you need them.