|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
#include "moab/Core.hpp"#include "moab_mpi.h"#include "moab/ParallelComm.hpp"#include "MBParallelConventions.h"#include "moab/iMOAB.h"#include "TestUtil.hpp"#include "moab/CpuTimer.hpp"#include "moab/ProgOptions.hpp"#include <iostream>#include <sstream>#include "imoab_coupler_utils.hpp"
Include dependency graph for imoab_phatm_ocn_coupler.cpp:Go to the source code of this file.
Defines | |
| #define | ENABLE_ATMOCN_COUPLING |
| #define | ENABLE_ATMLND_COUPLING |
Functions | |
| int | main (int argc, char *argv[]) |
| #define ENABLE_ATMLND_COUPLING |
Definition at line 46 of file imoab_phatm_ocn_coupler.cpp.
| #define ENABLE_ATMOCN_COUPLING |
Definition at line 42 of file imoab_phatm_ocn_coupler.cpp.
| int main | ( | int | argc, |
| char * | argv[] | ||
| ) |
Definition at line 48 of file imoab_phatm_ocn_coupler.cpp.
References ProgOptions::addOpt(), atmFilename, CHECKIERR, cmpatm, context, create_group_and_comm(), create_joint_comm_group(), DENSE_DOUBLE, endG1, endG2, fileWriteOptions(), ierr, iMOAB_AppID, iMOAB_DefineTagStorage(), iMOAB_DeregisterApplication(), iMOAB_Finalize(), iMOAB_GetMeshInfo(), iMOAB_Initialize(), iMOAB_LoadMesh(), iMOAB_RegisterApplication(), iMOAB_SetDoubleTagStorage(), iMOAB_WriteMesh(), jgroup, MPI_COMM_WORLD, nghlay, numProcesses, ProgOptions::parseCommandLine(), POP_TIMER, PUSH_TIMER, rankInGlobalComm, readopts(), setup_component_coupler_meshes(), startG1, startG2, typeA, and typeB.
{
int ierr;
int rankInGlobalComm, numProcesses;
MPI_Group jgroup;
std::string readopts( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS" );
std::string readoptsPhysAtm( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" );
// Timer data
moab::CpuTimer timer;
double timer_ops;
std::string opName;
int repartitioner_scheme = 0;
#ifdef MOAB_HAVE_ZOLTAN
repartitioner_scheme = 2; // use the graph partitioner in that case
#endif
MPI_Init( &argc, &argv );
MPI_Comm_rank( MPI_COMM_WORLD, &rankInGlobalComm );
MPI_Comm_size( MPI_COMM_WORLD, &numProcesses );
MPI_Comm_group( MPI_COMM_WORLD, &jgroup ); // all processes in jgroup
std::string atmFilename = TestDir + "unittest/wholeATM_T.h5m"; // we should use only mesh from here
std::string atmPhysMesh =
TestDir + "unittest/AtmPhys_01.h5m"; // it has some data associated to vertices, T_ph, u_ph, v_ph
// we will eventually project that data to ocean mesh, after intx atm/ocn
// on a regular case, 5 ATM, 6 CPLATM (ATMX), 17 OCN , 18 CPLOCN (OCNX) ;
// intx atm/ocn is not in e3sm yet, give a number
// 6 * 100+ 18 = 618 : atmocnid
// 9 LND, 10 CPLLND
// 6 * 100 + 10 = 610 atmlndid:
// cmpatm is for atm on atm pes
// cmpocn is for ocean, on ocean pe
// cplatm is for atm on coupler pes
// cplocn is for ocean on coupler pes
// atmocnid is for intx atm / ocn on coupler pes
//
int cmpatm = 5,
cplatm = 6; // component ids are unique over all pes, and established in advance;
int cmpPhysAtm = 105; // different from atm spectral ?
#ifdef ENABLE_ATMOCN_COUPLING
std::string ocnFilename = TestDir + "unittest/recMeshOcn.h5m";
int rankInOcnComm = -1;
int cmpocn = 17, cplocn = 18,
atmocnid = 618; // component ids are unique over all pes, and established in advance;
#endif
#ifdef ENABLE_ATMLND_COUPLING
std::string lndFilename = TestDir + "unittest/wholeLnd.h5m";
int cpllnd = 10,
cmplnd = 9; // component ids are unique over all pes, and established in advance;
#endif
int rankInCouComm = -1;
int nghlay = 0; // number of ghost layers for loading the file
int startG1 = 0, startG2 = 0, endG1 = numProcesses - 1, endG2 = numProcesses - 1, startG3 = startG1, endG3 = endG1;
int startG4 = startG1, endG4 = endG1; // these are for coupler layout
int context_id = -1; // used now for freeing buffers
// default: load atm on 2 proc, ocean on 2, land on 2; migrate to 2 procs, then compute intx
// later, we need to compute weight matrix with tempestremap
ProgOptions opts;
opts.addOpt< std::string >( "atmosphere,t", "atm mesh filename (source)", &atmFilename );
#ifdef ENABLE_ATMOCN_COUPLING
opts.addOpt< std::string >( "ocean,m", "ocean mesh filename (target)", &ocnFilename );
#endif
#ifdef ENABLE_ATMLND_COUPLING
opts.addOpt< std::string >( "land,l", "land mesh filename (target)", &lndFilename );
#endif
opts.addOpt< std::string >( "physgrid,q", "physics grid file", &atmPhysMesh );
opts.addOpt< int >( "startAtm,a", "start task for atmosphere layout", &startG1 );
opts.addOpt< int >( "endAtm,b", "end task for atmosphere layout", &endG1 );
#ifdef ENABLE_ATMOCN_COUPLING
opts.addOpt< int >( "startOcn,c", "start task for ocean layout", &startG2 );
opts.addOpt< int >( "endOcn,d", "end task for ocean layout", &endG2 );
#endif
#ifdef ENABLE_ATMLND_COUPLING
opts.addOpt< int >( "startLnd,e", "start task for land layout", &startG3 );
opts.addOpt< int >( "endLnd,f", "end task for land layout", &endG3 );
#endif
opts.addOpt< int >( "startCoupler,g", "start task for coupler layout", &startG4 );
opts.addOpt< int >( "endCoupler,j", "end task for coupler layout", &endG4 );
opts.addOpt< int >( "partitioning,p", "partitioning option for migration", &repartitioner_scheme );
opts.parseCommandLine( argc, argv );
char fileWriteOptions[] = "PARALLEL=WRITE_PART";
if( !rankInGlobalComm )
{
std::cout << " atm file: " << atmFilename << "\n on tasks : " << startG1 << ":" << endG1 <<
#ifdef ENABLE_ATMOCN_COUPLING
"\n ocn file: " << ocnFilename << "\n on tasks : " << startG2 << ":" << endG2 <<
#endif
#ifdef ENABLE_ATMLND_COUPLING
"\n land file: " << lndFilename << "\n on tasks : " << startG3 << ":" << endG3 <<
#endif
"\n atm phys file: " << atmPhysMesh << "\n on tasks : " << startG1 << ":" << endG1 <<
"\n partitioning (0 trivial, 1 graph, 2 geometry) " << repartitioner_scheme << "\n ";
}
// load files on 3 different communicators, groups
// first groups has task 0, second group tasks 0 and 1
// coupler will be on joint tasks, will be on a third group (0 and 1, again)
MPI_Group atmPEGroup;
MPI_Comm atmComm;
ierr = create_group_and_comm( startG1, endG1, jgroup, &atmPEGroup, &atmComm );
CHECKIERR( ierr, "Cannot create atm MPI group and communicator " )
#ifdef ENABLE_ATMOCN_COUPLING
MPI_Group ocnPEGroup;
MPI_Comm ocnComm;
ierr = create_group_and_comm( startG2, endG2, jgroup, &ocnPEGroup, &ocnComm );
CHECKIERR( ierr, "Cannot create ocn MPI group and communicator " )
#endif
#ifdef ENABLE_ATMLND_COUPLING
MPI_Group lndPEGroup;
MPI_Comm lndComm;
ierr = create_group_and_comm( startG3, endG3, jgroup, &lndPEGroup, &lndComm );
CHECKIERR( ierr, "Cannot create lnd MPI group and communicator " )
#endif
// we will always have a coupler
MPI_Group couPEGroup;
MPI_Comm couComm;
ierr = create_group_and_comm( startG4, endG4, jgroup, &couPEGroup, &couComm );
CHECKIERR( ierr, "Cannot create cpl MPI group and communicator " )
// now, create the joint communicators atm_coupler, ocn_coupler, lnd_coupler
// for each, we will have to create the group first, then the communicator
// atm_coupler
MPI_Group joinAtmCouGroup;
MPI_Comm atmCouComm;
ierr = create_joint_comm_group( atmPEGroup, couPEGroup, &joinAtmCouGroup, &atmCouComm );
CHECKIERR( ierr, "Cannot create joint atm cou communicator" )
#ifdef ENABLE_ATMOCN_COUPLING
// ocn_coupler
MPI_Group joinOcnCouGroup;
MPI_Comm ocnCouComm;
ierr = create_joint_comm_group( ocnPEGroup, couPEGroup, &joinOcnCouGroup, &ocnCouComm );
CHECKIERR( ierr, "Cannot create joint ocn cou communicator" )
#endif
#ifdef ENABLE_ATMLND_COUPLING
// lnd_coupler
MPI_Group joinLndCouGroup;
MPI_Comm lndCouComm;
ierr = create_joint_comm_group( lndPEGroup, couPEGroup, &joinLndCouGroup, &lndCouComm );
CHECKIERR( ierr, "Cannot create joint ocn cou communicator" )
#endif
ierr = iMOAB_Initialize( argc, argv ); // not really needed anything from argc, argv, yet; maybe we should
CHECKIERR( ierr, "Cannot initialize iMOAB" )
int cmpAtmAppID = -1;
iMOAB_AppID cmpAtmPID = &cmpAtmAppID; // atm
int cplAtmAppID = -1; // -1 means it is not initialized
iMOAB_AppID cplAtmPID = &cplAtmAppID; // atm on coupler PEs
#ifdef ENABLE_ATMOCN_COUPLING
int cmpOcnAppID = -1;
iMOAB_AppID cmpOcnPID = &cmpOcnAppID; // ocn
int cplOcnAppID = -1, cplAtmOcnAppID = -1; // -1 means it is not initialized
iMOAB_AppID cplOcnPID = &cplOcnAppID; // ocn on coupler PEs
iMOAB_AppID cplAtmOcnPID = &cplAtmOcnAppID; // intx atm -ocn on coupler PEs
#endif
#ifdef ENABLE_ATMLND_COUPLING
int cmpLndAppID = -1;
iMOAB_AppID cmpLndPID = &cmpLndAppID; // lnd
int cplLndAppID = -1; //, cplAtmLndAppID=-1;// -1 means it is not initialized
iMOAB_AppID cplLndPID = &cplLndAppID; // land on coupler PEs
// iMOAB_AppID cplAtmLndPID=&cplAtmLndAppID; // intx atm - lnd on coupler PEs not needed anymore
#endif
int cmpPhysAtmID = -1;
iMOAB_AppID cmpPhAtmPID = &cmpPhysAtmID; // phys atm; we do not need to move it to cpl!
if( couComm != MPI_COMM_NULL )
{
MPI_Comm_rank( couComm, &rankInCouComm );
// Register all the applications on the coupler PEs
ierr = iMOAB_RegisterApplication( "ATMX", &couComm, &cplatm,
cplAtmPID ); // atm on coupler pes
CHECKIERR( ierr, "Cannot register ATM over coupler PEs" )
#ifdef ENABLE_ATMOCN_COUPLING
ierr = iMOAB_RegisterApplication( "OCNX", &couComm, &cplocn,
cplOcnPID ); // ocn on coupler pes
CHECKIERR( ierr, "Cannot register OCN over coupler PEs" )
#endif
#ifdef ENABLE_ATMLND_COUPLING
ierr = iMOAB_RegisterApplication( "LNDX", &couComm, &cpllnd,
cplLndPID ); // lnd on coupler pes
CHECKIERR( ierr, "Cannot register LND over coupler PEs" )
#endif
}
if( atmCouComm != MPI_COMM_NULL )
{
ierr = iMOAB_RegisterApplication( "ATM1", &atmComm, &cmpatm, cmpAtmPID );
CHECKIERR( ierr, "Cannot register ATM App" )
}
#ifdef ENABLE_ATMOCN_COUPLING
if( ocnComm != MPI_COMM_NULL )
{
MPI_Comm_rank( ocnComm, &rankInOcnComm );
ierr = iMOAB_RegisterApplication( "OCN1", &ocnComm, &cmpocn, cmpOcnPID );
CHECKIERR( ierr, "Cannot register OCN App" )
}
#endif
// atm
ierr =
setup_component_coupler_meshes( cmpAtmPID, cmpatm, cplAtmPID, cplatm, &atmComm, &atmPEGroup, &couComm,
&couPEGroup, &atmCouComm, atmFilename, readopts, nghlay, repartitioner_scheme );
CHECKIERR( ierr, "Cannot load and migrate atm mesh" )
MPI_Barrier( MPI_COMM_WORLD );
#ifdef ENABLE_ATMOCN_COUPLING
// ocean
ierr =
setup_component_coupler_meshes( cmpOcnPID, cmpocn, cplOcnPID, cplocn, &ocnComm, &ocnPEGroup, &couComm,
&couPEGroup, &ocnCouComm, ocnFilename, readopts, nghlay, repartitioner_scheme );
CHECKIERR( ierr, "Cannot load and migrate ocn mesh" )
#ifdef VERBOSE
if( couComm != MPI_COMM_NULL )
{
char outputFileTgt3[] = "recvOcn2.h5m";
ierr = iMOAB_WriteMesh( cplOcnPID, outputFileTgt3, fileWriteOptions );
CHECKIERR( ierr, "cannot write ocn mesh after receiving" )
}
#endif
#endif
MPI_Barrier( MPI_COMM_WORLD );
// load phys atm mesh, with some data on it already
if( atmComm != MPI_COMM_NULL )
{
ierr = iMOAB_RegisterApplication( "PhysAtm", &atmComm, &cmpPhysAtm, cmpPhAtmPID );
CHECKIERR( ierr, "Cannot register Phys Atm App " )
// load the next component mesh
ierr = iMOAB_LoadMesh( cmpPhAtmPID, atmPhysMesh.c_str(), readoptsPhysAtm.c_str(), &nghlay );
CHECKIERR( ierr, "Cannot load Atm Phys mesh on atm pes" )
int nverts[3], nelem[3];
ierr = iMOAB_GetMeshInfo( cmpPhAtmPID, nverts, nelem, 0, 0, 0 );
CHECKIERR( ierr, "failed to get mesh info" );
printf( "Phys Atm Component Mesh: %d vertices and %d elements\n", nverts[0], nelem[0] );
}
MPI_Barrier( MPI_COMM_WORLD );
#ifdef ENABLE_ATMLND_COUPLING
// land
if( lndComm != MPI_COMM_NULL )
{
ierr = iMOAB_RegisterApplication( "LND1", &lndComm, &cmplnd, cmpLndPID );
CHECKIERR( ierr, "Cannot register LND App " )
}
ierr = setup_component_coupler_meshes( cmpLndPID, cmplnd, cplLndPID, cpllnd, &lndComm, &lndPEGroup, &couComm,
&couPEGroup, &lndCouComm, lndFilename, readoptsPhysAtm, nghlay,
repartitioner_scheme );
if( couComm != MPI_COMM_NULL )
{
char outputFileLnd[] = "recvLnd2.h5m";
ierr = iMOAB_WriteMesh( cplLndPID, outputFileLnd, fileWriteOptions );
CHECKIERR( ierr, "cannot write lnd mesh after receiving" )
}
#endif //
MPI_Barrier( MPI_COMM_WORLD );
#ifdef ENABLE_ATMOCN_COUPLING
if( couComm != MPI_COMM_NULL )
{
// now compute intersection between OCNx and ATMx on coupler PEs
ierr = iMOAB_RegisterApplication( "ATMOCN", &couComm, &atmocnid, cplAtmOcnPID );
CHECKIERR( ierr, "Cannot register ocn_atm intx over coupler pes " )
}
#endif
const char* weights_identifiers[2] = { "scalar", "scalar-pc" };
int disc_orders[3] = { 4, 1, 1 };
const char* disc_methods[3] = { "cgll", "fv", "pcloud" };
const char* dof_tag_names[3] = { "GLOBAL_DOFS", "GLOBAL_ID", "GLOBAL_ID" };
#ifdef ENABLE_ATMOCN_COUPLING
if( couComm != MPI_COMM_NULL )
{
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
CHECKIERR( ierr, "cannot compute intersection" )
}
if( atmCouComm != MPI_COMM_NULL )
{
// 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
CHECKIERR( ierr, "cannot recompute direct coverage graph for ocean" )
}
// need to compute graph between phys atm and atm/ocn intx coverage
if( atmCouComm != MPI_COMM_NULL )
{
int typeA = 2; // point cloud
int typeB = 1; // quads in coverage set
ierr = iMOAB_ComputeCommGraph( cmpPhAtmPID, cplAtmOcnPID, &atmCouComm, &atmPEGroup, &couPEGroup, &typeA, &typeB,
&cmpatm, &atmocnid );
}
#endif
#ifdef ENABLE_ATMLND_COUPLING
// we will compute comm graph between atm phys and land, directly; we do not need any
// intersection later, we will send data from atm phys to land on coupler; then back to land
// comp;
if( atmCouComm != MPI_COMM_NULL )
{
int typeA = 2; // point cloud
int typeB = 2; // point cloud for land on coupler, too
ierr = iMOAB_ComputeCommGraph( cmpPhAtmPID, cplLndPID, &atmCouComm, &atmPEGroup, &couPEGroup, &typeA, &typeB,
&cmpatm, &cpllnd );
}
#endif
MPI_Barrier( MPI_COMM_WORLD );
int fMonotoneTypeID = 0, fVolumetric = 0, fValidate = 1, fNoConserve = 0, fNoBubble = 1, fInverseDistanceMap = 0;
#ifdef ENABLE_ATMOCN_COUPLING
#ifdef VERBOSE
if( couComm != MPI_COMM_NULL )
{
char serialWriteOptions[] = ""; // for writing in serial
std::stringstream outf;
outf << "intxAtmOcn_" << rankInCouComm << ".h5m";
std::string intxfile = outf.str(); // write in serial the intx file, for debugging
ierr = iMOAB_WriteMesh( cplAtmOcnPID, intxfile.c_str(), serialWriteOptions );
CHECKIERR( ierr, "cannot write intx file result" )
}
#endif
if( couComm != MPI_COMM_NULL )
{
PUSH_TIMER( "Compute the projection weights with TempestRemap" )
ierr = iMOAB_ComputeScalarProjectionWeights( cplAtmOcnPID, weights_identifiers[0], disc_methods[0],
&disc_orders[0], disc_methods[1], &disc_orders[1], &fNoBubble,
&fMonotoneTypeID, &fVolumetric, &fInverseDistanceMap, &fNoConserve,
&fValidate, dof_tag_names[0], dof_tag_names[1] );
CHECKIERR( ierr, "cannot compute scalar projection weights" )
POP_TIMER( couComm, rankInCouComm )
}
#endif
MPI_Barrier( MPI_COMM_WORLD );
int tagIndex[2];
int tagTypes[2] = { DENSE_DOUBLE, DENSE_DOUBLE };
int atmCompNDoFs = disc_orders[0] * disc_orders[0], ocnCompNDoFs = 1 /*FV*/;
const char* bottomFields =
"Sa_z:Sa_topo:Sa_u:Sa_v:Sa_tbot:Sa_ptem:Sa_shum:Sa_pbot:Sa_dens:Sa_uovern:Sa_pslv:Sa_co2prog:Sa_co2diag:Faxa_"
"rainc:Faxa_rainl:Faxa_snowc:Faxa_snowl:Faxa_lwdn:Faxa_swndr:Faxa_swvdr:Faxa_swndf:Faxa_swvdf:Faxa_swnet:Faxa_"
"bcphidry:Faxa_bcphodry:Faxa_bcphiwet:Faxa_ocphidry:Faxa_ocphodry:Faxa_ocphiwet:Faxa_dstwet1:Faxa_dstwet2:Faxa_"
"dstwet3:Faxa_dstwet4:Faxa_dstdry1:Faxa_dstdry2:Faxa_dstdry3:Faxa_dstdry4";
const char* bottomFieldsExt =
"Sa_z_ext:Sa_topo_ext:Sa_u_ext:Sa_v_ext:Sa_tbot_ext:Sa_ptem_ext:Sa_shum_ext:Sa_pbot_ext:Sa_dens_ext:Sa_uovern_"
"ext:Sa_pslv_ext:Sa_co2prog_ext:Sa_co2diag_ext:Faxa_rainc_ext:Faxa_rainl_ext:Faxa_snowc_ext:Faxa_snowl_ext:"
"Faxa_lwdn_ext:Faxa_swndr_ext:Faxa_swvdr_ext:Faxa_swndf_ext:Faxa_swvdf_ext:Faxa_swnet_ext:Faxa_bcphidry_ext:"
"Faxa_bcphodry_ext:Faxa_bcphiwet_ext:Faxa_ocphidry_ext:Faxa_ocphodry_ext:Faxa_ocphiwet_ext:Faxa_dstwet1_ext:"
"Faxa_dstwet2_ext:Faxa_dstwet3_ext:Faxa_dstwet4_ext:Faxa_dstdry1_ext:Faxa_dstdry2_ext:Faxa_dstdry3_ext:Faxa_"
"dstdry4_ext";
if( couComm != MPI_COMM_NULL )
{
ierr = iMOAB_DefineTagStorage( cplAtmPID, bottomFieldsExt, &tagTypes[0], &atmCompNDoFs, &tagIndex[0] );
CHECKIERR( ierr, "failed to define the field tags T16_ph:u16_ph:v16_ph " );
#ifdef ENABLE_ATMOCN_COUPLING
ierr = iMOAB_DefineTagStorage( cplOcnPID, bottomFields, &tagTypes[1], &ocnCompNDoFs, &tagIndex[1] );
CHECKIERR( ierr, "failed to define the field tags T_proj:u_proj:v_proj" );
#endif
#ifdef ENABLE_ATMLND_COUPLING
// need to define tag storage for land; will use the same T_proj, u_proj, v_proj name,
// because it will be used to send between point clouds ! use the same ndof and same size as
// ocnCompNDoFs (1) !!
ierr = iMOAB_DefineTagStorage( cplLndPID, bottomFields, &tagTypes[1], &ocnCompNDoFs, &tagIndex[1] );
CHECKIERR( ierr, "failed to define the bottomFields tag on coupler land" );
#endif
}
// need to make sure that the coverage mesh (created during intx method) received the tag that
// need to be projected to target so far, the coverage mesh has only the ids and global dofs;
// need to change the migrate method to accommodate any GLL tag
// now send a tag from original atmosphere (cmpAtmPID) towards migrated coverage mesh
// (cplAtmPID), using the new coverage graph communicator
// make the tag 0, to check we are actually sending needed data
{
if( cplAtmAppID >= 0 )
{
int nverts[3], nelem[3], nblocks[3], nsbc[3], ndbc[3];
/*
* 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 );
CHECKIERR( ierr, "failed to get num primary elems" );
int numAllElem = nelem[2];
std::vector< double > vals;
int storLeng = atmCompNDoFs * numAllElem * 37;
vals.resize( storLeng );
for( int k = 0; k < storLeng; k++ )
vals[k] = 0.;
int eetype = 1;
ierr = iMOAB_SetDoubleTagStorage( cplAtmPID, bottomFieldsExt, &storLeng, &eetype, &vals[0] );
CHECKIERR( ierr, "cannot make tags nul" )
// set the tag to 0
}
}
#ifdef ENABLE_ATMOCN_COUPLING
PUSH_TIMER( "Send/receive data from atm component to coupler in ocn context" )
if( atmComm != MPI_COMM_NULL )
{
// as always, use nonblocking sends
// this is for projection to ocean:
ierr = iMOAB_SendElementTag( cmpPhAtmPID, bottomFields, &atmCouComm, &atmocnid );
CHECKIERR( ierr, "cannot send tag values" )
#ifdef VERBOSE
int is_sender = 1;
int context = atmocnid; // used to identity the parcommgraph in charge
iMOAB_DumpCommGraph( cmpPhAtmPID, &context, &is_sender, "PhysAtmA2OS", );
#endif
}
if( couComm != MPI_COMM_NULL )
{
// receive on atm on coupler pes, that was redistributed according to coverage
ierr = iMOAB_ReceiveElementTag( cplAtmOcnPID, bottomFieldsExt, &atmCouComm, &cmpatm );
CHECKIERR( ierr, "cannot receive tag values" )
#ifdef VERBOSE
int is_sender = 0;
int context = cmpatm; // used to identity the parcommgraph in charge
iMOAB_DumpCommGraph( cplAtmOcnPID, &context, &is_sender, "PhysAtmA2OR", );
#endif
}
POP_TIMER( MPI_COMM_WORLD, rankInGlobalComm )
// we can now free the sender buffers
if( atmComm != MPI_COMM_NULL )
{
ierr = iMOAB_FreeSenderBuffers( cmpPhAtmPID, &atmocnid ); // context is for ocean
CHECKIERR( ierr, "cannot free buffers used to resend atm tag towards the coverage mesh" )
}
/*
#ifdef VERBOSE
if (couComm != MPI_COMM_NULL) {
char outputFileRecvd[] = "recvAtmCoupOcn.h5m";
ierr = iMOAB_WriteMesh(cplAtmPID, outputFileRecvd, fileWriteOptions,
strlen(outputFileRecvd), strlen(fileWriteOptions) );
}
#endif
*/
if( couComm != MPI_COMM_NULL )
{
/* 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 */
PUSH_TIMER( "Apply Scalar projection weights" )
ierr =
iMOAB_ApplyScalarProjectionWeights( cplAtmOcnPID, weights_identifiers[0], bottomFieldsExt, bottomFields );
CHECKIERR( ierr, "failed to compute projection weight application" );
POP_TIMER( couComm, rankInCouComm )
char outputFileTgt[] = "fOcnOnCpl7.h5m";
ierr = iMOAB_WriteMesh( cplOcnPID, outputFileTgt, fileWriteOptions );
}
// send the projected tag back to ocean pes, with send/receive tag
if( ocnComm != MPI_COMM_NULL )
{
int tagIndexIn2;
ierr = iMOAB_DefineTagStorage( cmpOcnPID, bottomFields, &tagTypes[1], &ocnCompNDoFs, &tagIndexIn2 );
CHECKIERR( ierr, "failed to define the field tags for receiving back the atm tags on ocn pes" );
}
// 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( couComm != MPI_COMM_NULL )
{
context_id = cmpocn;
ierr = iMOAB_SendElementTag( cplOcnPID, bottomFields, &ocnCouComm, &context_id );
CHECKIERR( ierr, "cannot send tag values back to ocean pes" )
}
// receive on component 2, ocean
if( ocnComm != MPI_COMM_NULL )
{
context_id = cplocn;
ierr = iMOAB_ReceiveElementTag( cmpOcnPID, bottomFields, &ocnCouComm, &context_id );
CHECKIERR( ierr, "cannot receive tag values from ocean mesh on coupler pes" )
}
MPI_Barrier( MPI_COMM_WORLD );
if( couComm != MPI_COMM_NULL )
{
context_id = cmpocn;
ierr = iMOAB_FreeSenderBuffers( cplOcnPID, &context_id );
}
if( ocnComm != MPI_COMM_NULL )
{
char outputFileOcn[] = "OcnWithProj2.h5m";
ierr = iMOAB_WriteMesh( cmpOcnPID, outputFileOcn, fileWriteOptions );
}
#endif
#ifdef ENABLE_ATMLND_COUPLING
// start land proj:
// we used this to compute
// ierr = iMOAB_ComputeCommGraph(cmpPhAtmPID, cplLndPID, &atmCouComm, &atmPEGroup, &couPEGroup,
// &typeA, &typeB, &cmpatm, &cpllnd);
// end copy
PUSH_TIMER( "Send/receive data from phys comp atm to coupler land, using computed graph" )
if( atmComm != MPI_COMM_NULL )
{
// as always, use nonblocking sends
// this is for projection to land:
ierr = iMOAB_SendElementTag( cmpPhAtmPID, bottomFields, &atmCouComm, &cpllnd );
CHECKIERR( ierr, "cannot send tag values towards cpl on land" )
}
if( couComm != MPI_COMM_NULL )
{
// receive on lnd on coupler pes
ierr = iMOAB_ReceiveElementTag( cplLndPID, bottomFields, &atmCouComm, &cmpatm );
CHECKIERR( ierr, "cannot receive tag values on land on coupler" )
}
POP_TIMER( MPI_COMM_WORLD, rankInGlobalComm )
// we can now free the sender buffers
if( atmComm != MPI_COMM_NULL )
{
ierr = iMOAB_FreeSenderBuffers( cmpPhAtmPID, &cpllnd );
CHECKIERR( ierr, "cannot free buffers used to send atm tag towards the land on coupler" )
}
#ifdef VERBOSE
if( couComm != MPI_COMM_NULL )
{
char outputFileTgtLnd[] = "fLndOnCpl.h5m";
ierr = iMOAB_WriteMesh( cplLndPID, outputFileTgtLnd, fileWriteOptions );
}
#endif
// end land proj
// send the tags back to land pes, from land mesh on coupler pes
// send from cplLndPID to cmpLndPID, using common joint comm
// as always, use nonblocking sends
// original graph
// int context_id = -1;
// the land might not have these tags yet; it should be a different name for land
// in e3sm we do have different names
if( lndComm != MPI_COMM_NULL )
{
int tagIndexIn2;
ierr = iMOAB_DefineTagStorage( cmpLndPID, bottomFields, &tagTypes[1], &ocnCompNDoFs, &tagIndexIn2 );
CHECKIERR( ierr, "failed to define the field tag for receiving back the tag T_proj:u_proj:v_proj on lnd pes" );
}
if( couComm != MPI_COMM_NULL )
{
context_id = cmplnd;
ierr = iMOAB_SendElementTag( cplLndPID, bottomFields, &lndCouComm, &context_id );
CHECKIERR( ierr, "cannot send tag values back to land pes" )
}
// receive on component 3, land
if( lndComm != MPI_COMM_NULL )
{
context_id = cpllnd; // receive from coupler on land
ierr = iMOAB_ReceiveElementTag( cmpLndPID, bottomFields, &lndCouComm, &context_id );
CHECKIERR( ierr, "cannot receive tag values from land mesh on coupler pes" )
}
MPI_Barrier( MPI_COMM_WORLD );
if( couComm != MPI_COMM_NULL )
{
context_id = cmplnd;
ierr = iMOAB_FreeSenderBuffers( cplLndPID, &context_id );
}
if( lndComm != MPI_COMM_NULL )
{
char outputFileLnd[] = "LndWithProj2.h5m";
ierr = iMOAB_WriteMesh( cmpLndPID, outputFileLnd, fileWriteOptions );
}
#endif // ENABLE_ATMLND_COUPLING
#ifdef ENABLE_ATMOCN_COUPLING
if( couComm != MPI_COMM_NULL )
{
ierr = iMOAB_DeregisterApplication( cplAtmOcnPID );
CHECKIERR( ierr, "cannot deregister app intx AO" )
}
#endif
#ifdef ENABLE_ATMLND_COUPLING
if( lndComm != MPI_COMM_NULL )
{
ierr = iMOAB_DeregisterApplication( cmpLndPID );
CHECKIERR( ierr, "cannot deregister app LND1" )
}
#endif // ENABLE_ATMLND_COUPLING
if( atmComm != MPI_COMM_NULL )
{
ierr = iMOAB_DeregisterApplication( cmpPhAtmPID );
CHECKIERR( ierr, "cannot deregister app PhysAtm " )
}
#ifdef ENABLE_ATMOCN_COUPLING
if( ocnComm != MPI_COMM_NULL )
{
ierr = iMOAB_DeregisterApplication( cmpOcnPID );
CHECKIERR( ierr, "cannot deregister app OCN1" )
}
#endif // ENABLE_ATMOCN_COUPLING
if( atmComm != MPI_COMM_NULL )
{
ierr = iMOAB_DeregisterApplication( cmpAtmPID );
CHECKIERR( ierr, "cannot deregister app ATM1" )
}
#ifdef ENABLE_ATMLND_COUPLING
if( couComm != MPI_COMM_NULL )
{
ierr = iMOAB_DeregisterApplication( cplLndPID );
CHECKIERR( ierr, "cannot deregister app LNDX" )
}
#endif // ENABLE_ATMLND_COUPLING
#ifdef ENABLE_ATMOCN_COUPLING
if( couComm != MPI_COMM_NULL )
{
ierr = iMOAB_DeregisterApplication( cplOcnPID );
CHECKIERR( ierr, "cannot deregister app OCNX" )
}
#endif // ENABLE_ATMOCN_COUPLING
if( couComm != MPI_COMM_NULL )
{
ierr = iMOAB_DeregisterApplication( cplAtmPID );
CHECKIERR( ierr, "cannot deregister app ATMX" )
}
//#endif
ierr = iMOAB_Finalize();
CHECKIERR( ierr, "did not finalize iMOAB" )
// free atm coupler group and comm
if( MPI_COMM_NULL != atmCouComm ) MPI_Comm_free( &atmCouComm );
MPI_Group_free( &joinAtmCouGroup );
if( MPI_COMM_NULL != atmComm ) MPI_Comm_free( &atmComm );
#ifdef ENABLE_ATMOCN_COUPLING
if( MPI_COMM_NULL != ocnComm ) MPI_Comm_free( &ocnComm );
// free ocn - coupler group and comm
if( MPI_COMM_NULL != ocnCouComm ) MPI_Comm_free( &ocnCouComm );
MPI_Group_free( &joinOcnCouGroup );
#endif
#ifdef ENABLE_ATMLND_COUPLING
if( MPI_COMM_NULL != lndComm ) MPI_Comm_free( &lndComm );
// free land - coupler group and comm
if( MPI_COMM_NULL != lndCouComm ) MPI_Comm_free( &lndCouComm );
MPI_Group_free( &joinLndCouGroup );
#endif
if( MPI_COMM_NULL != couComm ) MPI_Comm_free( &couComm );
MPI_Group_free( &atmPEGroup );
#ifdef ENABLE_ATMOCN_COUPLING
MPI_Group_free( &ocnPEGroup );
#endif
#ifdef ENABLE_ATMLND_COUPLING
MPI_Group_free( &lndPEGroup );
#endif
MPI_Group_free( &couPEGroup );
MPI_Group_free( &jgroup );
MPI_Finalize();
// endif #if 0
return 0;
}
1.7.6.1