Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
2 changes: 2 additions & 0 deletions conda_package/mpas_tools/io.py
Original file line number Diff line number Diff line change
Expand Up @@ -135,6 +135,8 @@ def write_netcdf(
# make sure the Time dimension is unlimited because MPAS has trouble
# reading Time otherwise
ds.encoding['unlimited_dims'] = {'Time'}
else:
ds.encoding['unlimited_dims'] = None

# for performance, we have to handle this as a special case
convert = format == 'NETCDF3_64BIT_DATA'
Expand Down
105 changes: 58 additions & 47 deletions conda_package/mpas_tools/mesh/creation/sort_mesh.py
Original file line number Diff line number Diff line change
@@ -1,11 +1,13 @@
import argparse
import os

import numpy as np
import os
import xarray
import argparse
from scipy.sparse import csr_matrix
from scipy.sparse.csgraph import reverse_cuthill_mckee

from mpas_tools.io import write_netcdf


def sort_mesh(mesh):
"""
Expand All @@ -24,9 +26,9 @@ def sort_mesh(mesh):
"""
# Authors: Darren Engwirda

ncells = mesh.sizes["nCells"]
nedges = mesh.sizes["nEdges"]
nduals = mesh.sizes["nVertices"]
ncells = mesh.sizes['nCells']
nedges = mesh.sizes['nEdges']
nduals = mesh.sizes['nVertices']

cell_fwd = np.arange(0, ncells) + 1
cell_rev = np.arange(0, ncells) + 1
Expand All @@ -42,23 +44,20 @@ def sort_mesh(mesh):
cell_rev = np.zeros(ncells, dtype=np.int32)
cell_rev[cell_fwd - 1] = np.arange(ncells) + 1

mesh["cellsOnCell"][:] = \
_sort_rev(mesh["cellsOnCell"], cell_rev)
mesh["cellsOnEdge"][:] = \
_sort_rev(mesh["cellsOnEdge"], cell_rev)
mesh["cellsOnVertex"][:] = \
_sort_rev(mesh["cellsOnVertex"], cell_rev)
mesh['cellsOnCell'][:] = _sort_rev(mesh['cellsOnCell'], cell_rev)
mesh['cellsOnEdge'][:] = _sort_rev(mesh['cellsOnEdge'], cell_rev)
mesh['cellsOnVertex'][:] = _sort_rev(mesh['cellsOnVertex'], cell_rev)

for var in mesh.keys():
dims = mesh.variables[var].dims
if ("nCells" in dims):
if 'nCells' in dims:
mesh[var][:] = _sort_fwd(mesh[var], cell_fwd)

mesh["indexToCellID"][:] = np.arange(ncells) + 1
mesh['indexToCellID'][:] = np.arange(ncells) + 1

# sort duals via pseudo-linear cell-wise ordering

dual_fwd = np.ravel(mesh["verticesOnCell"].values)
dual_fwd = np.ravel(mesh['verticesOnCell'].values)
dual_fwd = dual_fwd[dual_fwd > 0]

__, imap = np.unique(dual_fwd, return_index=True)
Expand All @@ -68,22 +67,20 @@ def sort_mesh(mesh):
dual_rev = np.zeros(nduals, dtype=np.int32)
dual_rev[dual_fwd - 1] = np.arange(nduals) + 1

mesh["verticesOnCell"][:] = \
_sort_rev(mesh["verticesOnCell"], dual_rev)
mesh['verticesOnCell'][:] = _sort_rev(mesh['verticesOnCell'], dual_rev)

mesh["verticesOnEdge"][:] = \
_sort_rev(mesh["verticesOnEdge"], dual_rev)
mesh['verticesOnEdge'][:] = _sort_rev(mesh['verticesOnEdge'], dual_rev)

for var in mesh.keys():
dims = mesh.variables[var].dims
if ("nVertices" in dims):
if 'nVertices' in dims:
mesh[var][:] = _sort_fwd(mesh[var], dual_fwd)

mesh["indexToVertexID"][:] = np.arange(nduals) + 1
mesh['indexToVertexID'][:] = np.arange(nduals) + 1

# sort edges via pseudo-linear cell-wise ordering

edge_fwd = np.ravel(mesh["edgesOnCell"].values)
edge_fwd = np.ravel(mesh['edgesOnCell'].values)
edge_fwd = edge_fwd[edge_fwd > 0]

__, imap = np.unique(edge_fwd, return_index=True)
Expand All @@ -93,60 +90,75 @@ def sort_mesh(mesh):
edge_rev = np.zeros(nedges, dtype=np.int32)
edge_rev[edge_fwd - 1] = np.arange(nedges) + 1

mesh["edgesOnCell"][:] = \
_sort_rev(mesh["edgesOnCell"], edge_rev)
mesh['edgesOnCell'][:] = _sort_rev(mesh['edgesOnCell'], edge_rev)

mesh["edgesOnEdge"][:] = \
_sort_rev(mesh["edgesOnEdge"], edge_rev)
mesh['edgesOnEdge'][:] = _sort_rev(mesh['edgesOnEdge'], edge_rev)

mesh["edgesOnVertex"][:] = \
_sort_rev(mesh["edgesOnVertex"], edge_rev)
mesh['edgesOnVertex'][:] = _sort_rev(mesh['edgesOnVertex'], edge_rev)

for var in mesh.keys():
dims = mesh.variables[var].dims
if ("nEdges" in dims):
if 'nEdges' in dims:
mesh[var][:] = _sort_fwd(mesh[var], edge_fwd)

mesh["indexToEdgeID"][:] = np.arange(nedges) + 1
mesh['indexToEdgeID'][:] = np.arange(nedges) + 1

return mesh


def main():
parser = argparse.ArgumentParser(
description=__doc__,
formatter_class=argparse.RawTextHelpFormatter)
description=__doc__, formatter_class=argparse.RawTextHelpFormatter
)

parser.add_argument(
'--mesh-file',
dest='mesh_file',
type=str,
required=True,
help='Path+name to unsorted mesh file.',
)

parser.add_argument(
"--mesh-file", dest="mesh_file", type=str,
required=True, help="Path+name to unsorted mesh file.")
'--sort-file',
dest='sort_file',
type=str,
required=True,
help='Path+name to sorted output file.',
)

parser.add_argument(
"--sort-file", dest="sort_file", type=str,
required=True, help="Path+name to sorted output file.")
'--format',
dest='format',
type=str,
required=False,
default='NETCDF4',
help='Output format for the sorted mesh file.',
)

args = parser.parse_args()

mesh = xarray.open_dataset(args.mesh_file)

sort_mesh(mesh)

with open(os.path.join(os.path.dirname(
args.sort_file), "graph.info"), "w") as fptr:
cellsOnCell = mesh["cellsOnCell"].values
with open(
os.path.join(os.path.dirname(args.sort_file), 'graph.info'), 'w'
) as fptr:
cellsOnCell = mesh['cellsOnCell'].values

ncells = mesh.sizes["nCells"]
ncells = mesh.sizes['nCells']
nedges = np.count_nonzero(cellsOnCell) // 2

fptr.write(f"{ncells} {nedges}\n")
fptr.write(f'{ncells} {nedges}\n')
for cell in range(ncells):
data = cellsOnCell[cell, :]
data = data[data > 0]
for item in data:
fptr.write(f"{item} ")
fptr.write("\n")
fptr.write(f'{item} ')
fptr.write('\n')

mesh.to_netcdf(args.sort_file, format="NETCDF4")
write_netcdf(mesh, args.sort_file, format=args.format)


def _sort_fwd(data, fwd):
Expand Down Expand Up @@ -212,11 +224,10 @@ def _cell_del2(mesh):
ivec = np.array([], dtype=np.int32)
jvec = np.array([], dtype=np.int32)

topolOnCell = mesh["nEdgesOnCell"].values
cellsOnCell = mesh["cellsOnCell"].values
topolOnCell = mesh['nEdgesOnCell'].values
cellsOnCell = mesh['cellsOnCell'].values

for edge in range(np.max(topolOnCell)):

# cell-to-cell pairs, if edges exist
mask = topolOnCell > edge
idx_self = np.argwhere(mask).ravel()
Expand All @@ -237,5 +248,5 @@ def _cell_del2(mesh):
return csr_matrix((xvec, (ivec, jvec)))


if (__name__ == "__main__"):\
if __name__ == '__main__':
main()
Loading