8000 Add `meshDensity` in mesh conversion if not already present by xylar · Pull Request #583 · MPAS-Dev/MPAS-Tools · GitHub
[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Skip to content

Add meshDensity in mesh conversion if not already present #583

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.

Sign up for GitHub

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

Jump to bottom
Merged
merged 1 commit into from
Sep 26, 2024
Merged

Add meshDensity in mesh conversion if not already present #583

merged 1 commit into from
Sep 26, 2024

Conversation

xylar
Copy link
Collaborator
@xylar xylar commented Sep 26, 2024

@changliao1025 discovered that the new MpasMeshConverter.x (the one with NetCDF-C support) doesn't correctly handle the case where meshDensity is not passed in. In the older version of the tool, meshDensity was set to all ones if it isn't present, and I have added that here.

@xylar xylar added the bug label Sep 26, 2024
@xylar xylar self-assigned this Sep 26, 2024
@xylar
Copy link
Collaborator Author
xylar commented Sep 26, 2024

Testing

I was able to convert an MPAS mesh from JIGSAW without meshDensity:

$ ncdump -h mesh_in.nc 
netcdf mesh_in {
dimensions:
	nCells = 237355 ;
	nVertices = 474706 ;
	vertexDegree = 3 ;
variables:
	double xCell(nCells) ;
	double yCell(nCells) ;
	double zCell(nCells) ;
	int featureTagCell(nCells) ;
	double xVertex(nVertices) ;
	double yVertex(nVertices) ;
	double zVertex(nVertices) ;
	int featureTagVertex(nVertices) ;
	int cellsOnVertex(nVertices, vertexDegree) ;

// global attributes:
		:on_a_sphere = "YES" ;
		:sphere_radius = 6371000. ;
		:history = "Wed Sep 25 13:33:41 2024: /qfs/people/liao313/workspace/python/pyflowline_icom/examples/susquehanna/run_simulation_mpas.py" ;
}

to a full MPAS mesh:

$ MpasMeshConverter.x mesh_in.nc mesh_out.nc
$ ncdump -h mesh_out.nc 
netcdf mesh_out {
dimensions:
	nCells = 237355 ;
	nEdges = 712059 ;
	nVertices = 474706 ;
	maxEdges = 8 ;
	maxEdges2 = 16 ;
	TWO = 2 ;
	vertexDegree = 3 ;
	Time = UNLIMITED ; // (0 currently)
variables:
	double latCell(nCells) ;
		latCell:long_name = "latitudes of cell centres" ;
	double lonCell(nCells) ;
		lonCell:long_name = "longitudes of cell centres" ;
	double xCell(nCells) ;
		xCell:long_name = "x-coordinates of cell centres" ;
	double yCell(nCells) ;
		yCell:long_name = "y-coordinates of cell centres" ;
	double zCell(nCells) ;
		zCell:long_name = "z-coordinates of cell centres" ;
	int indexToCellID(nCells) ;
		indexToCellID:long_name = "index to cell ID mapping" ;
	double latEdge(nEdges) ;
		latEdge:long_name = "latitudes of edge centres" ;
	double lonEdge(nEdges) ;
		lonEdge:long_name = "longitudes of edge centres" ;
	double xEdge(nEdges) ;
		xEdge:long_name = "x-coordinates of edge centres" ;
	double yEdge(nEdges) ;
		yEdge:long_name = "y-coordinates of edge centres" ;
	double zEdge(nEdges) ;
		zEdge:long_name = "z-coordinates of edge centres" ;
	int indexToEdgeID(nEdges) ;
		indexToEdgeID:long_name = "index to edge ID mapping" ;
	double latVertex(nVertices) ;
		latVertex:long_name = "latitudes of vertices" ;
	double lonVertex(nVertices) ;
		lonVertex:long_name = "longitudes of vertices" ;
	double xVertex(nVertices) ;
		xVertex:long_name = "x-coordinates of vertices" ;
	double yVertex(nVertices) ;
		yVertex:long_name = "y-coordinates of vertices" ;
	double zVertex(nVertices) ;
		zVertex:long_name = "z-coordinates of vertices" ;
	int indexToVertexID(nVertices) ;
		indexToVertexID:long_name = "index to vertex ID mapping" ;
	int cellsOnCell(nCells, maxEdges) ;
		cellsOnCell:long_name = "cells adj. to each cell" ;
	int edgesOnCell(nCells, maxEdges) ;
		edgesOnCell:long_name = "edges on each cell" ;
	int verticesOnCell(nCells, maxEdges) ;
		verticesOnCell:long_name = "vertices on each cell" ;
	int nEdgesOnCell(nCells) ;
		nEdgesOnCell:long_name = "number of edges on each cell" ;
	int edgesOnEdge(nEdges, maxEdges2) ;
		edgesOnEdge:long_name = "edges adj. to each edge" ;
	int cellsOnEdge(nEdges, TWO) ;
		cellsOnEdge:long_name = "cells adj. to each edge" ;
	int verticesOnEdge(nEdges, TWO) ;
		verticesOnEdge:long_name = "vertices on each edge" ;
	int nEdgesOnEdge(nEdges) ;
		nEdgesOnEdge:long_name = "number of edges on each edge" ;
	int cellsOnVertex(nVertices, vertexDegree) ;
		cellsOnVertex:long_name = "vertices adj. to each vertex" ;
	int edgesOnVertex(nVertices, vertexDegree) ;
		edgesOnVertex:long_name = "edges adj. to each vertex" ;
	int boundaryVertex(nVertices) ;
		boundaryVertex:long_name = "non-zero for each vertex on mesh boundary" ;
	double areaCell(nCells) ;
		areaCell:long_name = "surface areas of cells" ;
	double angleEdge(nEdges) ;
		angleEdge:long_name = "angle to edges" ;
	double dcEdge(nEdges) ;
		dcEdge:long_name = "length of arc between centres" ;
	double dvEdge(nEdges) ;
		dvEdge:long_name = "length of arc between vertices" ;
	double weightsOnEdge(nEdges, maxEdges2) ;
		weightsOnEdge:long_name = "tangential flux reconstruction weights" ;
	double areaTriangle(nVertices) ;
		areaTriangle:long_name = "surface areas of dual cells" ;
	double kiteAreasOnVertex(nVertices, vertexDegree) ;
		kiteAreasOnVertex:long_name = "surface areas of overlap between cells and dual cells" ;
	double cellQuality(nCells) ;
		cellQuality:long_name = "quality of mesh cells" ;
	double gridSpacing(nCells) ;
		gridSpacing:long_name = "grid spacing distribution" ;
	double triangleQuality(nVertices) ;
		triangleQuality:long_name = "quality of mesh dual cells" ;
	double triangleAngleQuality(nVertices) ;
		triangleAngleQuality:long_name = "quality of mesh dual cells" ;
	int obtuseTriangle(nVertices) ;
		obtuseTriangle:long_name = "non-zero for any dual cell containing obtuse angles" ;
	double meshDensity(nCells) ;
		meshDensity:long_name = "mesh density distribution" ;

// global attributes:
		:on_a_sphere = "YES" ;
		:sphere_radius = 6371000. ;
		:is_periodic = "NO" ;
		:history = "MpasMeshConverter.x mesh_in.nc mesh_out.nc\nWed Sep 25 13:33:41 2024: /qfs/people/liao313/workspace/python/pyflowline_icom/examples/susquehanna/run_simulation_mpas.py" ;
		:mesh_spec = "1.0" ;
		:Conventions = "MPAS" ;
		:source = "MpasMeshConverter.x" ;
		:file_id = "42zjng5yki" ;
}

Without this fix, an error was thrown about meshDensity not being present:

terminate called after throwing an instance of 'std::invalid_argument'
  what():  Error getting variable meshDensity : -49
Aborted (core dumped)

@changliao1025
Copy link

It would be helpful to provide a link to the documentation explaining why meshDensity is necessary and how it can be generated.

@xylar
Copy link
Collaborator Author
xylar commented Sep 26, 2024

I don't believe meshDensity is necessary, expect that it is part of the MPAS mesh spec:
https://mpas-dev.github.io/files/documents/MPAS-MeshSpec.pdf
It is a relatively complicated function of resolution that I don't have handy. It used to be part of certain MPAS-Ocean paraemterizations but is no longer used for anything as far as I know. While we could probably just remove it from the converter, I think the safest thing to do for now is to create a field of all ones (as to tool did before) and let folks remove it again if they don't need it.

@xylar
Copy link
Collaborator Author
xylar commented Sep 26, 2024

Here is the definition of meshDensity:

MPAS-Tools/conda_package/mpas_tools/ocean/inject_meshDensity.py

Lines 54 to 59 in 2a46cbe

def inject_spherical_meshDensity(cellWidth, lon, lat, mesh_filename):
"""
Add a ``meshDensity`` field into a spherical MPAS mesh. The mesh density is
defined as:
meshDensity = (minCellWidth / cellWidth)**4

@xylar xylar merged commit 20b3ca4 into MPAS-Dev:master Sep 26, 2024
9 checks passed
@xylar xylar deleted the fix-mesh-density-in-converter branch September 26, 2024 17:37
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers
No reviews
Assignees

@xylar xylar

Labels
bug
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.
2 participants
@xylar @changliao1025
0