Inconsistency in vertices and edges

[changliao1025]

I ran a mesh generation simulation successfully. But a script crashed when I tried to convert the invert_mesh.nc to a different format. I debugged and found something unexpected:

It appears cell 60469 has six vertices but only five edges.

[changliao1025]

At the same time, the nEdgeonCell shows seven:

[xylar]

@changliao1025, the indexing is 1-based (Fortran) indexing. Are you sure that is what you are using? I think such fundamental issues would have broken our workflows long ago if they were caused by MPAS-Tools.

[xylar]

If you'd like us to take a closer look, you're going to need to point to the mesh somewhere so we can have a look.

[xylar]

Here is the workflow that @changliao1025 shared with me:

def _land_mask_from_geojson( mesh_filename=sFilename_base_mesh,
                                geojson_filename=sFilename_land_coverage,
                                mask_filename=sFilename_land_mask):
    dsLandMask = xr.open_dataset(sFilename_land_mask)

    sFilename_culled_mesh = os.path.join(sWorkspace_jigsaw_out, "out", "culled_graph.info")

    dsCulledMesh = cull(
        dsBaseMesh, dsMask=dsLandMask,
        graphInfoFileName=sFilename_culled_mesh)

    sFilename_invert_mesh = os.path.join(sWorkspace_jigsaw_out, "out", "invert_graph.info")
    dsInvertMesh = cull(
        dsBaseMesh, dsInverse=dsLandMask,
        graphInfoFileName= sFilename_invert_mesh)

    sFilename_culled_mesh = os.path.join(sWorkspace_jigsaw_out, "out", "culled_mesh.nc")
    write_netcdf(
        dsCulledMesh, sFilename_culled_mesh, netcdfFormat)

    sFilename_invert_mesh = os.path.join(sWorkspace_jigsaw_out, "out", "invert_mesh.nc")
    write_netcdf(
        dsInvertMesh, sFilename_invert_mesh, netcdfFormat)

Here is my suggested fix, calling convert():

def _land_mask_from_geojson( mesh_filename=sFilename_base_mesh,
                                geojson_filename=sFilename_land_coverage,
                                mask_filename=sFilename_land_mask):
    dsLandMask = xr.open_dataset(sFilename_land_mask)

    sFilename_culled_mesh = os.path.join(sWorkspace_jigsaw_out, "out", "culled_graph.info")

    dsCulledMesh = cull(
        dsBaseMesh, dsMask=dsLandMask)
    dsCulledMesh = convert(
        dsCulledMesh, graphInfoFileName=sFilename_culled_mesh)

    sFilename_invert_mesh = os.path.join(sWorkspace_jigsaw_out, "out", "invert_graph.info")
    dsInvertMesh = cull(
        dsBaseMesh, dsInverse=dsLandMask)
    dsInvertMesh = convert(
        dsInvertMesh, graphInfoFileName=sFilename_invert_mesh)

    sFilename_culled_mesh = os.path.join(sWorkspace_jigsaw_out, "out", "culled_mesh.nc")
    write_netcdf(
        dsCulledMesh, sFilename_culled_mesh, netcdfFormat)

    sFilename_invert_mesh = os.path.join(sWorkspace_jigsaw_out, "out", "invert_mesh.nc")
    write_netcdf(
        dsInvertMesh, sFilename_invert_mesh, netcdfFormat)

[xylar]

I was able to visualize the original dataset (invert_mesh.nc) on edges and vertices but not on cells, so that suggests that there may be a problem with the connectivity array verticesOnCell, which is used to create cell visualizations. This is a very large mesh so not super easy to debug. If the convert() calls work, I'm probably going to have to let that be the solution.

[changliao1025]

There is indeed an issue with some code. I used the new code but got a different error caused by this:
In verticesOnCell in the invert_mesh file: