diff --git a/Makefile b/Makefile index baceb05cb7..3693ce0726 100644 --- a/Makefile +++ b/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)" diff --git a/src/Makefile.in.ACME b/src/Makefile.in.ACME index 866ee1fedf..b8ca28e513 100644 --- a/src/Makefile.in.ACME +++ b/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))) diff --git a/src/driver/mpas_subdriver.F b/src/driver/mpas_subdriver.F index fe80338ceb..4631c8536b 100644 --- a/src/driver/mpas_subdriver.F +++ b/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 diff --git a/src/framework/Makefile b/src/framework/Makefile index 4e635bbbc9..4a05b1c320 100644 --- a/src/framework/Makefile +++ b/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) diff --git a/src/framework/mpas_moabmesh.F b/src/framework/mpas_moabmesh.F new file mode 100644 index 0000000000..dc7ef18f83 --- /dev/null +++ b/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 + +