imoab_coupler_fortran.F90

Go to the documentation of this file.

! This program shows how to do an imoab_coupler type test in fortran 90
! the atm and ocean files are read from repo, meshes are migrated to coupler pes,
! the intx is caried on, weight generation,; tag migration and projection,
! migrate back to ocean pes, and compare against baseline1.txt
! cannot change source/target files, or the layout, just test the imoab
! calls, that need to give the same results as C=++ equivalent test (imoab_coupler)
! between atm SE and ocn FV
 
SUBROUTINE errorout(ierr, message)
 integer ierr
 character*(*) message
 if (ierr .ne. 0) then
 print *, message
 call exit(1)
 end if
 return
end
!
#include "moab/MOABConfig.h"
 
#ifndef MOAB_MESH_DIR
#error Specify MOAB_MESH_DIR path
#endif
 
program imoab_coupler_fortran
 
 use iso_c_binding
 use imoab
 
#include "mpif.h"
#include "moab/MOABConfig.h"
 !implicit none
 integer :: m ! for number of arguments ; if less than 1, exit
 integer :: global_comm
 integer :: rankinglobalcomm, numprocesses
 integer :: ierr
 integer :: my_id, num_procs
 integer :: jgroup ! group for global comm
 character(:), allocatable :: atmfilename
 character(:), allocatable :: ocnfilename
 character(:), allocatable :: baselinefilename
 character(:), allocatable :: readopts, filewriteoptions
 character :: appname*128
 character(:), allocatable :: weights_identifier1
 character(:), allocatable :: disc_methods1, disc_methods2, dof_tag_names1, dof_tag_names2
 integer :: disc_orders1, disc_orders2
 character(:), allocatable :: atmocn_map_file_name, intx_from_file_identifier
 
 ! all groups and comms are on the same number of processes, for simplicity
 integer :: atmgroup, atmcomm, ocngroup, ocncomm, cplgroup, cplcomm
 integer :: atmcoucomm, ocncoucomm
 
 integer :: cmpatm, cplatm, cmpocn, cplocn, atmocnid ! global IDs
 integer :: cmpatmpid, cplatmpid ! iMOAB app ids
 integer :: cmpocnpid, cplocnpid ! iMOAB app ids
 integer :: cplatmocnpid ! intx pid
 integer :: nghlay, partscheme, context_id
 integer :: fnobubble, fmonotonetypeid, fvolumetric, fnoconserve, fvalidate, finversedistancemap
 
 integer, dimension(2) :: tagindex
 integer, dimension (2) :: tagtypes! { DENSE_DOUBLE, DENSE_DOUBLE }
 integer :: atmcompndofs ! = disc_orders[0] * disc_orders[0],
 integer :: ocncompndofs ! = 1 /*FV*/
 character(:), allocatable :: bottomfields, bottomprojectedfields
 integer, dimension(3) :: nverts, nelem, nblocks, nsbc, ndbc
 double precision, allocatable :: vals(:) ! to set the double values to 0
 integer :: i ! for loops
 integer :: storleng, eetype ! for tags defs
 character(:), allocatable :: concat_fieldname, concat_fieldnamet, outputfileocn
 integer :: tagindexin2 ! not really needed
 integer :: dummycpl, dummyrc, dummytype
 real*8 :: boxeps
 integer :: gnomonic
 
 cmpatm = 5
 cplatm = 6
 cmpocn = 17
 cplocn = 18
 atmocnid = 618
 
 call mpi_init(ierr)
 call errorout(ierr,'mpi_init')
 call mpi_comm_rank (mpi_comm_world, my_id, ierr)
 call errorout(ierr, 'fail to get MPI rank')
 
 call mpi_comm_size (mpi_comm_world, num_procs, ierr)
 call errorout(ierr, 'fail to get MPI size')
 call mpi_comm_dup(mpi_comm_world, global_comm, ierr)
 call errorout(ierr, 'fail to get global comm duplicate')
 
 call mpi_comm_group( global_comm, jgroup, ierr ); ! all processes in jgroup
 call errorout(ierr, 'fail to get joint group')
 ! readopts( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS" )
 ! readoptsLnd( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" )
 atmfilename = &
 moab_mesh_dir &
 //'unittest/wholeATM_T.h5m'//c_null_char
 ocnfilename = &
 moab_mesh_dir &
 //'unittest/recMeshOcn.h5m'//c_null_char
 baselinefilename = &
 moab_mesh_dir &
 //'unittest/baseline1.txt'//c_null_char
 
 ! all comms span the whole world, for simplicity
 atmcomm = mpi_comm_null
 ocncomm = mpi_comm_null
 cplcomm = mpi_comm_null
 atmcoucomm = mpi_comm_null
 ocncoucomm = mpi_comm_null
 call mpi_comm_dup(global_comm, atmcomm, ierr)
 call mpi_comm_dup(global_comm, ocncomm, ierr)
 call mpi_comm_dup(global_comm, cplcomm, ierr)
 call mpi_comm_dup(global_comm, atmcoucomm, ierr)
 call mpi_comm_dup(global_comm, ocncoucomm, ierr)
 
 ! all groups of interest are easy breezy
 call mpi_comm_group( atmcomm, atmgroup, ierr )
 call mpi_comm_group( ocncomm, ocngroup, ierr )
 call mpi_comm_group( cplcomm, cplgroup, ierr )
 
 readopts ='PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS'//c_null_char
 nghlay = 0 ! no ghost layers
#ifdef MOAB_HAVE_ZOLTAN
 partscheme = 2 ! RCB with zoltan
#else
 partscheme = 0 ! Trivial partitioner
#endif
 
 if (my_id .eq. 0) then
 print *, ' number of tasks: ', num_procs
 print *, ' Atm file: ', atmfilename
 print *, ' Ocn file: ', ocnfilename
 print *, ' baseline file: ', baselinefilename
 print *, ' using partitioner: ', partscheme
 end if
 
 ierr = imoab_initialize()
 appname = 'ATM'//c_null_char
 ierr = imoab_registerapplication(appname, atmcomm, cmpatm, cmpatmpid)
 appname = 'ATMX'//c_null_char
 ierr = imoab_registerapplication(appname, cplcomm, cplatm, cplatmpid)
 appname = 'OCN'//c_null_char
 ierr = imoab_registerapplication(appname, ocncomm, cmpocn, cmpocnpid)
 appname = 'OCNX'//c_null_char
 ierr = imoab_registerapplication(appname, cplcomm, cplocn, cplocnpid)
 
 appname = 'ATMOCN'//c_null_char
 ierr = imoab_registerapplication(appname, cplcomm, atmocnid, cplatmocnpid)
 
 ! read atm and migrate
 if (atmcomm .NE. mpi_comm_null) then
 ierr = imoab_loadmesh(cmpatmpid, atmfilename, readopts, nghlay)
 call errorout(ierr, 'fail to load atm')
 ierr = imoab_sendmesh(cmpatmpid, atmcoucomm, cplgroup, cplatm, partscheme)
 call errorout(ierr, 'fail to send atm')
 end if
 if (cplcomm .NE. mpi_comm_null) then
 ierr = imoab_receivemesh(cplatmpid, atmcoucomm, atmgroup, cmpatm) !
 call errorout(ierr, 'fail to receive atm')
 end if
 
 ! we can now free the sender buffers
 if (atmcomm .NE. mpi_comm_null) then
 context_id = cplatm
 ierr = imoab_freesenderbuffers(cmpatmpid, context_id)
 call errorout(ierr, 'fail to free atm buffers')
 end if
 
 ! read ocn and migrate
 if (ocncomm .NE. mpi_comm_null) then
 ierr = imoab_loadmesh(cmpocnpid, ocnfilename, readopts, nghlay)
 call errorout(ierr, 'fail to load ocn')
 ierr = imoab_sendmesh(cmpocnpid, ocncoucomm, cplgroup, cplocn, partscheme)
 call errorout(ierr, 'fail to send ocn')
 end if
 if (cplcomm .NE. mpi_comm_null) then
 ierr = imoab_receivemesh(cplocnpid, ocncoucomm, ocngroup, cmpocn) !
 call errorout(ierr, 'fail to receive ocn')
 end if
 
 ! we can now free the sender buffers
 if (ocncomm .NE. mpi_comm_null) then
 context_id = cplocn
 ierr = imoab_freesenderbuffers(cmpocnpid, context_id)
 call errorout(ierr, 'fail to free ocn buffers')
 end if
 
 if (cplcomm .NE. mpi_comm_null) then
 ierr = imoab_computemeshintersectiononsphere(cplatmpid, cplocnpid, cplatmocnpid)
 ! coverage mesh was computed here, for cplAtmPID, atm on coupler pes
 ! basically, atm was redistributed according to target (ocean) partition, to "cover" the
 !ocean partitions check if intx valid, write some h5m intx file
 call errorout(ierr, 'cannot compute intersection')
 end if
 
 if (atmcoucomm .NE. mpi_comm_null) then
 ! the new graph will be for sending data from atm comp to coverage mesh.
 ! it involves initial atm app; cmpAtmPID; also migrate atm mesh on coupler pes, cplAtmPID
 ! results are in cplAtmOcnPID, intx mesh; remapper also has some info about coverage mesh
 ! after this, the sending of tags from atm pes to coupler pes will use the new par comm
 ! graph, that has more precise info about what to send for ocean cover ; every time, we
 ! will use the element global id, which should uniquely identify the element
 ierr = imoab_coveragegraph(atmcoucomm, cmpatmpid, cplatmpid, cplatmocnpid, cmpatm, cplatm, cplocn) ! it happens over joint communicator
 call errorout(ierr, 'cannot recompute direct coverage graph for ocean')
 end if
 
 weights_identifier1 = 'scalar'//c_null_char
 disc_methods1 = 'cgll'//c_null_char
 disc_methods2 = 'fv'//c_null_char
 disc_orders1 = 4
 disc_orders2 = 1
 dof_tag_names1 = 'GLOBAL_DOFS'//c_null_char
 dof_tag_names2 = 'GLOBAL_ID'//c_null_char
 ! fMonotoneTypeID = 0, fVolumetric = 0, fValidate = 1, fNoConserve = 0, fNoBubble = 1
 fnobubble = 1
 fmonotonetypeid = 0
 fvolumetric = 0
 fnoconserve = 0
 fvalidate = 0
 finversedistancemap = 0
 
 if (cplcomm .NE. mpi_comm_null) then
 
 ierr = imoab_computescalarprojectionweights( &
 cplatmocnpid, weights_identifier1, disc_methods1, disc_orders1, &
 disc_methods2, disc_orders2, ""//c_null_char, fnobubble, fmonotonetypeid, fvolumetric, &
 finversedistancemap, fnoconserve, &
 fvalidate, dof_tag_names1, dof_tag_names2)
 call errorout(ierr, 'cannot compute scalar projection weights')
 
 ierr = imoab_computescalarprojectionweights( &
 cplatmocnpid, "bilinear"//c_null_char, "fv"//c_null_char, 1, &
 "fv"//c_null_char, 1, "bilin"//c_null_char, fnobubble, fmonotonetypeid, fvolumetric, &
 finversedistancemap, fnoconserve, &
 fvalidate, "GLOBAL_ID"//c_null_char, "GLOBAL_ID"//c_null_char)
 call errorout(ierr, 'cannot compute scalar projection weights')
 
#ifdef MOAB_HAVE_NETCDF
 atmocn_map_file_name = 'atm_ocn_map_f.nc'//c_null_char
 ierr = imoab_writemappingweightstofile( cplatmocnpid, weights_identifier1, atmocn_map_file_name)
 call errorout(ierr, 'failed to write map file to disk')
 intx_from_file_identifier = 'map-from-file'//c_null_char
 dummycpl = -1
 dummyrc = -1
 dummytype = 0
 ierr = imoab_loadmappingweightsfromfile( cplatmocnpid, dummycpl, dummyrc, dummytype, &
 intx_from_file_identifier, atmocn_map_file_name)
 call errorout(ierr, 'failed to load map file from disk')
#endif
 end if
 
 ! start copy
 tagtypes(1) = 1 ! somehow, DENSE_DOUBLE give 0, while it should be 1; maybe moab::DENSE_DOUBLE ?
 tagtypes(2) = 1 ! ! DENSE_DOUBLE
 atmcompndofs = disc_orders1*disc_orders1
 ocncompndofs = 1 ! /*FV*/
 
 bottomfields = 'a2oTbot:a2oUbot:a2oVbot'//c_null_char
 bottomprojectedfields = 'a2oTbot_proj:a2oUbot_proj:a2oVbot_proj'//c_null_char
 
 if (cplcomm .NE. mpi_comm_null) then
 ierr = imoab_definetagstorage(cplatmpid, bottomfields, tagtypes(1), atmcompndofs, tagindex(1))
 call errorout(ierr, 'failed to define the field tags a2oTbot:a2oUbot:a2oVbot ')
 ierr = imoab_definetagstorage(cplocnpid, bottomprojectedfields, tagtypes(2), ocncompndofs, tagindex(2))
 call errorout(ierr, 'failed to define the field tags a2oTbot_proj:a2oUbot_proj:a2oVbot_proj')
 end if
 
 ! make the tag 0, to check we are actually sending needed data
 if (cplatmpid .ge. 0) then
 
 ! Each process in the communicator will have access to a local mesh instance, which
 ! will contain the original cells in the local partition and ghost entities. Number of
 ! vertices, primary cells, visible blocks, number of sidesets and nodesets boundary
 ! conditions will be returned in numProcesses 3 arrays, for local, ghost and total
 ! numbers.
 
 ierr = imoab_getmeshinfo(cplatmpid, nverts, nelem, nblocks, nsbc, ndbc)
 call errorout(ierr, 'failed to get num primary elems')
 storleng = nelem(3)*atmcompndofs*3 ! 3 tags
 allocate (vals(storleng))
 eetype = 1 ! double type
 
 do i = 1, storleng
 vals(:) = 0.
 end do
 
 ! set the tag values to 0.0
 ierr = imoab_setdoubletagstorage(cplatmpid, bottomfields, storleng, eetype, vals)
 call errorout(ierr, 'cannot make tag nul')
 
 end if
 
 ! Define the field variables to project
 concat_fieldname = 'a2oTbot:a2oUbot:a2oVbot'//c_null_char
 concat_fieldnamet = 'a2oTbot_proj:a2oUbot_proj:a2oVbot_proj'//c_null_char
 
 if (atmcomm .NE. mpi_comm_null) then
 
 ! As always, use nonblocking sends
 ! this is for projection to ocean:
 ierr = imoab_sendelementtag(cmpatmpid, concat_fieldname, atmcoucomm, cplocn)
 call errorout(ierr, 'cannot send tag values')
 
 end if
 
 if (cplcomm .NE. mpi_comm_null) then
 !// receive on atm on coupler pes, that was redistributed according to coverage
 ierr = imoab_receiveelementtag(cplatmpid, concat_fieldname, atmcoucomm, cplocn)
 call errorout(ierr, 'cannot receive tag values')
 end if
 
 ! we can now free the sender buffers
 if (atmcomm .NE. mpi_comm_null) then
 
 ierr = imoab_freesenderbuffers(cmpatmpid, cplocn) !context is for ocean
 call errorout(ierr, 'cannot free buffers used to resend atm tag towards the coverage mesh')
 
 end if
 if (cplcomm .ne. mpi_comm_null) then
 
 outputfileocn = "AtmOnCplF.h5m"//c_null_char
 filewriteoptions = 'PARALLEL=WRITE_PART'//c_null_char
 ierr = imoab_writemesh(cplatmpid, outputfileocn, filewriteoptions)
 call errorout(ierr, 'could not write AtmOnCpl.h5m to disk')
 
 end if
 if (cplcomm .ne. mpi_comm_null) then
 
 ! We have the remapping weights now. Let us apply the weights onto the tag we defined
 ! on the source mesh and get the projection on the target mesh
 ierr = imoab_applyscalarprojectionweights(cplatmocnpid, weights_identifier1, concat_fieldname, &
 concat_fieldnamet)
 call errorout(ierr, 'failed to compute projection weight application')
 
 outputfileocn = "OcnOnCplF.h5m"//c_null_char
 filewriteoptions = 'PARALLEL=WRITE_PART'//c_null_char
 ierr = imoab_writemesh(cplocnpid, outputfileocn, filewriteoptions)
 call errorout(ierr, 'could not write OcnOnCpl.h5m to disk')
 
 end if
 ! send the projected tag back to ocean pes, with send/receive tag
 ! first makje sure the tags are defined, otherwise they cannot be received
 if (ocncomm .ne. mpi_comm_null) then
 
 ierr = imoab_definetagstorage(cmpocnpid, bottomprojectedfields, tagtypes(2), ocncompndofs, tagindexin2)
 call errorout(ierr, 'failed to define the field tag for receiving back the tag a2oTbot_proj, on ocn pes')
 
 end if
 
 ! send the tag to ocean pes, from ocean mesh on coupler pes
 ! from couComm, using common joint comm ocn_coupler
 ! as always, use nonblocking sends
 ! original graph (context is -1_
 if (cplcomm .ne. mpi_comm_null) then
 context_id = cmpocn
 ierr = imoab_sendelementtag(cplocnpid, concat_fieldnamet, ocncoucomm, context_id)
 call errorout(ierr, 'cannot send tag values back to ocean pes')
 end if
 
 if (ocncomm .ne. mpi_comm_null) then
 context_id = cplocn
 ierr = imoab_receiveelementtag(cmpocnpid, concat_fieldnamet, ocncoucomm, context_id)
 call errorout(ierr, 'cannot receive tag values from ocean mesh on coupler pes')
 end if
 
 if (cplcomm .ne. mpi_comm_null) then
 context_id = cmpocn
 ierr = imoab_freesenderbuffers(cplocnpid, context_id)
 call errorout(ierr, 'cannot free sender buffers on coupler')
 end if
 
 if (ocncomm .ne. mpi_comm_null) then
 
 outputfileocn = "OcnWithProjF.h5m"//c_null_char
 filewriteoptions = 'PARALLEL=WRITE_PART'//c_null_char
 if (my_id .eq. 0) then
 print *, ' Writing ocean mesh file with projected solution to disk: ', outputfileocn
 end if
 ierr = imoab_writemesh(cmpocnpid, outputfileocn, filewriteoptions)
 call errorout(ierr, 'could not write OcnWithProjF.h5m to disk')
 
 end if
 
 ! free up resources
 ierr = imoab_deregisterapplication(cplatmocnpid)
 call errorout(ierr, 'could not de-register OCN component')
 
 ierr = imoab_deregisterapplication(cplocnpid)
 call errorout(ierr, 'could not de-register OCN component')
 ierr = imoab_deregisterapplication(cmpocnpid)
 call errorout(ierr, 'could not de-register OCN component')
 ierr = imoab_deregisterapplication(cplatmpid)
 call errorout(ierr, 'could not de-register OCN component')
 ierr = imoab_deregisterapplication(cmpatmpid)
 call errorout(ierr, 'could not de-register OCN component')
 
 ! Free all MPI datastructures
 call mpi_comm_free(atmcomm, ierr)
 call mpi_group_free(atmgroup, ierr)
 call mpi_comm_free(ocncomm, ierr)
 call mpi_group_free(ocngroup, ierr)
 call mpi_comm_free(cplcomm, ierr)
 call mpi_group_free(cplgroup, ierr)
 
 call mpi_comm_free(atmcoucomm, ierr)
 call mpi_comm_free(ocncoucomm, ierr)
 call mpi_comm_free(global_comm, ierr)
 call mpi_group_free(jgroup, ierr)
 call mpi_finalize(ierr)
 
end program imoab_coupler_fortran