16 #error mbtempest tool requires MPI configuration
37 #ifndef MOAB_HAVE_TEMPESTREMAP
38 #error The climate coupler test example requires MOAB configuration with TempestRemap
41 #define ENABLE_ATMOCN_COUPLING
42 #define ENABLE_ATMLND_COUPLING
44 #if( !defined( ENABLE_ATMOCN_COUPLING ) && !defined( ENABLE_ATMLND_COUPLING ) )
45 #error Enable either OCN (ENABLE_ATMOCN_COUPLING) and/or LND (ENABLE_ATMLND_COUPLING) for coupling
48 int main(
int argc,
char* argv[] )
53 std::string
readopts(
"PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS" );
54 std::string readoptsLnd(
"PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" );
61 int repartitioner_scheme = 0;
62 #ifdef MOAB_HAVE_ZOLTAN
63 repartitioner_scheme = 2;
66 MPI_Init( &argc, &argv );
72 std::string
atmFilename = TestDir +
"unittest/wholeATM_T.h5m";
87 #ifdef ENABLE_ATMOCN_COUPLING
88 std::string ocnFilename = TestDir +
"unittest/recMeshOcn.h5m";
89 std::string baseline = TestDir +
"unittest/baseline1.txt";
90 int rankInOcnComm = -1;
91 int cmpocn = 17, cplocn = 18,
95 #ifdef ENABLE_ATMLND_COUPLING
96 std::string lndFilename = TestDir +
"unittest/wholeLnd.h5m";
97 int cpllnd = 10, cmplnd = 9,
101 int rankInCouComm = -1;
114 opts.
addOpt< std::string >(
"atmosphere,t",
"atm mesh filename (source)", &
atmFilename );
115 #ifdef ENABLE_ATMOCN_COUPLING
116 opts.
addOpt< std::string >(
"ocean,m",
"ocean mesh filename (target)", &ocnFilename );
118 #ifdef ENABLE_ATMLND_COUPLING
119 opts.
addOpt< std::string >(
"land,l",
"land mesh filename (target)", &lndFilename );
121 opts.
addOpt<
int >(
"startAtm,a",
"start task for atmosphere layout", &
startG1 );
122 opts.
addOpt<
int >(
"endAtm,b",
"end task for atmosphere layout", &
endG1 );
123 #ifdef ENABLE_ATMOCN_COUPLING
124 opts.
addOpt<
int >(
"startOcn,c",
"start task for ocean layout", &
startG2 );
125 opts.
addOpt<
int >(
"endOcn,d",
"end task for ocean layout", &
endG2 );
127 #ifdef ENABLE_ATMLND_COUPLING
128 opts.
addOpt<
int >(
"startLnd,e",
"start task for land layout", &startG3 );
129 opts.
addOpt<
int >(
"endLnd,f",
"end task for land layout", &endG3 );
132 opts.
addOpt<
int >(
"startCoupler,g",
"start task for coupler layout", &startG4 );
133 opts.
addOpt<
int >(
"endCoupler,j",
"end task for coupler layout", &endG4 );
135 opts.
addOpt<
int >(
"partitioning,p",
"partitioning option for migration", &repartitioner_scheme );
138 opts.
addOpt<
int >(
"iterations,n",
"number of iterations for coupler", &n );
140 bool no_regression_test =
false;
141 opts.
addOpt<
void >(
"no_regression,r",
"do not do regression test against baseline 1", &no_regression_test );
149 #ifdef ENABLE_ATMOCN_COUPLING
150 "\n ocn file: " << ocnFilename <<
"\n on tasks : " <<
startG2 <<
":" <<
endG2 <<
152 #ifdef ENABLE_ATMLND_COUPLING
153 "\n lnd file: " << lndFilename <<
"\n on tasks : " << startG3 <<
":" << endG3 <<
155 "\n partitioning (0 trivial, 1 graph, 2 geometry) " << repartitioner_scheme <<
"\n ";
163 MPI_Group atmPEGroup;
166 CHECKIERR(
ierr,
"Cannot create atm MPI group and communicator " )
168 #ifdef ENABLE_ATMOCN_COUPLING
169 MPI_Group ocnPEGroup;
172 CHECKIERR(
ierr,
"Cannot create ocn MPI group and communicator " )
175 #ifdef ENABLE_ATMLND_COUPLING
176 MPI_Group lndPEGroup;
179 CHECKIERR(
ierr,
"Cannot create lnd MPI group and communicator " )
183 MPI_Group couPEGroup;
186 CHECKIERR(
ierr,
"Cannot create cpl MPI group and communicator " )
189 MPI_Group joinAtmCouGroup;
192 CHECKIERR(
ierr,
"Cannot create joint atm cou communicator" )
194 #ifdef ENABLE_ATMOCN_COUPLING
196 MPI_Group joinOcnCouGroup;
199 CHECKIERR(
ierr,
"Cannot create joint ocn cou communicator" )
202 #ifdef ENABLE_ATMLND_COUPLING
204 MPI_Group joinLndCouGroup;
207 CHECKIERR(
ierr,
"Cannot create joint ocn cou communicator" )
213 int cmpAtmAppID = -1;
215 int cplAtmAppID = -1;
217 #ifdef ENABLE_ATMOCN_COUPLING
218 int cmpOcnAppID = -1;
220 int cplOcnAppID = -1, cplAtmOcnAppID = -1;
225 #ifdef ENABLE_ATMLND_COUPLING
226 int cmpLndAppID = -1;
228 int cplLndAppID = -1, cplAtmLndAppID = -1;
233 if( couComm != MPI_COMM_NULL )
235 MPI_Comm_rank( couComm, &rankInCouComm );
240 #ifdef ENABLE_ATMOCN_COUPLING
245 #ifdef ENABLE_ATMLND_COUPLING
252 if( atmComm != MPI_COMM_NULL )
259 #ifdef ENABLE_ATMOCN_COUPLING
260 if( ocnComm != MPI_COMM_NULL )
262 MPI_Comm_rank( ocnComm, &rankInOcnComm );
274 int global_atm[2] = { 0, 0 };
275 if( atmComm != MPI_COMM_NULL )
280 std::cout <<
" atm mesh: vertices: " << global_atm[0] <<
" cells: " << global_atm[1] <<
"\n";
283 int global_atm_red[2] = { 0, 0 };
284 MPI_Allreduce( global_atm, global_atm_red, 2, MPI_INT, MPI_MAX,
MPI_COMM_WORLD );
285 if( couComm != MPI_COMM_NULL )
287 int global_atm_cou[2];
292 std::cout <<
" atm mesh on coupler: vertices: " << global_atm[0] <<
" cells: " << global_atm[1] <<
"\n";
296 if( couComm != MPI_COMM_NULL && 1 == n )
298 char outputFileTgt3[] =
"recvAtm.h5m";
305 if( atmComm != MPI_COMM_NULL )
310 iMOAB_DumpCommGraph( cmpAtmPID, &
context, &is_sender,
"AtmMigS" );
312 if( couComm != MPI_COMM_NULL )
316 iMOAB_DumpCommGraph( cplAtmPID, &
context, &is_sender,
"AtmMigR" );
321 #ifdef ENABLE_ATMOCN_COUPLING
325 &couPEGroup, &ocnCouComm, ocnFilename,
readopts,
nghlay, repartitioner_scheme );
331 if( couComm != MPI_COMM_NULL && 1 == n )
333 char outputFileTgt3[] =
"recvOcn.h5m";
340 #ifdef ENABLE_ATMLND_COUPLING
342 if( lndComm != MPI_COMM_NULL )
348 &couPEGroup, &lndCouComm, lndFilename, readoptsLnd,
nghlay,
349 repartitioner_scheme );
351 if( couComm != MPI_COMM_NULL && 1 == n )
353 char outputFileLnd[] =
"recvLnd.h5m";
362 #ifdef ENABLE_ATMOCN_COUPLING
363 if( couComm != MPI_COMM_NULL )
367 CHECKIERR(
ierr,
"Cannot register ocn_atm intx over coupler pes " )
370 #ifdef ENABLE_ATMLND_COUPLING
371 if( couComm != MPI_COMM_NULL )
375 CHECKIERR(
ierr,
"Cannot register ocn_atm intx over coupler pes " )
379 int disc_orders[3] = { 4, 1, 1 };
380 const std::string weights_identifiers[2] = {
"scalar",
"scalar-pc" };
381 const std::string disc_methods[3] = {
"cgll",
"fv",
"pcloud" };
382 const std::string dof_tag_names[3] = {
"GLOBAL_DOFS",
"GLOBAL_ID",
"GLOBAL_ID" };
383 #ifdef ENABLE_ATMOCN_COUPLING
384 if( couComm != MPI_COMM_NULL )
386 PUSH_TIMER(
"Compute ATM-OCN mesh intersection" )
387 ierr = iMOAB_ComputeMeshIntersectionOnSphere( cplAtmPID, cplOcnPID, cplAtmOcnPID );
394 char prefix[] =
"intx_atmocn";
400 if( atmCouComm != MPI_COMM_NULL )
410 PUSH_TIMER(
"Compute OCN coverage graph for ATM mesh" )
411 ierr = iMOAB_CoverageGraph( &atmCouComm, cmpAtmPID, cplAtmPID, cplAtmOcnPID, &
cmpatm, &cplatm,
413 CHECKIERR(
ierr,
"cannot recompute direct coverage graph for ocean" )
418 #ifdef ENABLE_ATMLND_COUPLING
419 if( couComm != MPI_COMM_NULL )
421 PUSH_TIMER(
"Compute ATM-LND mesh intersection" )
422 ierr = iMOAB_ComputePointDoFIntersection( cplAtmPID, cplLndPID, cplAtmLndPID );
426 if( atmCouComm != MPI_COMM_NULL )
435 PUSH_TIMER(
"Compute LND coverage graph for ATM mesh" )
436 ierr = iMOAB_CoverageGraph( &atmCouComm, cmpAtmPID, cplAtmPID, cplAtmLndPID, &
cmpatm, &cplatm,
438 CHECKIERR(
ierr,
"cannot recompute direct coverage graph for land" )
445 int fMonotoneTypeID = 0, fVolumetric = 0, fValidate = 0, fNoConserve = 0, fNoBubble = 1, fInverseDistanceMap = 0;
447 #ifdef ENABLE_ATMOCN_COUPLING
449 if( couComm != MPI_COMM_NULL && 1 == n )
451 char serialWriteOptions[] =
"";
452 std::stringstream outf;
453 outf <<
"intxAtmOcn_" << rankInCouComm <<
".h5m";
454 std::string intxfile = outf.str();
460 if( couComm != MPI_COMM_NULL )
462 PUSH_TIMER(
"Compute the projection weights with TempestRemap" )
463 ierr = iMOAB_ComputeScalarProjectionWeights( cplAtmOcnPID, weights_identifiers[0].c_str(),
464 disc_methods[0].c_str(), &disc_orders[0], disc_methods[1].c_str(),
465 &disc_orders[1],
nullptr, &fNoBubble, &fMonotoneTypeID,
466 &fVolumetric, &fInverseDistanceMap, &fNoConserve, &fValidate,
467 dof_tag_names[0].c_str(), dof_tag_names[1].c_str() );
468 CHECKIERR(
ierr,
"cannot compute scalar projection weights" )
472 #ifdef MOAB_HAVE_NETCDF
474 const std::string atmocn_map_file_name =
"atm_ocn_map.nc";
475 ierr = iMOAB_WriteMappingWeightsToFile( cplAtmOcnPID, weights_identifiers[0].c_str(),
476 atmocn_map_file_name.c_str() );
479 const std::string intx_from_file_identifier =
"map-from-file";
481 int dummy_rowcol = -1;
483 ierr = iMOAB_LoadMappingWeightsFromFile( cplAtmOcnPID, &dummyCpl, &dummy_rowcol, &dummyType,
484 intx_from_file_identifier.c_str(), atmocn_map_file_name.c_str() );
494 #ifdef ENABLE_ATMLND_COUPLING
495 if( couComm != MPI_COMM_NULL )
499 PUSH_TIMER(
"Compute ATM-LND remapping weights" )
500 ierr = iMOAB_ComputeScalarProjectionWeights( cplAtmLndPID, weights_identifiers[1].c_str(),
501 disc_methods[0].c_str(), &disc_orders[0], disc_methods[2].c_str(),
502 &disc_orders[2],
nullptr, &fNoBubble, &fMonotoneTypeID,
503 &fVolumetric, &fInverseDistanceMap, &fNoConserve, &fValidate,
504 dof_tag_names[0].c_str(), dof_tag_names[2].c_str() );
505 CHECKIERR(
ierr,
"failed to compute remapping projection weights for ATM-LND scalar "
506 "non-conservative field" );
528 int atmCompNDoFs = disc_orders[0] * disc_orders[0], ocnCompNDoFs = 1 ;
530 const char* bottomFields =
"a2oTbot:a2oUbot:a2oVbot";
531 const char* bottomProjectedFields =
"a2oTbot_proj:a2oUbot_proj:a2oVbot_proj";
533 if( couComm != MPI_COMM_NULL )
537 #ifdef ENABLE_ATMOCN_COUPLING
539 CHECKIERR(
ierr,
"failed to define the field tag a2oTbot_proj" );
551 if( cplAtmAppID >= 0 )
553 int nverts[3], nelem[3], nblocks[3], nsbc[3], ndbc[3];
563 int numAllElem = nelem[2];
564 std::vector< double > vals;
565 int storLeng = atmCompNDoFs * numAllElem * 3;
568 vals.resize( storLeng );
569 for(
int k = 0; k < storLeng; k++ )
578 const char* concat_fieldname =
"a2oTbot:a2oUbot:a2oVbot";
579 const char* concat_fieldnameT =
"a2oTbot_proj:a2oUbot_proj:a2oVbot_proj";
582 for(
int iters = 0; iters < n; iters++ )
584 #ifdef ENABLE_ATMOCN_COUPLING
585 PUSH_TIMER(
"Send/receive data from atm component to coupler in ocn context" )
586 if( atmComm != MPI_COMM_NULL )
590 ierr = iMOAB_SendElementTag( cmpAtmPID,
"a2oTbot:a2oUbot:a2oVbot", &atmCouComm, &cplocn );
595 iMOAB_DumpCommGraph( cmpAtmPID, &
context, &is_sender,
"AtmCovOcnS" );
598 if( couComm != MPI_COMM_NULL )
601 ierr = iMOAB_ReceiveElementTag( cplAtmPID,
"a2oTbot:a2oUbot:a2oVbot", &atmCouComm, &cplocn );
606 iMOAB_DumpCommGraph( cmpAtmPID, &
context, &is_sender,
"AtmCovOcnR" );
612 if( atmComm != MPI_COMM_NULL )
614 ierr = iMOAB_FreeSenderBuffers( cmpAtmPID, &cplocn );
615 CHECKIERR(
ierr,
"cannot free buffers used to resend atm tag towards the coverage mesh" )
618 if( couComm != MPI_COMM_NULL && 1 == n )
621 char outputFileRecvd[] =
"recvAtmCoupOcn.h5m";
623 CHECKIERR(
ierr,
"could not write recvAtmCoupOcn.h5m to disk" )
627 if( couComm != MPI_COMM_NULL )
631 PUSH_TIMER(
"Apply Scalar projection weights" )
632 ierr = iMOAB_ApplyScalarProjectionWeights( cplAtmOcnPID, weights_identifiers[0].c_str(), concat_fieldname,
634 CHECKIERR(
ierr,
"failed to compute projection weight application" );
638 char outputFileTgt[] =
"fOcnOnCpl.h5m";
645 if( ocnComm != MPI_COMM_NULL )
650 CHECKIERR(
ierr,
"failed to define the field tag for receiving back the tags "
651 "a2oTbot_proj, a2oUbot_proj, a2oVbot_proj on ocn pes" );
657 if( couComm != MPI_COMM_NULL )
662 iMOAB_SendElementTag( cplOcnPID,
"a2oTbot_proj:a2oUbot_proj:a2oVbot_proj", &ocnCouComm, &context_id );
663 CHECKIERR(
ierr,
"cannot send tag values back to ocean pes" )
667 if( ocnComm != MPI_COMM_NULL )
670 ierr = iMOAB_ReceiveElementTag( cmpOcnPID,
"a2oTbot_proj:a2oUbot_proj:a2oVbot_proj", &ocnCouComm,
672 CHECKIERR(
ierr,
"cannot receive tag values from ocean mesh on coupler pes" )
677 if( couComm != MPI_COMM_NULL )
680 ierr = iMOAB_FreeSenderBuffers( cplOcnPID, &context_id );
681 CHECKIERR(
ierr,
"cannot free send/receive buffers for OCN context" )
683 if( ocnComm != MPI_COMM_NULL && 1 == n )
685 char outputFileOcn[] =
"OcnWithProj.h5m";
689 if( !no_regression_test )
694 int nverts[3], nelem[3];
697 std::vector< int > gidElems;
698 gidElems.resize( nelem[2] );
699 std::vector< double > tempElems;
700 tempElems.resize( nelem[2] );
702 const std::string GidStr =
"GLOBAL_ID";
713 check_baseline_file( baseline, gidElems, tempElems, 1.e-9, err_code );
715 std::cout <<
" passed baseline test atm2ocn on ocean task " << rankInOcnComm <<
"\n";
720 #ifdef ENABLE_ATMLND_COUPLING
722 PUSH_TIMER(
"Send/receive data from component atm to coupler, in land context" )
723 if( atmComm != MPI_COMM_NULL )
727 ierr = iMOAB_SendElementTag( cmpAtmPID,
"a2oTbot:a2oUbot:a2oVbot", &atmCouComm, &cpllnd );
730 if( couComm != MPI_COMM_NULL )
734 ierr = iMOAB_ReceiveElementTag( cplAtmPID,
"a2oTbot:a2oUbot:a2oVbot", &atmCouComm, &cpllnd );
740 if( atmComm != MPI_COMM_NULL )
742 ierr = iMOAB_FreeSenderBuffers( cmpAtmPID, &cpllnd );
743 CHECKIERR(
ierr,
"cannot free buffers used to resend atm tag towards the coverage mesh "
747 if( couComm != MPI_COMM_NULL && 1 == n )
749 char outputFileRecvd[] =
"recvAtmCoupLnd.h5m";
751 CHECKIERR(
ierr,
"could not write recvAtmCoupLnd.h5m to disk" )
757 if( couComm != MPI_COMM_NULL )
759 PUSH_TIMER(
"Apply Scalar projection weights for land" )
760 ierr = iMOAB_ApplyScalarProjectionWeights( cplAtmLndPID, weights_identifiers[1].c_str(), concat_fieldname,
762 CHECKIERR(
ierr,
"failed to compute projection weight application" );
767 if( couComm != MPI_COMM_NULL && 1 == n )
769 char outputFileTgtLnd[] =
"fLndOnCpl.h5m";
783 if( lndComm != MPI_COMM_NULL )
788 CHECKIERR(
ierr,
"failed to define the field tag for receiving back the tag "
789 "a2oTbot_proj, a2oUbot_proj, a2oVbot_proj on lnd pes" );
791 if( couComm != MPI_COMM_NULL )
795 iMOAB_SendElementTag( cplLndPID,
"a2oTbot_proj:a2oUbot_proj:a2oVbot_proj", &lndCouComm, &context_id );
799 if( lndComm != MPI_COMM_NULL )
802 ierr = iMOAB_ReceiveElementTag( cmpLndPID,
"a2oTbot_proj:a2oUbot_proj:a2oVbot_proj", &lndCouComm,
804 CHECKIERR(
ierr,
"cannot receive tag values from land mesh on coupler pes" )
808 if( couComm != MPI_COMM_NULL )
811 ierr = iMOAB_FreeSenderBuffers( cplLndPID, &context_id );
812 CHECKIERR(
ierr,
"cannot free buffers used to resend land tag towards the coverage mesh "
815 if( lndComm != MPI_COMM_NULL && 1 == n )
817 char outputFileLnd[] =
"LndWithProj.h5m";
824 #ifdef ENABLE_ATMLND_COUPLING
825 if( lndComm != MPI_COMM_NULL )
832 #ifdef ENABLE_ATMOCN_COUPLING
833 if( couComm != MPI_COMM_NULL )
838 if( ocnComm != MPI_COMM_NULL )
845 if( atmComm != MPI_COMM_NULL )
851 #ifdef ENABLE_ATMLND_COUPLING
852 if( couComm != MPI_COMM_NULL )
859 #ifdef ENABLE_ATMOCN_COUPLING
860 if( couComm != MPI_COMM_NULL )
867 if( couComm != MPI_COMM_NULL )
878 if( MPI_COMM_NULL != atmCouComm ) MPI_Comm_free( &atmCouComm );
879 MPI_Group_free( &joinAtmCouGroup );
880 if( MPI_COMM_NULL != atmComm ) MPI_Comm_free( &atmComm );
882 #ifdef ENABLE_ATMOCN_COUPLING
883 if( MPI_COMM_NULL != ocnComm ) MPI_Comm_free( &ocnComm );
885 if( MPI_COMM_NULL != ocnCouComm ) MPI_Comm_free( &ocnCouComm );
886 MPI_Group_free( &joinOcnCouGroup );
889 #ifdef ENABLE_ATMLND_COUPLING
890 if( MPI_COMM_NULL != lndComm ) MPI_Comm_free( &lndComm );
892 if( MPI_COMM_NULL != lndCouComm ) MPI_Comm_free( &lndCouComm );
893 MPI_Group_free( &joinLndCouGroup );
896 if( MPI_COMM_NULL != couComm ) MPI_Comm_free( &couComm );
898 MPI_Group_free( &atmPEGroup );
899 #ifdef ENABLE_ATMOCN_COUPLING
900 MPI_Group_free( &ocnPEGroup );
902 #ifdef ENABLE_ATMLND_COUPLING
903 MPI_Group_free( &lndPEGroup );
905 MPI_Group_free( &couPEGroup );
906 MPI_Group_free( &
jgroup );