Moab instance for mpas_ocean

Makefile

@@ -450,6 +450,13 @@ ifneq "$(NETCDF)" ""
 	LIBS += $(NCLIB)
 endif
 
+ifneq "$(MOAB_PATH)" ""
+	CPPINCLUDES +=  -DHAVE_MOAB -I$(MOAB_PATH)/include
+	FCINCLUDES += -DHAVE_MOAB -I$(MOAB_PATH)/include
+	include $(MOAB_PATH)/lib/moab.make
+	LIBS +=  ${MOAB_LIBS_LINK} -lstdc++
+endif
+
 RM = rm -f
 CPP = cpp -P -traditional
 RANLIB = ranlib
@@ -539,6 +546,12 @@ else # USE_PIO2 IF
 	PIO_MESSAGE="Using the PIO 1.x library."
 endif # USE_PIO2 IF
 
+ifneq "$(MOAB_PATH)" ""
+	MOAB_MESSAGE="Using MOAB library"
+else # 
+	MOAB_MESSAGE="Not using MOAB library"
+endif 
+
 ifdef TIMER_LIB
 ifeq "$(TIMER_LIB)" "tau"
 	override TAU=true
@@ -729,6 +742,7 @@ endif
 	@echo $(GEN_F90_MESSAGE)
 	@echo $(TIMER_MESSAGE)
 	@echo $(PIO_MESSAGE)
+	@echo $(MOAB_MESSAGE)
 	@echo "*******************************************************************************"
 clean:
 	cd src; $(MAKE) clean RM="$(RM)" CORE="$(CORE)"

src/Makefile.in.ACME

@@ -37,6 +37,16 @@ else
    ESMFDIR = noesmf
 endif
 
+# Set MOAB info if it is being used
+ifeq ($(strip $(USE_MOAB)), TRUE)
+  ifdef MOAB_PATH
+    CPPDEFS += -DHAVE_MOAB -DMOAB_MPAS_INSIDE_E3SM
+    MOAB_INCLUDES = -I$(MOAB_PATH)/include
+  else
+    $(error MOAB_PATH must be defined when USE_MOAB is TRUE)
+  endif
+endif
+
 RM = rm -f
 CPP = cpp -P -traditional
 FC=$(MPIFC)
@@ -48,9 +58,9 @@ FILE_OFFSET = -DOFFSET64BIT
 override CFLAGS += -DMPAS_NO_LOG_REDIRECT -DMPAS_NO_ESMF_INIT -DMPAS_ESM_SHR_CONST -DMPAS_PERF_MOD_TIMERS
 override FFLAGS += -DMPAS_NO_LOG_REDIRECT -DMPAS_NO_ESMF_INIT -DMPAS_ESM_SHR_CONST -DMPAS_PERF_MOD_TIMERS
 override CPPFLAGS += $(CPPDEFS) $(MODEL_FORMULATION) $(FILE_OFFSET) $(ZOLTAN_DEFINE) -DMPAS_NO_LOG_REDIRECT -DMPAS_NO_ESMF_INIT -DMPAS_ESM_SHR_CONST -D_MPI -DMPAS_NAMELIST_SUFFIX=$(NAMELIST_SUFFIX) -DMPAS_EXE_NAME=$(EXE_NAME) -DMPAS_PERF_MOD_TIMERS
-override CPPINCLUDES += -I$(EXEROOT)/$(COMPONENT)/source/inc -I$(INSTALL_SHAREDPATH)/include -I$(INSTALL_SHAREDPATH)/$(COMP_INTERFACE)/$(ESMFDIR)/$(NINST_VALUE)/csm_share -I$(NETCDF)/include -I$(PIO) -I$(PNETCDF)/include
+override CPPINCLUDES += -I$(EXEROOT)/$(COMPONENT)/source/inc -I$(INSTALL_SHAREDPATH)/include -I$(INSTALL_SHAREDPATH)/$(COMP_INTERFACE)/$(ESMFDIR)/$(NINST_VALUE)/csm_share -I$(NETCDF)/include -I$(PIO) -I$(PNETCDF)/include $(MOAB_INCLUDES)
 override FCINCLUDES += -I$(EXEROOT)/$(COMPONENT)/source/inc -I$(INSTALL_SHAREDPATH)/include -I$(INSTALL_SHAREDPATH)/$(COMP_INTERFACE)/$(ESMFDIR)/$(NINST_VALUE)/csm_share -I$(NETCDF)/include -I$(PIO) -I$(PNETCDF)/include
-LIBS += -L$(PIO) -L$(PNETCDF)/lib -L$(NETCDF)/lib -L$(LIBROOT) -L$(INSTALL_SHAREDPATH)/lib -lpio -lpnetcdf -lnetcdf
+LIBS += $(IMESH_LIBS) -L$(PIO) -L$(PNETCDF)/lib -L$(NETCDF)/lib -L$(LIBROOT) -L$(INSTALL_SHAREDPATH)/lib -lpio -lpnetcdf -lnetcdf 
 
 ifneq (,$(findstring FORTRANUNDERSCORE, $(CPPFLAGS)))
 ifeq (,$(findstring DUNDERSCORE, $(CPPFLAGS)))

src/driver/mpas_subdriver.F

@@ -35,6 +35,10 @@ module mpas_subdriver
    use test_core_interface
 #endif
 
+#ifdef HAVE_MOAB
+   use mpas_moabmesh
+#endif
+
    type (core_type), pointer :: corelist => null()
    type (dm_info), pointer :: dminfo
    type (domain_type), pointer :: domain_ptr
@@ -81,6 +85,10 @@ subroutine mpas_init()
       character(len=StrKIND) :: iotype
       logical :: streamsExists
       integer :: mesh_iotype
+#ifdef HAVE_MOAB
+      integer , external :: iMOAB_InitializeFortran
+#endif
+
 
       interface
          subroutine xml_stream_parser(xmlname, mgr_p, comm, ierr) bind(c)
@@ -289,6 +297,16 @@ end subroutine xml_stream_get_attributes
           call mpas_set_timeInterval(filename_interval, timeString=filename_interval_temp, ierr=ierr)
           call mpas_build_stream_filename(ref_time, start_time, filename_interval, mesh_filename_temp, blockID, mesh_filename, ierr)
       end if
+
+#ifdef HAVE_MOAB
+      ierr = iMOAB_InitializeFortran()
+      if ( ierr /= 0 ) then
+         call mpas_log_write('cannot initialize MOAB', messageType=MPAS_LOG_CRIT)
+      else
+         call mpas_log_write(' initialized MOAB', messageType=MPAS_LOG_WARN)
+      end if
+#endif
+
       call mpas_log_write(' ** Attempting to bootstrap MPAS framework using stream: ' // trim(mesh_stream))
       call mpas_bootstrap_framework_phase1(domain_ptr, mesh_filename, mesh_iotype)
 
@@ -333,6 +351,10 @@ end subroutine xml_stream_get_attributes
          call mpas_log_write('Core init failed for core '//trim(domain_ptr % core % coreName), messageType=MPAS_LOG_CRIT)
       end if
 
+#ifdef HAVE_MOAB
+      call mpas_moab_instance(domain_ptr)
+#endif
+
       call mpas_timer_stop('initialize')
 
    end subroutine mpas_init

src/framework/Makefile

@@ -36,6 +36,10 @@ OBJS = mpas_kind_types.o \
        mpas_field_accessor.o \
        mpas_log.o
 
+ifneq "$(MOAB_PATH)" ""
+  OBJS +=  mpas_moabmesh.o 
+endif
+
 all: framework $(DEPS)
 
 framework: $(OBJS)

src/framework/mpas_moabmesh.F

@@ -0,0 +1,221 @@
+
+
+module mpas_moabmesh
+!  use, intrinsic :: ISO_C_BINDING
+
+  use mpas_log     
+  use mpas_derived_types, only: dm_info, domain_type
+  use mpas_field_routines
+  use mpas_sort
+  use mpas_stream_manager
+  use mpas_pool_routines
+  !use mpas_vector_operations
+#include "moab/MOABConfig.h"
+
+#ifdef MOAB_MPAS_INSIDE_E3SM
+  use seq_comm_mct, only: MPOID, OCNID
+#else
+  integer, public :: MPOID !   app id on moab side, for MPAS ocean
+#endif
+  implicit none
+
+  contains
+
+  SUBROUTINE errorout(ierr, message)
+  integer ierr
+  character*(*) message
+  if (ierr.ne.0) then
+    print *, message
+    call exit (1)
+  end if
+  return
+  end subroutine
+
+  subroutine mpas_moab_instance(domain)
+
+    type (domain_type), intent(inout) :: domain
+
+    type (block_type), pointer :: block
+    type (mpas_pool_type), pointer :: meshPool
+    integer, pointer :: nCells, nVertices, maxEdges
+    integer  :: pid, nblocks
+    integer, dimension(:,:), pointer :: verticesOnCell
+    integer, dimension(:), pointer :: nEdgesOnCell
+    integer, dimension(:), pointer :: indexToVertexID, indexToCellID
+    real(kind=RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex
+    real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
+    logical, pointer :: on_a_sphere, is_periodic
+    real(kind=RKIND), pointer :: x_period, y_period
+    integer, pointer :: nCellsSolve, nEdgesSolve, nVerticesSolve
+
+    integer , external :: iMOAB_RegisterFortranApplication, &
+        iMOAB_CreateVertices, iMOAB_WriteMesh, iMOAB_CreateElements, &
+        iMOAB_ResolveSharedEntities, iMOAB_DetermineGhostEntities, &
+        iMOAB_DefineTagStorage,  iMOAB_SetIntTagStorage , &
+        iMOAB_UpdateMeshInfo
+    integer :: c_comm, i1, j1, ic, lastvertex
+    character*12 appname
+    integer :: ierr, num_verts_in_cells, ext_comp_id
+    real(kind=RKIND), allocatable, target :: moab_vert_coords(:)
+    integer, allocatable, target :: indexUsed(:), invMap(:), localIds(:)
+    integer dimcoord, dimen, mbtype, block_ID, proc_id
+    integer  ,allocatable , target :: all_connects(:)
+    character*100 outfile, wopts, localmeshfile, lnum, tagname
+    integer tagtype, numco, tag_sto_len, ent_type, tagindex, currentVertex
+
+    c_comm = domain % dminfo % comm
+    appname = 'MPAS_MESH'//CHAR(0)
+    MPOID = -1 ! initialized negative , so we know if this register passes
+#ifdef MOAB_MPAS_INSIDE_E3SM
+    ext_comp_id = OCNID(1) ! should be 17
+#else
+    ext_comp_id = 17  ! some number less than 33 ! maybe this should be input
+#endif
+    ierr = iMOAB_RegisterFortranApplication(appname, c_comm, ext_comp_id, pid)
+    MPOID = pid ! this is exported, need for send to work
+    call errorout(ierr, 'fail to register MPAS_MOAB mesh')
+    proc_id = domain % dminfo % my_proc_id
+    call mpas_log_write('MOAB MPAS app pid:  $i  task  $i ', intArgs=(/pid, proc_id/) )
+
+!   blocks should be merged if there is more than one block per task
+    nblocks = 0
+    block => domain % blocklist
+    do while (associated(block)) !{{{
+       nblocks = nblocks + 1
+        ! allocate scratch memory 
+       call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
+       call mpas_pool_get_array(meshPool, 'xVertex', xVertex)
+       call mpas_pool_get_array(meshPool, 'yVertex', yVertex)
+       call mpas_pool_get_array(meshPool, 'zVertex', zVertex)
+       call mpas_pool_get_dimension(meshPool, 'nVertices', nVertices)
+       call mpas_pool_get_dimension(meshPool, 'maxEdges', maxEdges)
+       call mpas_pool_get_config(meshPool, 'on_a_sphere', on_a_sphere)
+       call mpas_pool_get_config(meshPool, 'is_periodic', is_periodic)
+       call mpas_pool_get_config(meshPool, 'x_period', x_period)
+       call mpas_pool_get_config(meshPool, 'y_period', y_period)
+       call mpas_pool_get_array(meshPool, 'verticesOnCell', verticesOnCell)
+       call mpas_pool_get_array(meshPool, 'indexToVertexID', indexToVertexID)
+       call mpas_pool_get_array(meshPool, 'indexToCellID', indexToCellID)
+       call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
+       call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
+       call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve)
+       call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve)
+       call mpas_pool_get_dimension(meshPool, 'nVerticesSolve', nVerticesSolve)
+!       call mpas_pool_get_array(meshPool, 'xCell', xCell)
+!       call mpas_pool_get_array(meshPool, 'yCell', yCell)
+!       call mpas_pool_get_array(meshPool, 'zCell', zCell)
+
+       call mpas_log_write(' MOAB instance: number of vertices:: $i number of cells:: $i solve: v:$i c:$i', intArgs=(/nVertices, nCells, nVerticesSolve, nCellsSolve/) )
+        !!
+       allocate(indexUsed(nVertices), invMap(nVertices) ) ! conservative, invMap should be smaller
+       indexUsed = 0
+       invMap = 0
+    ! fill now connectivity array, nCellsSolve; fist pad to max nc
+       num_verts_in_cells = nCellsSolve * maxEdges
+       allocate(all_connects(num_verts_in_cells))
+!    collect all vertices, and also pad 
+       j1 = 0
+       do ic=1, nCellsSolve
+         do i1 = 1, nEdgesOnCell(ic)
+           j1 = j1 + 1
+           all_connects(j1) = verticesOnCell( i1, ic)
+           indexUsed(all_connects(j1)) = 1
+         enddo
+         lastvertex = verticesOnCell( nEdgesOnCell (ic), ic)
+         ! pad the rest with the last vertex
+         do i1 =   nEdgesOnCell (ic) + 1, maxEdges
+           j1 = j1 + 1
+           all_connects(j1) = lastvertex ! repeat the last vertex (pad)
+         enddo
+         ! call mpas_log_write('cell: $i  v:: $i $i $i $i $i $i', intArgs=(/ic, all_connects(j1-5), all_connects(j1-4), all_connects(j1-3), all_connects(j1-2), all_connects(j1-1), all_connects(j1)/) )
+       enddo
+
+       currentVertex = 0
+       do i1 = 1, nVertices
+         if (indexUsed(i1) > 0) then
+           currentVertex = currentVertex + 1
+           indexUsed(i1) = currentVertex
+           invMap(currentVertex) = i1
+         endif
+       enddo
+       !! convert all_connects to indexUsed
+       do i1 = 1, num_verts_in_cells
+         all_connects(i1) = indexUsed( all_connects(i1) )
+       enddo
+       allocate(moab_vert_coords(3*currentVertex))
+       do i1 =1, currentVertex
+         moab_vert_coords(3*i1-2) = xVertex(invMap(i1))
+         moab_vert_coords(3*i1-1) = yVertex(invMap(i1))
+         moab_vert_coords(3*i1  ) = zVertex(invMap(i1))
+         ! call mpas_log_write('i:: $i coords:: $r $r $r $r', intArgs=(/i1/), realArgs=(/moab_vert_coords(3*i1-2),moab_vert_coords(3*i1-1), moab_vert_coords(3*i1)/) )
+       enddo
+       dimcoord = 3*currentVertex
+       dimen = 3
+       ierr = iMOAB_CreateVertices(pid, dimcoord, dimen, moab_vert_coords)
+       call errorout(ierr, 'fail to create vertices')
+       call mpas_log_write(' MOAB instance: created $i vertices on local proc $i ',intArgs=(/currentVertex, proc_id/))
+!  so we know we used only currentvertex vertices in the pool (the rest are in halo)
+       mbtype = 4 !  polygon
+
+       block_ID = 100*proc_id + nblocks ! we should have only one block right now
+
+       ierr = iMOAB_CreateElements( pid, nCellsSolve, mbtype, maxEdges, all_connects, block_ID );
+       call errorout(ierr, 'fail to create polygons')
+!      set the global id for vertices
+!   first, retrieve the tag
+       tagname='GLOBAL_ID'//CHAR(0)
+       tagtype = 0  ! dense, integer
+       numco = 1
+       ierr = iMOAB_DefineTagStorage(pid, tagname, tagtype, numco,  tagindex )
+       call errorout(ierr, 'fail to get global id tag')
+! now set the values
+       ent_type = 0 ! vertex type
+       allocate(localIds(currentVertex))
+       do i1 = 1, currentVertex
+         localIds(i1) = indexToVertexID (invMap(i1))
+       enddo
+       ierr = iMOAB_SetIntTagStorage ( pid, tagname, currentVertex , ent_type, localIds )
+       call errorout(ierr, 'fail to set global id tag for vertices')
+       ! set global id tag for elements
+       ent_type = 1 ! now set the global id tag on elements
+       ierr = iMOAB_SetIntTagStorage ( pid, tagname, nCellsSolve, ent_type, indexToCellID)
+       call errorout(ierr, 'fail to set global id tag for polygons')
+       ! get next block
+#ifdef MPAS_DEBUG
+       if (proc_id.lt. 5) then
+         write(lnum,"(I0.2)")proc_id
+         localmeshfile = 'ownedOcn_'//trim(lnum)// '.h5m' // CHAR(0)
+         wopts = CHAR(0)
+         ierr = iMOAB_WriteMesh(pid, localmeshfile, wopts)
+         call errorout(ierr, 'fail to write local mesh file')
+       endif
+#endif
+       ierr = iMOAB_ResolveSharedEntities( pid, currentVertex, localIds );
+       call errorout(ierr, 'fail to resolve shared entities')
+
+       deallocate (moab_vert_coords)
+       deallocate (all_connects)
+       deallocate (indexUsed)
+       deallocate (invMap)
+       deallocate (localIds)
+       block => block % next
+
+    end do !}}}
+
+
+    if (nblocks .ne. 1) then
+      call errorout(1, 'more than one block per task')
+    endif 
+    ierr = iMOAB_UpdateMeshInfo(pid)
+    call errorout(ierr, 'fail to update mesh info')
+    !     write out the mesh file to disk, in parallel
+    outfile = 'wholeMPAS.h5m'//CHAR(0)
+    wopts   = 'PARALLEL=WRITE_PART'//CHAR(0)
+    ierr = iMOAB_WriteMesh(pid, outfile, wopts)
+    call errorout(ierr, 'fail to write the mesh file')
+
+
+  end subroutine mpas_moab_instance
+end module mpas_moabmesh
+
+