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

Face Area Fix, Face Dimension Variable, & Data Type Fix #283

Merged
merged 7 commits into from
Apr 19, 2023
Merged
Show file tree
Hide file tree
Changes from 6 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
1 change: 1 addition & 0 deletions docs/internal_api/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@ Grid Methods
grid.Grid.__init_grid_var_attrs__
grid.Grid._populate_cartesian_xyz_coord
grid.Grid._populate_lonlat_coord
grid.Grid._build_face_dimension


Grid Helper Modules
Expand Down
1 change: 1 addition & 0 deletions test/constants.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,3 +11,4 @@
MESH30_AREA = 12.566
PSI_INTG = 12.566
VAR2_INTG = 12.566
UNIT_SPHERE_AREA = 4 * np.pi
59 changes: 47 additions & 12 deletions test/test_grid.py
Original file line number Diff line number Diff line change
Expand Up @@ -296,6 +296,20 @@ def test_read_scrip(self):
xr_grid_u30 = xr.open_dataset(ug_30)
ux_grid_u30 = ux.Grid(xr_grid_u30) # tests from ugrid

def test_build_face_dimension(self):
"""Tests the construction of the ``Mesh2_face_dimension`` variable."""
grids = [self.tgrid1, self.tgrid2, self.tgrid3]

for grid in grids:
# highest possible dimension dimension for a face
max_dimension = grid.nMaxMesh2_face_nodes

# face must be at least a triangle
min_dimension = 3

assert grid.Mesh2_face_dimension.min() >= min_dimension
assert grid.Mesh2_face_dimension.max() <= max_dimension


class TestIntegrate(TestCase):
mesh_file30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug"
Expand All @@ -307,16 +321,16 @@ def test_calculate_total_face_area_triangle(self):
verts = [[[0.57735027, -5.77350269e-01, -0.57735027],
[0.57735027, 5.77350269e-01, -0.57735027],
[-0.57735027, 5.77350269e-01, -0.57735027]]]
vgrid = ux.Grid(verts)
vgrid = ux.Grid(verts, vertices=True, islatlon=False, concave=False)

# get node names for each grid object
x_var = vgrid.ds_var_names["Mesh2_node_x"]
y_var = vgrid.ds_var_names["Mesh2_node_y"]
z_var = vgrid.ds_var_names["Mesh2_node_z"]

vgrid.ds[x_var].attrs["units"] = "m"
vgrid.ds[y_var].attrs["units"] = "m"
vgrid.ds[z_var].attrs["units"] = "m"
vgrid.Mesh2_node_x.attrs["units"] = "m"
vgrid.Mesh2_node_y.attrs["units"] = "m"
vgrid.Mesh2_node_z.attrs["units"] = "m"

area_gaussian = vgrid.calculate_total_face_area(
quadrature_rule="gaussian", order=5)
Expand All @@ -336,6 +350,26 @@ def test_calculate_total_face_area_file(self):

nt.assert_almost_equal(area, constants.MESH30_AREA, decimal=3)

def test_calculate_total_face_area_sphere(self):
"""Computes the total face area of an MPAS mesh that lies on a unit
sphere, with an expected total face area of 4pi."""
mpas_grid_path = current_path / 'meshfiles' / "mpas" / "QU" / 'mesh.QU.1920km.151026.nc'

ds = xr.open_dataset(mpas_grid_path)
primal_grid = ux.Grid(ds, use_dual=False)
dual_grid = ux.Grid(ds, use_dual=True)

primal_face_area = primal_grid.calculate_total_face_area()
dual_face_area = dual_grid.calculate_total_face_area()

nt.assert_almost_equal(primal_face_area,
constants.UNIT_SPHERE_AREA,
decimal=3)

nt.assert_almost_equal(dual_face_area,
constants.UNIT_SPHERE_AREA,
decimal=3)

erogluorhan marked this conversation as resolved.
Show resolved Hide resolved
def test_integrate(self):
xr_grid = xr.open_dataset(self.mesh_file30)
xr_psi = xr.open_dataset(self.data_file30)
Expand All @@ -359,14 +393,15 @@ def test_compute_face_areas_geoflow_small(self):
grid_1 = ux.Grid(grid_1_ds)
grid_1.compute_face_areas()

def test_compute_face_areas_fesom(self):
"""Checks if the FESOM PI-Grid Output can generate a face areas
output."""

fesom_grid_small = current_path / "meshfiles" / "ugrid" / "fesom" / "fesom.mesh.diag.nc"
grid_2_ds = xr.open_dataset(fesom_grid_small)
grid_2 = ux.Grid(grid_2_ds)
grid_2.compute_face_areas()
# removed test until fix to tranposed face nodes
# def test_compute_face_areas_fesom(self):
# """Checks if the FESOM PI-Grid Output can generate a face areas
# output."""
#
# fesom_grid_small = current_path / "meshfiles" / "ugrid" / "fesom" / "fesom.mesh.diag.nc"
# grid_2_ds = xr.open_dataset(fesom_grid_small)
# grid_2 = ux.Grid(grid_2_ds)
# grid_2.compute_face_areas()
rajeeja marked this conversation as resolved.
Show resolved Hide resolved


class TestPopulateCoordinates(TestCase):
Expand Down
15 changes: 15 additions & 0 deletions test/test_helpers.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,10 +34,12 @@ def test_face_area_coords(self):
y = np.array([-5.77350269e-01, 5.77350269e-01, 5.77350269e-01])
z = np.array([-0.57735027, -0.57735027, -0.57735027])
face_nodes = np.array([[0, 1, 2]]).astype(INT_DTYPE)
face_dimension = np.array([3], dtype=INT_DTYPE)
area = ux.get_all_face_area_from_coords(x,
y,
z,
face_nodes,
face_dimension,
3,
coords_type="cartesian")
nt.assert_almost_equal(area, constants.TRI_AREA, decimal=1)
Expand Down Expand Up @@ -131,6 +133,19 @@ class TestConstants(TestCase):
# INT_FILL_VALUE as set in constants.py
fv = INT_FILL_VALUE

def test_invalid_indexing(self):
"""Tests if the current INT_DTYPE and INT_FILL_VALUE throw the correct
errors when indexing."""
dummy_data = np.array([1, 2, 3, 4])

invalid_indices = np.array([self.fv, self.fv], dtype=INT_DTYPE)
invalid_index = self.fv

# invalid index/indices should throw an Index Error
with self.assertRaises(IndexError):
dummy_data[invalid_indices]
dummy_data[invalid_index]

def test_replace_fill_values(self):
"""Tests _replace_fill_values() helper function across multiple
different dtype arrays used as face_nodes."""
Expand Down
4 changes: 2 additions & 2 deletions test/test_mpas.py
Original file line number Diff line number Diff line change
Expand Up @@ -98,8 +98,8 @@ def test_add_fill_values(self):
assert np.array_equal(verticesOnCell, gold_output)

def test_set_attrs(self):
"""Tests the execution of ``_set_global_attrs``, checking for
attributes being correctly stored in ``Grid.ds``"""
"""Tests the execution of "_set_global_attrs", checking for attributes
being correctly stored in "Grid.ds"."""

# full set of expected mpas attributes
expected_attrs = [
Expand Down
5 changes: 3 additions & 2 deletions uxarray/constants.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
import numpy as np

INT_DTYPE = np.uintp
INT_FILL_VALUE = np.iinfo(INT_DTYPE).max
# numpy indexing code is written for np.intp
INT_DTYPE = np.intp
INT_FILL_VALUE = np.iinfo(INT_DTYPE).min
33 changes: 32 additions & 1 deletion uxarray/grid.py
Original file line number Diff line number Diff line change
Expand Up @@ -99,6 +99,9 @@ def __init__(self, dataset, **kwargs):
# initialize convenience attributes
self.__init_grid_var_attrs__()

# build face dimension, possibly safeguard for large datasets
self._build_face_dimension()
philipc2 marked this conversation as resolved.
Show resolved Hide resolved

def __init_ds_var_names__(self):
"""Populates a dictionary for storing uxarray's internal representation
of xarray object.
Expand Down Expand Up @@ -352,6 +355,7 @@ def compute_face_areas(self, quadrature_rule="triangular", order=4):
coords_type = "cartesian"

face_nodes = self.Mesh2_face_nodes.data
face_dimension = self.Mesh2_face_dimension.data
dim = self.Mesh2.attrs['topology_dimension']

# initialize z
Expand All @@ -366,7 +370,8 @@ def compute_face_areas(self, quadrature_rule="triangular", order=4):

# call function to get area of all the faces as a np array
self._face_areas = get_all_face_area_from_coords(
x, y, z, face_nodes, dim, quadrature_rule, order, coords_type)
x, y, z, face_nodes, face_dimension, dim, quadrature_rule,
order, coords_type)

return self._face_areas

Expand Down Expand Up @@ -554,3 +559,29 @@ def _populate_lonlat_coord(self):
"long_name": "latitude of mesh nodes",
"units": "degrees_north",
})

def _build_face_dimension(self):
"""Constructs ``Mesh2_face_dimension``, which calculates the dimension
of each face in ``Mesh2_face_nodes``"""

# Triangular Mesh
if not hasattr(self, "nMaxMesh2_face_nodes"):
nMaxMesh2_face_nodes = self.Mesh2_face_nodes.shape[1]
setattr(self, "nMaxMesh2_face_nodes", nMaxMesh2_face_nodes)

# padding to shape [nMesh2_face, nMaxMesh2_face_nodes + 1]
closed = np.ones((self.nMesh2_face, self.nMaxMesh2_face_nodes + 1),
dtype=INT_DTYPE) * INT_FILL_VALUE

closed[:, :-1] = self.Mesh2_face_nodes.copy()

face_dimension = np.argmax(closed == INT_FILL_VALUE, axis=1)

# add to internal dataset
self.ds["Mesh2_face_dimension"] = xr.DataArray(
data=face_dimension,
dims=["nMesh2_face"],
attrs={"long_name": "number of non-fill value nodes for each face"})

# standardized attribute
setattr(self, "Mesh2_face_dimension", self.ds["Mesh2_face_dimension"])
24 changes: 14 additions & 10 deletions uxarray/helpers.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
from numba import njit, config
import math

from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
from .constants import INT_DTYPE

config.DISABLE_JIT = False
Expand Down Expand Up @@ -188,6 +189,7 @@ def get_all_face_area_from_coords(x,
y,
z,
face_nodes,
face_geometry,
dim,
quadrature_rule="triangular",
order=4,
Expand Down Expand Up @@ -226,25 +228,27 @@ def get_all_face_area_from_coords(x,
-------
area of all faces : ndarray
"""
num_faces = face_nodes.shape[0]
area = np.zeros(num_faces) # set area of each face to 0

face_nodes = face_nodes[:].astype(INT_DTYPE)
n_face, n_max_face_nodes = face_nodes.shape

for i in range(num_faces):
face_z = np.zeros(len(face_nodes[i]))
# set initial area of each face to 0
area = np.zeros(n_face)

for face_idx, max_nodes in enumerate(face_geometry):
face_x = x[face_nodes[face_idx, 0:max_nodes]]
face_y = y[face_nodes[face_idx, 0:max_nodes]]

face_x = x[face_nodes[i]]
face_y = y[face_nodes[i]]
# check if z dimension
if dim > 2:
face_z = z[face_nodes[i]]
face_z = z[face_nodes[face_idx, 0:max_nodes]]
else:
face_z = np.zeros_like(face_y)

# After getting all the nodes of a face assembled call the cal. face area routine
face_area = calculate_face_area(face_x, face_y, face_z, quadrature_rule,
order, coords_type)

area[i] = face_area
# store current face area
area[face_idx] = face_area

return area

Expand Down