MOAB: Mesh Oriented datABase  (version 5.4.0)
imoab_remap.cpp
Go to the documentation of this file.
00001 #include "moab/MOABConfig.h"
00002 
00003 #ifdef MOAB_HAVE_MPI
00004 #include "moab_mpi.h"
00005 #endif
00006 #include "moab/iMOAB.h"
00007 #include "TestUtil.hpp"
00008 #include "moab/CpuTimer.hpp"
00009 #include "moab/ProgOptions.hpp"
00010 
00011 // for malloc, free:
00012 #include <iostream>
00013 #include <iomanip>
00014 
00015 #define CHECKIERR( ierr, message ) \
00016     if( 0 != ( ierr ) )            \
00017     {                              \
00018         printf( "%s", message );   \
00019         return 1;                  \
00020     }
00021 
00022 #define ENABLE_ATMLND_COUPLING
00023 // this test will be run in serial only
00024 int main( int argc, char* argv[] )
00025 {
00026     ErrCode ierr;
00027     std::string atmFilename, ocnFilename, lndFilename, mapFilename;
00028 
00029     atmFilename = TestDir + "unittest/wholeATM_T.h5m";
00030     ocnFilename = TestDir + "unittest/recMeshOcn.h5m";
00031 #ifdef ENABLE_ATMLND_COUPLING
00032     lndFilename = TestDir + "unittest/wholeLnd.h5m";
00033 #endif
00034     // mapFilename = TestDir + "unittest/outCS5ICOD5_map.nc";
00035 
00036     ProgOptions opts;
00037     opts.addOpt< std::string >( "atmosphere,t", "atm mesh filename (source)", &atmFilename );
00038     opts.addOpt< std::string >( "ocean,m", "ocean mesh filename (target)", &ocnFilename );
00039 #ifdef ENABLE_ATMLND_COUPLING
00040     opts.addOpt< std::string >( "land,l", "land mesh filename (target)", &lndFilename );
00041 #endif
00042     bool gen_baseline = false;
00043     opts.addOpt< void >( "genbase,g", "generate baseline 1", &gen_baseline );
00044 
00045     bool no_test_against_baseline = false;
00046     opts.addOpt< void >( "no_testbase,t", "do not test against baseline 1", &no_test_against_baseline );
00047     std::string baseline = TestDir + "unittest/baseline1.txt";
00048 
00049     opts.parseCommandLine( argc, argv );
00050 
00051     {
00052         std::cout << " atm file: " << atmFilename << "\n ocn file: " << ocnFilename
00053 #ifdef ENABLE_ATMLND_COUPLING
00054                   << "\n lnd file: " << lndFilename
00055 #endif
00056                   << "\n";
00057     }
00058 
00059 #ifdef MOAB_HAVE_MPI
00060     MPI_Init( &argc, &argv );
00061     MPI_Comm comm = MPI_COMM_WORLD;
00062 #endif
00063     /*
00064      * MOAB needs to be initialized; A MOAB instance will be created, and will be used by each
00065      * application in this framework. There is no parallel context yet.
00066      */
00067     ierr = iMOAB_Initialize( argc, argv );
00068     CHECKIERR( ierr, "failed to initialize MOAB" );
00069 
00070     int atmAppID, ocnAppID, lndAppID, atmocnAppID, atmlndAppID, lndatmAppID;
00071     iMOAB_AppID atmPID    = &atmAppID;
00072     iMOAB_AppID ocnPID    = &ocnAppID;
00073     iMOAB_AppID lndPID    = &lndAppID;
00074     iMOAB_AppID atmocnPID = &atmocnAppID;
00075     iMOAB_AppID atmlndPID = &atmlndAppID;
00076     iMOAB_AppID lndatmPID = &lndatmAppID;
00077 
00078     /*
00079      * Each application has to be registered once. A mesh set and a parallel communicator will be
00080      * associated with each application. A unique application id will be returned, and will be used
00081      * for all future mesh operations/queries.
00082      */
00083     int atmCompID = 10, ocnCompID = 20, lndCompID = 30, atmocnCompID = 100, atmlndCompID = 102, lndatmCompID = 103;
00084     ierr = iMOAB_RegisterApplication( "ATM-APP",
00085 #ifdef MOAB_HAVE_MPI
00086                                       &comm,
00087 #endif
00088                                       &atmCompID, atmPID );
00089     CHECKIERR( ierr, "failed to register application1" );
00090 
00091     ierr = iMOAB_RegisterApplication( "OCN-APP",
00092 #ifdef MOAB_HAVE_MPI
00093                                       &comm,
00094 #endif
00095                                       &ocnCompID, ocnPID );
00096     CHECKIERR( ierr, "failed to register application2" );
00097 #ifdef ENABLE_ATMLND_COUPLING
00098     ierr = iMOAB_RegisterApplication( "LND-APP",
00099 #ifdef MOAB_HAVE_MPI
00100                                       &comm,
00101 #endif
00102                                       &lndCompID, lndPID );
00103     CHECKIERR( ierr, "failed to register application3" );
00104 #endif
00105     ierr = iMOAB_RegisterApplication( "ATM-OCN-CPL",
00106 #ifdef MOAB_HAVE_MPI
00107                                       &comm,
00108 #endif
00109                                       &atmocnCompID, atmocnPID );
00110     CHECKIERR( ierr, "failed to register application4" );
00111 #ifdef ENABLE_ATMLND_COUPLING
00112     ierr = iMOAB_RegisterApplication( "ATM-LND-CPL",
00113 #ifdef MOAB_HAVE_MPI
00114                                       &comm,
00115 #endif
00116                                       &atmlndCompID, atmlndPID );
00117     CHECKIERR( ierr, "failed to register application5" );
00118 
00119     ierr = iMOAB_RegisterApplication( "LND-ATM-CPL",
00120 #ifdef MOAB_HAVE_MPI
00121                                       &comm,
00122 #endif
00123                                       &lndatmCompID, lndatmPID );
00124     CHECKIERR( ierr, "failed to register application5" );
00125 #endif
00126     const char* read_opts = "";
00127     int num_ghost_layers  = 0;
00128     /*
00129      * Loading the mesh is a parallel IO operation. Ghost layers can be exchanged too, and default
00130      * MOAB sets are augmented with ghost elements. By convention, blocks correspond to MATERIAL_SET
00131      * sets, side sets with NEUMANN_SET sets, node sets with DIRICHLET_SET sets. Each element and
00132      * vertex entity should have a GLOBAL ID tag in the file, which will be available for visible
00133      * entities
00134      */
00135     ierr = iMOAB_LoadMesh( atmPID, atmFilename.c_str(), read_opts, &num_ghost_layers );
00136     CHECKIERR( ierr, "failed to load mesh" );
00137     ierr = iMOAB_LoadMesh( ocnPID, ocnFilename.c_str(), read_opts, &num_ghost_layers );
00138     CHECKIERR( ierr, "failed to load mesh" );
00139 #ifdef ENABLE_ATMLND_COUPLING
00140     ierr = iMOAB_LoadMesh( lndPID, lndFilename.c_str(), read_opts, &num_ghost_layers );
00141     CHECKIERR( ierr, "failed to load mesh" );
00142 #endif
00143 
00144     int nverts[3], nelem[3];
00145     /*
00146      * Each process in the communicator will have access to a local mesh instance, which will
00147      * contain the original cells in the local partition and ghost entities. Number of vertices,
00148      * primary cells, visible blocks, number of sidesets and nodesets boundary conditions will be
00149      * returned in size 3 arrays, for local, ghost and total numbers.
00150      */
00151     ierr = iMOAB_GetMeshInfo( atmPID, nverts, nelem, 0, 0, 0 );
00152     CHECKIERR( ierr, "failed to get mesh info" );
00153     printf( "Atmosphere Component Mesh: %d vertices and %d elements\n", nverts[0], nelem[0] );
00154 
00155     ierr = iMOAB_GetMeshInfo( ocnPID, nverts, nelem, 0, 0, 0 );
00156     CHECKIERR( ierr, "failed to get mesh info" );
00157     printf( "Ocean Component Mesh: %d vertices and %d elements\n", nverts[0], nelem[0] );
00158 #ifdef ENABLE_ATMLND_COUPLING
00159     ierr = iMOAB_GetMeshInfo( lndPID, nverts, nelem, 0, 0, 0 );
00160     CHECKIERR( ierr, "failed to get mesh info" );
00161     printf( "Land Component Mesh: %d vertices and %d elements\n", nverts[0], nelem[0] );
00162 #endif
00163     /*
00164      * The 2 tags used in this example exist in the file, already.
00165      * If a user needs a new tag, it can be defined with the same call as this one
00166      * this method, iMOAB_DefineTagStorage, will return a local index for the tag.
00167      * The name of the tag is case sensitive.
00168      * This method is collective.
00169      */
00170     // int disc_orders[2] = {1, 1};
00171     // const char* disc_methods[2] = {"fv", "fv"};
00172     // const char* dof_tag_names[2] = {"GLOBAL_ID", "GLOBAL_ID"};
00173     int disc_orders[3]  = { 4, 1, 1 };
00174     int fMonotoneTypeID = 0, fVolumetric = 0, fValidate = 1, fNoConserve = 0, fNoBubble = 1, fInverseDistanceMap = 0;
00175 
00176     const std::string disc_methods[3]        = { "cgll", "fv", "pcloud" };
00177     const std::string dof_tag_names[3]       = { "GLOBAL_DOFS", "GLOBAL_ID", "GLOBAL_ID" };
00178     const std::string weights_identifiers[3] = { "scalar", "scalar_pointcloud", "scalar_conservative" };
00179 
00180     const std::string bottomTempField            = "a2oTbot";
00181     const std::string bottomTempFieldATM         = "a2oTbotATM";
00182     const std::string bottomTempProjectedField   = "a2oTbot_proj";
00183     const std::string bottomTempProjectedNCField = "a2oTbot_projnocons";
00184 
00185     // Attributes for tag data storage
00186     int tagIndex[4];
00187     // int entTypes[2] = {1, 1}; /* both on elements; */
00188     int tagTypes[2]  = { DENSE_DOUBLE, DENSE_DOUBLE };
00189     int atmCompNDoFs = disc_orders[0] * disc_orders[0], ocnCompNDoFs = 1 /*FV*/;
00190 
00191     ierr = iMOAB_DefineTagStorage( atmPID, bottomTempField.c_str(), &tagTypes[0], &atmCompNDoFs, &tagIndex[0] );
00192     CHECKIERR( ierr, "failed to define the field tag on ATM" );
00193 
00194     ierr = iMOAB_DefineTagStorage( atmPID, bottomTempFieldATM.c_str(), &tagTypes[0], &atmCompNDoFs, &tagIndex[0] );
00195     CHECKIERR( ierr, "failed to define the field tag on ATM" );
00196 
00197     ierr =
00198         iMOAB_DefineTagStorage( ocnPID, bottomTempProjectedField.c_str(), &tagTypes[1], &ocnCompNDoFs, &tagIndex[1] );
00199     CHECKIERR( ierr, "failed to define the field tag on OCN" );
00200 
00201     ierr =
00202         iMOAB_DefineTagStorage( ocnPID, bottomTempProjectedNCField.c_str(), &tagTypes[1], &ocnCompNDoFs, &tagIndex[2] );
00203     CHECKIERR( ierr, "failed to define the field tag on OCN" );
00204 
00205 #ifdef ENABLE_ATMLND_COUPLING
00206     ierr =
00207         iMOAB_DefineTagStorage( lndPID, bottomTempProjectedField.c_str(), &tagTypes[1], &ocnCompNDoFs, &tagIndex[3] );
00208     CHECKIERR( ierr, "failed to define the field tag on LND" );
00209 #endif
00210 
00211     /* Next compute the mesh intersection on the sphere between the source (ATM) and target (OCN)
00212      * meshes */
00213     ierr = iMOAB_ComputeMeshIntersectionOnSphere( atmPID, ocnPID, atmocnPID );
00214     CHECKIERR( ierr, "failed to compute mesh intersection between ATM and OCN" );
00215 
00216 #ifdef ENABLE_ATMLND_COUPLING
00217     /* Next compute the mesh intersection on the sphere between the source (ATM) and target (LND)
00218      * meshes */
00219     ierr = iMOAB_ComputePointDoFIntersection( atmPID, lndPID, atmlndPID );
00220     CHECKIERR( ierr, "failed to compute point-cloud mapping ATM-LND" );
00221 
00222     /* Next compute the mesh intersection on the sphere between the source (LND) and target (ATM)
00223      * meshes */
00224     ierr = iMOAB_ComputePointDoFIntersection( lndPID, atmPID, lndatmPID );
00225     CHECKIERR( ierr, "failed to compute point-cloud mapping LND-ATM" );
00226 #endif
00227 
00228     /* We have the mesh intersection now. Let us compute the remapping weights */
00229     fNoConserve = 0;
00230 
00231     /* Compute the weights to preoject the solution from ATM component to OCN compoenent */
00232     ierr = iMOAB_ComputeScalarProjectionWeights( atmocnPID, weights_identifiers[0].c_str(), disc_methods[0].c_str(),
00233                                                  &disc_orders[0], disc_methods[1].c_str(), &disc_orders[1], &fNoBubble,
00234                                                  &fMonotoneTypeID, &fVolumetric, &fInverseDistanceMap, &fNoConserve,
00235                                                  &fValidate, dof_tag_names[0].c_str(), dof_tag_names[1].c_str() );
00236     CHECKIERR( ierr, "failed to compute remapping projection weights for ATM-OCN scalar "
00237                      "non-conservative field" );
00238 
00239     // Let us now write the map file to disk and then read it back to test the I/O API in iMOAB
00240 #ifdef MOAB_HAVE_NETCDF
00241     {
00242         const std::string atmocn_map_file_name = "atm_ocn_map.nc";
00243         ierr =
00244             iMOAB_WriteMappingWeightsToFile( atmocnPID, weights_identifiers[0].c_str(), atmocn_map_file_name.c_str() );
00245         CHECKIERR( ierr, "failed to write map file to disk" );
00246 
00247         const std::string intx_from_file_identifier = "map-from-file";
00248         // use old reader, coupler PID is -1
00249         int dummyCpl     = -1;
00250         int dummy_rowcol = -1;
00251         int dummyType    = 0;
00252         ierr             = iMOAB_LoadMappingWeightsFromFile( atmocnPID, &dummyCpl, &dummy_rowcol, &dummyType,
00253                                                              intx_from_file_identifier.c_str(), atmocn_map_file_name.c_str() );
00254         CHECKIERR( ierr, "failed to load map file from disk" );
00255     }
00256 #endif
00257 #ifdef ENABLE_ATMLND_COUPLING
00258     fValidate = 0;
00259     /* Compute the weights to preoject the solution from ATM component to LND compoenent */
00260     ierr = iMOAB_ComputeScalarProjectionWeights( atmlndPID, weights_identifiers[1].c_str(), disc_methods[0].c_str(),
00261                                                  &disc_orders[0], disc_methods[2].c_str(), &disc_orders[2], &fNoBubble,
00262                                                  &fMonotoneTypeID, &fVolumetric, &fInverseDistanceMap, &fNoConserve,
00263                                                  &fValidate, dof_tag_names[0].c_str(), dof_tag_names[2].c_str() );
00264     CHECKIERR( ierr, "failed to compute remapping projection weights for ATM-LND scalar "
00265                      "non-conservative field" );
00266 
00267     /* Compute the weights to preoject the solution from ATM component to LND compoenent */
00268     ierr = iMOAB_ComputeScalarProjectionWeights( lndatmPID, weights_identifiers[1].c_str(), disc_methods[2].c_str(),
00269                                                  &disc_orders[2], disc_methods[0].c_str(), &disc_orders[0], &fNoBubble,
00270                                                  &fMonotoneTypeID, &fVolumetric, &fInverseDistanceMap, &fNoConserve,
00271                                                  &fValidate, dof_tag_names[2].c_str(), dof_tag_names[0].c_str() );
00272     CHECKIERR( ierr, "failed to compute remapping projection weights for LND-ATM scalar "
00273                      "non-conservative field" );
00274 #endif
00275     /* We have the mesh intersection now. Let us compute the remapping weights */
00276     fNoConserve = 0;
00277     ierr = iMOAB_ComputeScalarProjectionWeights( atmocnPID, weights_identifiers[2].c_str(), disc_methods[0].c_str(),
00278                                                  &disc_orders[0], disc_methods[1].c_str(), &disc_orders[1], &fNoBubble,
00279                                                  &fMonotoneTypeID, &fVolumetric, &fInverseDistanceMap, &fNoConserve,
00280                                                  &fValidate, dof_tag_names[0].c_str(), dof_tag_names[1].c_str() );
00281     CHECKIERR( ierr, "failed to compute remapping projection weights for scalar conservative field" );
00282 
00283     /* We have the remapping weights now. Let us apply the weights onto the tag we defined
00284        on the srouce mesh and get the projection on the target mesh */
00285     ierr = iMOAB_ApplyScalarProjectionWeights( atmocnPID, weights_identifiers[0].c_str(), bottomTempField.c_str(),
00286                                                bottomTempProjectedNCField.c_str() );
00287     CHECKIERR( ierr, "failed to apply projection weights for scalar non-conservative field" );
00288 
00289     /* We have the remapping weights now. Let us apply the weights onto the tag we defined
00290        on the srouce mesh and get the projection on the target mesh */
00291     ierr = iMOAB_ApplyScalarProjectionWeights( atmocnPID, weights_identifiers[2].c_str(), bottomTempField.c_str(),
00292                                                bottomTempProjectedField.c_str() );
00293     CHECKIERR( ierr, "failed to apply projection weights for scalar conservative field" );
00294 
00295     if( gen_baseline )
00296     {
00297         // get temp field on ocean, from conservative, the global ids, and dump to the baseline file
00298         // first get GlobalIds from ocn, and fields:
00299         ierr = iMOAB_GetMeshInfo( ocnPID, nverts, nelem, 0, 0, 0 );
00300         CHECKIERR( ierr, "failed to get ocn mesh info" );
00301         std::vector< int > gidElems;
00302         gidElems.resize( nelem[2] );
00303         std::vector< double > tempElems;
00304         tempElems.resize( nelem[2] );
00305         // get global id storage
00306         const std::string GidStr = "GLOBAL_ID";  // hard coded too
00307 
00308         int tag_type = DENSE_INTEGER, ncomp = 1, tagInd = 0;
00309         // we should not have to define global id tag, it is always defined
00310         ierr = iMOAB_DefineTagStorage( ocnPID, GidStr.c_str(), &tag_type, &ncomp, &tagInd );
00311         CHECKIERR( ierr, "failed to define global id tag" );
00312         int ent_type = 1;
00313         ierr         = iMOAB_GetIntTagStorage( ocnPID, GidStr.c_str(), &nelem[2], &ent_type, &gidElems[0] );
00314         CHECKIERR( ierr, "failed to get global ids" );
00315         ierr =
00316             iMOAB_GetDoubleTagStorage( ocnPID, bottomTempProjectedField.c_str(), &nelem[2], &ent_type, &tempElems[0] );
00317         CHECKIERR( ierr, "failed to get temperature field" );
00318         std::fstream fs;
00319         fs.open( baseline.c_str(), std::fstream::out );
00320         fs << std::setprecision( 15 );  // maximum precision for doubles
00321         for( int i = 0; i < nelem[2]; i++ )
00322             fs << gidElems[i] << " " << tempElems[i] << "\n";
00323         fs.close();
00324     }
00325     if( !no_test_against_baseline )
00326     {
00327         // get temp field on ocean, from conservative, the global ids, and dump to the baseline file
00328         // first get GlobalIds from ocn, and fields:
00329         ierr = iMOAB_GetMeshInfo( ocnPID, nverts, nelem, 0, 0, 0 );
00330         CHECKIERR( ierr, "failed to get ocn mesh info" );
00331         std::vector< int > gidElems;
00332         gidElems.resize( nelem[2] );
00333         std::vector< double > tempElems;
00334         tempElems.resize( nelem[2] );
00335         // get global id storage
00336         const std::string GidStr = "GLOBAL_ID";  // hard coded too
00337 
00338         int tag_type = DENSE_INTEGER, ncomp = 1, tagInd = 0;
00339         ierr = iMOAB_DefineTagStorage( ocnPID, GidStr.c_str(), &tag_type, &ncomp, &tagInd );
00340         CHECKIERR( ierr, "failed to define global id tag" );
00341 
00342         int ent_type = 1;
00343         ierr         = iMOAB_GetIntTagStorage( ocnPID, GidStr.c_str(), &nelem[2], &ent_type, &gidElems[0] );
00344         CHECKIERR( ierr, "failed to get global ids" );
00345         ierr =
00346             iMOAB_GetDoubleTagStorage( ocnPID, bottomTempProjectedField.c_str(), &nelem[2], &ent_type, &tempElems[0] );
00347         CHECKIERR( ierr, "failed to get temperature field" );
00348         int err_code = 1;
00349         check_baseline_file( baseline, gidElems, tempElems, 1.e-9, err_code );
00350         if( 0 == err_code ) std::cout << " passed baseline test 1\n";
00351     }
00352 #ifdef ENABLE_ATMLND_COUPLING
00353     /* We have the remapping weights now. Let us apply the weights onto the tag we defined
00354        on the srouce mesh and get the projection on the target mesh */
00355     ierr = iMOAB_ApplyScalarProjectionWeights( atmlndPID, weights_identifiers[1].c_str(), bottomTempField.c_str(),
00356                                                bottomTempProjectedField.c_str() );
00357     CHECKIERR( ierr, "failed to apply projection weights for ATM-LND scalar field" );
00358 
00359     /* We have the remapping weights now. Let us apply the weights onto the tag we defined
00360        on the srouce mesh and get the projection on the target mesh */
00361     ierr = iMOAB_ApplyScalarProjectionWeights( lndatmPID, weights_identifiers[1].c_str(),
00362                                                bottomTempProjectedField.c_str(), bottomTempFieldATM.c_str() );
00363     CHECKIERR( ierr, "failed to apply projection weights for LND-ATM scalar field" );
00364 #endif
00365     /*
00366      * the file can be written in parallel, and it will contain additional tags defined by the user
00367      * we may extend the method to write only desired tags to the file
00368      */
00369     {
00370         // free allocated data
00371 #ifdef ENABLE_ATMLND_COUPLING
00372         char outputFileAtmTgt[] = "fIntxAtmTarget.h5m";
00373 #endif
00374         char outputFileOcnTgt[] = "fIntxOcnTarget.h5m";
00375 #ifdef ENABLE_ATMLND_COUPLING
00376         char outputFileLndTgt[] = "fIntxLndTarget.h5m";
00377 #endif
00378         char writeOptions[] = "";
00379 
00380         // Write out all the component meshes to disk for verification
00381         ierr = iMOAB_WriteMesh( ocnPID, outputFileOcnTgt, writeOptions );
00382         CHECKIERR( ierr, "failed to write OCN mesh and data" );
00383 #ifdef ENABLE_ATMLND_COUPLING
00384         ierr = iMOAB_WriteMesh( lndPID, outputFileLndTgt, writeOptions );
00385         CHECKIERR( ierr, "failed to write LND mesh and data" );
00386         ierr = iMOAB_WriteMesh( atmPID, outputFileAtmTgt, writeOptions );
00387         CHECKIERR( ierr, "failed to write ATM mesh and data" );
00388 #endif
00389     }
00390 
00391     /*
00392      * deregistering application will delete all mesh entities associated with the application and
00393      * will free allocated tag storage.
00394      */
00395 #ifdef ENABLE_ATMLND_COUPLING
00396     ierr = iMOAB_DeregisterApplication( lndPID );
00397     CHECKIERR( ierr, "failed to de-register application3" );
00398 #endif
00399     ierr = iMOAB_DeregisterApplication( ocnPID );
00400     CHECKIERR( ierr, "failed to de-register application2" );
00401     ierr = iMOAB_DeregisterApplication( atmPID );
00402     CHECKIERR( ierr, "failed to de-register application1" );
00403 
00404     /*
00405      * this method will delete MOAB instance
00406      */
00407     ierr = iMOAB_Finalize();
00408     CHECKIERR( ierr, "failed to finalize MOAB" );
00409 
00410 #ifdef MOAB_HAVE_MPI
00411     MPI_Finalize();
00412 #endif
00413 
00414     return 0;
00415 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines