MPAS to Omegah mesh converter

CMakeLists.txt

@@ -31,6 +31,7 @@ bob_option(Omega_h_USE_CUDA_AWARE_MPI "Assume MPI is CUDA-aware, make use of tha
 bob_option(Omega_h_USE_SimModSuite "Enable reading Simmetrix MeshSim meshes" OFF)
 bob_input(Omega_h_VALGRIND "" STRING "Valgrind plus arguments for testing")
 bob_option(Omega_h_EXAMPLES "Compile examples" OFF)
+bob_option(Omega_h_USE_MPAS "Support reading MPAS meshes" OFF)
 
 bob_option(Omega_h_USE_MPI "Use MPI for parallelism" OFF)
 bob_option(Omega_h_ENABLE_DEMANGLED_STACKTRACE "Add linker options to enable human-readable stacktraces on gnu-c platforms." OFF)
@@ -141,6 +142,12 @@ if (Omega_h_USE_CUDA)
   endif()
 endif()
 
+
+if (Omega_h_USE_MPAS)
+  find_package(PkgConfig REQUIRED)
+  pkg_check_modules(NetCDF REQUIRED IMPORTED_TARGET netcdf>=4.7.3)
+endif()
+
 include(cmake/osh_use_dolfin.cmake)
 osh_use_dolfin()
 
@@ -161,6 +168,7 @@ set(Omega_h_KEY_BOOLS
     Omega_h_USE_ZLIB
     Omega_h_USE_libMeshb
     Omega_h_USE_EGADS
+    Omega_h_USE_MPAS
     Omega_h_USE_SEACASExodus
     Omega_h_USE_SimModSuite
     Omega_h_USE_DOLFIN

src/CMakeLists.txt

@@ -43,6 +43,7 @@ set(Omega_h_SOURCES
   Omega_h_ghost.cpp
   Omega_h_globals.cpp
   Omega_h_gmsh.cpp
+  Omega_h_mpas.cpp
   Omega_h_grammar.cpp
   Omega_h_graph.cpp
   Omega_h_hilbert.cpp
@@ -215,6 +216,10 @@ if(Omega_h_USE_EGADS)
   target_link_libraries(omega_h PUBLIC "${EGADS_LIBRARY}")
 endif()
 
+if(Omega_h_USE_MPAS)
+  target_link_libraries(omega_h PUBLIC PkgConfig::NetCDF)
+endif()
+
 if(Omega_h_USE_SimModSuite)
   target_include_directories(omega_h PUBLIC "${SIMMODSUITE_INCLUDE_DIR}")
   target_link_libraries(omega_h PUBLIC "${SIMMODSUITE_LIBS}")
@@ -272,6 +277,9 @@ macro(osh_add_util EXE_NAME)
   bob_export_target(${EXE_NAME})
 endmacro(osh_add_util)
 
+if(Omega_h_USE_MPAS)
+  osh_add_util(mpas2osh)
+endif()
 osh_add_util(msh2osh)
 osh_add_util(osh2vtk)
 osh_add_util(oshdiff)
@@ -374,6 +382,15 @@ if(BUILD_TESTING)
   set(TEST_EXES ${TEST_EXES} unit_math)
   test_basefunc(run_unit_math 1 ./unit_math)
 
+  if (Omega_h_USE_MPAS)
+    if (Omega_h_MPAS_MESH_DIR)
+      test_func(run_mpas2osh 1 ./mpas2osh
+        ${Omega_h_MPAS_MESH_DIR}/greenland.nc
+        ${Omega_h_MPAS_MESH_DIR}/fieldList.txt
+        ${Omega_h_MPAS_MESH_DIR}/greenland.osh)
+    endif()
+  endif()
+
   if(Omega_h_USE_SimModSuite)
     set(d3dMesh ${CMAKE_SOURCE_DIR}/meshes/d3d)
     test_func(run_tri_convert 1 ./meshsim2osh

src/Omega_h_file.cpp

@@ -698,6 +698,8 @@ OMEGA_H_DLL Mesh read_mesh_file(filesystem::path const& path, CommPtr comm) {
 #endif
   } else if (extension == ".msh") {
     return gmsh::read(path, comm);
+  } else if (extension == ".nc") {
+    return mpas::read(path, comm);
   } else if (extension == ".pvtu") {
     Mesh mesh(comm->library());
     vtk::read_parallel(path, comm, &mesh);

src/Omega_h_file.hpp

@@ -106,6 +106,12 @@ void write_parallel(filesystem::path const& filename, Mesh& mesh);
 
 }  // namespace gmsh
 
+namespace mpas {
+Mesh read(filesystem::path const& filename, const std::vector<std::string>& vtxFieldList, CommPtr comm);
+Mesh read(filesystem::path const& filename, filesystem::path const& vtxFieldListFile, CommPtr comm);
+Mesh read(filesystem::path const& filename, CommPtr comm);
+}
+
 namespace vtk {
 static constexpr bool do_compress = true;
 static constexpr bool dont_compress = false;

src/Omega_h_mpas.cpp

@@ -0,0 +1,220 @@
+#include "Omega_h_file.hpp"
+
+#include <algorithm>
+#include <fstream>
+#include <iomanip>
+#include <cctype>
+#include <sstream>
+#include <unordered_map>
+
+#include "Omega_h_build.hpp"
+#include "Omega_h_class.hpp"
+#include "Omega_h_element.hpp"
+#include "Omega_h_map.hpp"
+#include "Omega_h_vector.hpp"
+
+#ifdef OMEGA_H_USE_MPAS
+#include <netcdf.h>
+#endif  // OMEGA_H_USE_MPAS
+
+namespace {
+
+enum FieldType { integer, real };
+
+typedef std::map<FieldType, std::string> FieldMap;
+
+//TODO read from file
+FieldMap readVtxFieldList(Omega_h::filesystem::path const& filename) {
+  return {
+    {real, "observedSurfaceVelocityX"},
+    {real, "observedSurfaceVelocityY"},
+    {real, "observedSurfaceVelocityUncertainty"},
+    {real, "sfcMassBal"},
+    {real, "surfaceAirTemperature"},
+    {real, "basalHeatFlux"},
+    {real, "observedThicknessTendency"},
+    {real, "floatingBasalMassBal"},
+    {integer, "idCell"}
+  };
+}
+
+} //end anonymous namespace
+
+namespace Omega_h {
+
+namespace mpas {
+
+#ifndef OMEGA_H_USE_MPAS
+
+Mesh read(filesystem::path const&, const FieldMap&, CommPtr) {
+  Omega_h_fail("recompile with Omega_h_USE_MPAS enabled!\n");
+  return Mesh(comm->library()); //silence warning
+}
+
+Mesh read(filesystem::path const&, filesystem::path const&, CommPtr) {
+  Omega_h_fail("recompile with Omega_h_USE_MPAS enabled!\n");
+  return Mesh(comm->library()); //silence warning
+}
+
+Mesh read(filesystem::path const&, filesystem::path const&, CommPtr comm) {
+  Omega_h_fail("recompile with Omega_h_USE_MPAS enabled!\n");
+  return Mesh(comm->library()); //silence warning
+}
+
+#else
+
+inline void checkErr(int err) {
+  if(err != NC_NOERR) {
+    printf("Error: %d %s\n", err, nc_strerror(err));
+  }
+  OMEGA_H_CHECK(err==NC_NOERR);
+}
+
+size_t readDimension(int ncid, const std::string name) {
+  int dimId;
+  int err = nc_inq_dimid(ncid, name.c_str(), &dimId);
+  checkErr(err);
+  size_t dim;
+  err = nc_inq_dimlen(ncid, dimId, &dim);
+  checkErr(err);
+  return dim;
+}
+
+template<class T>
+HostWrite<T> readArray(int ncid, const std::string name, const size_t len) {
+  int arrId;
+  int err = nc_inq_varid(ncid, name.c_str(), &arrId);
+  checkErr(err);
+  auto arr = HostWrite<T>(len);
+  err = nc_get_var(ncid, arrId, arr.data());
+  return arr;
+}
+
+HostWrite<Real> createCoordinates(int ncid, bool useCartesianCoords,
+    const size_t dim, const size_t nCells) {
+  OMEGA_H_CHECK(dim==2); //TODO support 3d
+  auto x = useCartesianCoords ? readArray<Real>(ncid, "xCell", nCells) :
+                                      readArray<Real>(ncid, "latCell", nCells);
+  auto y = useCartesianCoords ? readArray<Real>(ncid, "yCell", nCells) :
+                                      readArray<Real>(ncid, "lonCell", nCells);
+  HostWrite<Real> coords(nCells*dim);
+  for(size_t i=0; i<nCells; i++) {
+    coords[i*2]   = x[i];
+    coords[i*2+1] = y[i];
+  }
+  return coords;
+}
+
+template <class T>
+bool anyZeros(T& arr, const size_t start, const size_t end) {
+  for(size_t i=start; i<end; i++) {
+    if(arr[i] <= 0) return true;
+  }
+  return false;
+}
+
+template <class T>
+void copySubArray(T& src, T& dest,
+    const size_t srcStart, const size_t destStart, const size_t len) {
+  for(size_t i=0; i<len; i++) {
+    dest[destStart+i] = src[srcStart+i];
+  }
+}
+
+HostWrite<LO> createElemConnectivity(HostWrite<LO>& mpasConn, const size_t nPrimalVtx) {
+  OMEGA_H_CHECK(nPrimalVtx*3 == mpasConn.size());
+  //Primal vertices that don't bound three primal cells
+  // don't form valid dual triangles.
+  //This handles boundaries of the mesh.
+  //See figure 5.3 of MPAS v1.0 spec
+  int nTri=0;
+  for(int i=0; i<nPrimalVtx; i++) {
+    if( ! anyZeros(mpasConn, i*3, (i+1)*3) )
+      nTri++;
+  }
+  printf("nTri %d\n", nTri);
+  //create omegah connectivity array
+  int triIdx=0;
+  HostWrite<LO> elm2vtx(nTri*3);
+  for(int i=0; i<nPrimalVtx; i++) {
+    if( ! anyZeros(mpasConn, i*3, (i+1)*3) ) {
+      copySubArray(mpasConn, elm2vtx, i*3, triIdx*3, 3);
+      triIdx++;
+    }
+  }
+  //decrement the indices - mpas uses 1-based indexing
+  for(int i=0; i<nTri*3; i++) {
+    --elm2vtx[i];
+    if(elm2vtx[i] < 0) printf("foo %d %d\n", i, elm2vtx[i]);
+    OMEGA_H_CHECK(elm2vtx[i]>=0);
+  }
+  return elm2vtx;
+}
+
+void read_internal(int ncid, bool useCartesianCoords,
+    const FieldMap& vtxFieldNames, Mesh* mesh) {
+  const size_t dim = 2;
+  const size_t nCells = readDimension(ncid, "nCells");
+  const size_t nVertices = readDimension(ncid, "nVertices");
+  auto coords = createCoordinates(ncid, useCartesianCoords, dim, nCells);
+  printf("primal cells verts: %u %u\n", nCells, nVertices);
+  auto cellsOnVertices = readArray<LO>(ncid, "cellsOnVertex", nVertices*3);
+  std::map<std::string, HostWrite<Real> > fieldValsReal;
+  std::map<std::string, HostWrite<LO> > fieldValsInt;
+  for(auto& [type,name] : vtxFieldNames) { //primal (polygonal) cell = dual (triangle) vertex
+    if(type == real) {
+      fieldValsReal[name] = readArray<Real>(ncid, name, nCells);
+    } else if(type == integer) {
+      fieldValsInt[name] = readArray<LO>(ncid, name, nCells);
+    } else {
+      Omega_h_fail("Unrecognized field type!\n");
+    }
+  }
+  auto elm2vtx = createElemConnectivity(cellsOnVertices, nVertices);
+  Read<LO> elm2vtx_d(elm2vtx.write());
+  Read<Real> coords_d(coords.write());
+  build_from_elems_and_coords(mesh, OMEGA_H_SIMPLEX, dim, elm2vtx_d, coords_d);
+  for(auto& [name,vals] : fieldValsReal)
+    mesh->add_tag<Real>(OMEGA_H_VERT, name, 1, vals.write());
+  for(auto& [name,vals] : fieldValsInt)
+    mesh->add_tag<LO>(OMEGA_H_VERT, name, 1, vals.write());
+}
+
+Mesh read(int ncid, const FieldMap& vtxFieldList, CommPtr comm) {
+  auto mesh = Mesh(comm->library());
+  bool useCartesianCoords = true;
+  if (comm->rank() == 0) {
+    read_internal(ncid, useCartesianCoords, vtxFieldList, &mesh);
+  }
+  fprintf(stderr, "done read\n");
+  mesh.set_comm(comm);
+  mesh.balance();
+  return mesh;
+}
+
+Mesh read(filesystem::path const& filename, const FieldMap& vtxFieldList, CommPtr comm) {
+  int retval;
+  int ncid;
+  retval = nc_open(filename.c_str(), NC_NOWRITE, &ncid);
+  checkErr(retval);
+  auto mesh = mpas::read(ncid, vtxFieldList, comm);
+  retval = nc_close(ncid);
+  checkErr(retval);
+  return mesh;
+}
+
+Mesh read(filesystem::path const& filename, filesystem::path const& vtxFieldListFile, CommPtr comm) {
+  const auto vtxFieldList = readVtxFieldList(vtxFieldListFile);
+  return mpas::read(filename, vtxFieldList, comm);
+}
+
+Mesh read(filesystem::path const& filename, CommPtr comm) {
+  const FieldMap vtxFieldList;
+  return mpas::read(filename, vtxFieldList, comm);
+}
+
+#endif  // OMEGA_H_USE_MPAS
+
+}  // namespace mpas
+
+}  // end namespace Omega_h

src/mpas2osh.cpp

@@ -0,0 +1,17 @@
+#include <Omega_h_file.hpp>
+#include <Omega_h_cmdline.hpp>
+
+int main(int argc, char** argv) {
+  auto lib = Omega_h::Library(&argc, &argv);
+  auto world = lib.world();
+  Omega_h::CmdLine cmdline;
+  cmdline.add_arg<std::string>("input.nc");
+  cmdline.add_arg<std::string>("vtxFieldList.txt");
+  cmdline.add_arg<std::string>("output.osh");
+  if (!cmdline.parse_final(world, &argc, argv)) return -1;
+  auto inpath = cmdline.get<std::string>("input.nc");
+  auto inVtxFieldList = cmdline.get<std::string>("vtxFieldList.txt");
+  auto outpath = cmdline.get<std::string>("output.osh");
+  auto mesh = Omega_h::mpas::read(inpath, inVtxFieldList, world);
+  Omega_h::binary::write(outpath, &mesh);
+}