diff --git a/conda_package/mpas_tools/io.py b/conda_package/mpas_tools/io.py index f4244b4ab..b6ebec839 100644 --- a/conda_package/mpas_tools/io.py +++ b/conda_package/mpas_tools/io.py @@ -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' diff --git a/conda_package/mpas_tools/mesh/creation/sort_mesh.py b/conda_package/mpas_tools/mesh/creation/sort_mesh.py index 9b7ec3b9f..3729d3f72 100644 --- a/conda_package/mpas_tools/mesh/creation/sort_mesh.py +++ b/conda_package/mpas_tools/mesh/creation/sort_mesh.py @@ -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): """ @@ -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 @@ -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) @@ -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) @@ -93,37 +90,51 @@ 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() @@ -131,22 +142,23 @@ def main(): 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): @@ -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() @@ -237,5 +248,5 @@ def _cell_del2(mesh): return csr_matrix((xvec, (ivec, jvec))) -if (__name__ == "__main__"):\ +if __name__ == '__main__': main()