MOAB: Mesh Oriented datABase  (version 5.2.1)
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 
00014 #define CHECKIERR( ierr, message ) \
00015     if( 0 != ierr )                \
00016     {                              \
00017         printf( "%s", message );   \
00018         return 1;                  \
00019     }
00020 
00021 // this test will be run in serial only
00022 int main( int argc, char* argv[] )
00023 {
00024     ErrCode ierr;
00025     std::string atmFilename, ocnFilename, lndFilename, mapFilename;
00026 
00027     atmFilename = TestDir + "/wholeATM_T.h5m";
00028     ocnFilename = TestDir + "/recMeshOcn.h5m";
00029     lndFilename = TestDir + "/wholeLnd.h5m";
00030     // mapFilename = TestDir + "/outCS5ICOD5_map.nc";
00031 
00032     ProgOptions opts;
00033     opts.addOpt< std::string >( "atmosphere,t", "atm mesh filename (source)", &atmFilename );
00034     opts.addOpt< std::string >( "ocean,m", "ocean mesh filename (target)", &ocnFilename );
00035     opts.addOpt< std::string >( "land,l", "land mesh filename (target)", &lndFilename );
00036 
00037     opts.parseCommandLine( argc, argv );
00038 
00039     {
00040         std::cout << " atm file: " << atmFilename << "\n ocn file: " << ocnFilename << "\n lnd file: " << lndFilename
00041                   << "\n";
00042     }
00043 
00044 #ifdef MOAB_HAVE_MPI
00045     MPI_Init( &argc, &argv );
00046     MPI_Comm comm = MPI_COMM_WORLD;
00047 #endif
00048     /*
00049      * MOAB needs to be initialized; A MOAB instance will be created, and will be used by each
00050      * application in this framework. There is no parallel context yet.
00051      */
00052     ierr = iMOAB_Initialize( argc, argv );
00053     CHECKIERR( ierr, "failed to initialize MOAB" );
00054 
00055     int atmAppID, ocnAppID, lndAppID, atmocnAppID, atmlndAppID, lndatmAppID;
00056     iMOAB_AppID atmPID    = &atmAppID;
00057     iMOAB_AppID ocnPID    = &ocnAppID;
00058     iMOAB_AppID lndPID    = &lndAppID;
00059     iMOAB_AppID atmocnPID = &atmocnAppID;
00060     iMOAB_AppID atmlndPID = &atmlndAppID;
00061     iMOAB_AppID lndatmPID = &lndatmAppID;
00062 
00063     /*
00064      * Each application has to be registered once. A mesh set and a parallel communicator will be
00065      * associated with each application. A unique application id will be returned, and will be used
00066      * for all future mesh operations/queries.
00067      */
00068     int atmCompID = 10, ocnCompID = 20, lndCompID = 30, atmocnCompID = 100, atmlndCompID = 102, lndatmCompID = 103;
00069     ierr = iMOAB_RegisterApplication( "ATM-APP",
00070 #ifdef MOAB_HAVE_MPI
00071                                       &comm,
00072 #endif
00073                                       &atmCompID, atmPID );
00074     CHECKIERR( ierr, "failed to register application1" );
00075 
00076     ierr = iMOAB_RegisterApplication( "OCN-APP",
00077 #ifdef MOAB_HAVE_MPI
00078                                       &comm,
00079 #endif
00080                                       &ocnCompID, ocnPID );
00081     CHECKIERR( ierr, "failed to register application2" );
00082 
00083     ierr = iMOAB_RegisterApplication( "LND-APP",
00084 #ifdef MOAB_HAVE_MPI
00085                                       &comm,
00086 #endif
00087                                       &lndCompID, lndPID );
00088     CHECKIERR( ierr, "failed to register application3" );
00089 
00090     ierr = iMOAB_RegisterApplication( "ATM-OCN-CPL",
00091 #ifdef MOAB_HAVE_MPI
00092                                       &comm,
00093 #endif
00094                                       &atmocnCompID, atmocnPID );
00095     CHECKIERR( ierr, "failed to register application4" );
00096 
00097     ierr = iMOAB_RegisterApplication( "ATM-LND-CPL",
00098 #ifdef MOAB_HAVE_MPI
00099                                       &comm,
00100 #endif
00101                                       &atmlndCompID, atmlndPID );
00102     CHECKIERR( ierr, "failed to register application5" );
00103 
00104     ierr = iMOAB_RegisterApplication( "LND-ATM-CPL",
00105 #ifdef MOAB_HAVE_MPI
00106                                       &comm,
00107 #endif
00108                                       &lndatmCompID, lndatmPID );
00109     CHECKIERR( ierr, "failed to register application5" );
00110 
00111     const char* read_opts = "";
00112     int num_ghost_layers  = 0;
00113     /*
00114      * Loading the mesh is a parallel IO operation. Ghost layers can be exchanged too, and default
00115      * MOAB sets are augmented with ghost elements. By convention, blocks correspond to MATERIAL_SET
00116      * sets, side sets with NEUMANN_SET sets, node sets with DIRICHLET_SET sets. Each element and
00117      * vertex entity should have a GLOBAL ID tag in the file, which will be available for visible
00118      * entities
00119      */
00120     ierr = iMOAB_LoadMesh( atmPID, atmFilename.c_str(), read_opts, &num_ghost_layers, atmFilename.size(),
00121                            strlen( read_opts ) );
00122     CHECKIERR( ierr, "failed to load mesh" );
00123     ierr = iMOAB_LoadMesh( ocnPID, ocnFilename.c_str(), read_opts, &num_ghost_layers, ocnFilename.size(),
00124                            strlen( read_opts ) );
00125     CHECKIERR( ierr, "failed to load mesh" );
00126     ierr = iMOAB_LoadMesh( lndPID, lndFilename.c_str(), read_opts, &num_ghost_layers, lndFilename.size(),
00127                            strlen( read_opts ) );
00128     CHECKIERR( ierr, "failed to load mesh" );
00129 
00130     int nverts[3], nelem[3];
00131     /*
00132      * Each process in the communicator will have access to a local mesh instance, which will
00133      * contain the original cells in the local partition and ghost entities. Number of vertices,
00134      * primary cells, visible blocks, number of sidesets and nodesets boundary conditions will be
00135      * returned in size 3 arrays, for local, ghost and total numbers.
00136      */
00137     ierr = iMOAB_GetMeshInfo( atmPID, nverts, nelem, 0, 0, 0 );
00138     CHECKIERR( ierr, "failed to get mesh info" );
00139     printf( "Atmosphere Component Mesh: %d vertices and %d elements\n", nverts[0], nelem[0] );
00140 
00141     ierr = iMOAB_GetMeshInfo( ocnPID, nverts, nelem, 0, 0, 0 );
00142     CHECKIERR( ierr, "failed to get mesh info" );
00143     printf( "Ocean Component Mesh: %d vertices and %d elements\n", nverts[0], nelem[0] );
00144 
00145     ierr = iMOAB_GetMeshInfo( lndPID, nverts, nelem, 0, 0, 0 );
00146     CHECKIERR( ierr, "failed to get mesh info" );
00147     printf( "Land Component Mesh: %d vertices and %d elements\n", nverts[0], nelem[0] );
00148 
00149     /*
00150      * The 2 tags used in this example exist in the file, already.
00151      * If a user needs a new tag, it can be defined with the same call as this one
00152      * this method, iMOAB_DefineTagStorage, will return a local index for the tag.
00153      * The name of the tag is case sensitive.
00154      * This method is collective.
00155      */
00156     // int disc_orders[2] = {1, 1};
00157     // const char* disc_methods[2] = {"fv", "fv"};
00158     // const char* dof_tag_names[2] = {"GLOBAL_ID", "GLOBAL_ID"};
00159     int disc_orders[3]           = { 4, 1, 1 };
00160     int fMonotoneTypeID = 0, fVolumetric = 0, fValidate = 1, fNoConserve = 0;
00161 
00162     const std::string disc_methods[3]        = { "cgll", "fv", "pcloud" };
00163     const std::string dof_tag_names[3]       = { "GLOBAL_DOFS", "GLOBAL_ID", "GLOBAL_ID" };
00164     const std::string weights_identifiers[3] = { "scalar", "scalar_pointcloud", "scalar_conservative" };
00165 
00166     const std::string bottomTempField            = "a2oTbot";
00167     const std::string bottomTempFieldATM         = "a2oTbotATM";
00168     const std::string bottomTempProjectedField   = "a2oTbot_proj";
00169     const std::string bottomTempProjectedNCField = "a2oTbot_projnocons";
00170 
00171     // Attributes for tag data storage
00172     int tagIndex[4];
00173     // int entTypes[2] = {1, 1}; /* both on elements; */
00174     int tagTypes[2]  = { DENSE_DOUBLE, DENSE_DOUBLE };
00175     int atmCompNDoFs = disc_orders[0] * disc_orders[0], ocnCompNDoFs = 1 /*FV*/;
00176 
00177     ierr = iMOAB_DefineTagStorage( atmPID, bottomTempField.c_str(), &tagTypes[0], &atmCompNDoFs, &tagIndex[0],
00178                                    bottomTempField.size() );
00179     CHECKIERR( ierr, "failed to define the field tag on ATM" );
00180 
00181     ierr = iMOAB_DefineTagStorage( atmPID, bottomTempFieldATM.c_str(), &tagTypes[0], &atmCompNDoFs, &tagIndex[0],
00182                                    bottomTempFieldATM.size() );
00183     CHECKIERR( ierr, "failed to define the field tag on ATM" );
00184 
00185     ierr = iMOAB_DefineTagStorage( ocnPID, bottomTempProjectedField.c_str(), &tagTypes[1], &ocnCompNDoFs, &tagIndex[1],
00186                                    bottomTempProjectedField.size() );
00187     CHECKIERR( ierr, "failed to define the field tag on OCN" );
00188 
00189     ierr = iMOAB_DefineTagStorage( ocnPID, bottomTempProjectedNCField.c_str(), &tagTypes[1], &ocnCompNDoFs, &tagIndex[2],
00190                                    bottomTempProjectedNCField.size() );
00191     CHECKIERR( ierr, "failed to define the field tag on OCN" );
00192 
00193     ierr = iMOAB_DefineTagStorage( lndPID, bottomTempProjectedField.c_str(), &tagTypes[1], &ocnCompNDoFs, &tagIndex[3],
00194                                    bottomTempProjectedField.size() );
00195     CHECKIERR( ierr, "failed to define the field tag on LND" );
00196 
00197     /* Next compute the mesh intersection on the sphere between the source (ATM) and target (OCN)
00198      * meshes */
00199     ierr = iMOAB_ComputeMeshIntersectionOnSphere( atmPID, ocnPID, atmocnPID );
00200     CHECKIERR( ierr, "failed to compute mesh intersection between ATM and OCN" );
00201 
00202     /* Next compute the mesh intersection on the sphere between the source (ATM) and target (LND)
00203      * meshes */
00204     ierr = iMOAB_ComputePointDoFIntersection( atmPID, lndPID, atmlndPID );
00205     CHECKIERR( ierr, "failed to compute point-cloud mapping ATM-LND" );
00206 
00207     /* Next compute the mesh intersection on the sphere between the source (LND) and target (ATM)
00208      * meshes */
00209     ierr = iMOAB_ComputePointDoFIntersection( lndPID, atmPID, lndatmPID );
00210     CHECKIERR( ierr, "failed to compute point-cloud mapping LND-ATM" );
00211 
00212     /* We have the mesh intersection now. Let us compute the remapping weights */
00213     fNoConserve = 0;
00214     /* Compute the weights to preoject the solution from ATM component to OCN compoenent */
00215     ierr = iMOAB_ComputeScalarProjectionWeights( atmocnPID, weights_identifiers[0].c_str(), disc_methods[0].c_str(), &disc_orders[0],
00216                                                  disc_methods[1].c_str(), &disc_orders[1], &fMonotoneTypeID, &fVolumetric,
00217                                                  &fNoConserve, &fValidate, dof_tag_names[0].c_str(), dof_tag_names[1].c_str(),
00218                                                  weights_identifiers[0].size(), disc_methods[0].size(),
00219                                                  disc_methods[1].size(), dof_tag_names[0].size(),
00220                                                  dof_tag_names[1].size() );
00221     CHECKIERR( ierr, "failed to compute remapping projection weights for ATM-OCN scalar "
00222                      "non-conservative field" );
00223 
00224     // Let us now write the map file to disk and then read it back to test the I/O API in iMOAB
00225 #ifdef MOAB_HAVE_NETCDF
00226     {
00227         const std::string atmocn_map_file_name = "atm_ocn_map.nc";
00228         ierr = iMOAB_WriteMappingWeightsToFile( atmocnPID, weights_identifiers[0].c_str(), atmocn_map_file_name.c_str(),
00229                                                 weights_identifiers[0].size(), atmocn_map_file_name.size() );
00230         CHECKIERR( ierr, "failed to write map file to disk" );
00231 
00232         const std::string intx_from_file_identifier = "map-from-file";
00233         ierr = iMOAB_LoadMappingWeightsFromFile( atmocnPID, intx_from_file_identifier.c_str(),
00234                                                  atmocn_map_file_name.c_str(), NULL, NULL, NULL,
00235                                                  intx_from_file_identifier.size(), atmocn_map_file_name.size() );
00236         CHECKIERR( ierr, "failed to load map file from disk" );
00237     }
00238 #endif
00239 
00240     /* Compute the weights to preoject the solution from ATM component to LND compoenent */
00241     ierr = iMOAB_ComputeScalarProjectionWeights( atmlndPID, weights_identifiers[1].c_str(), disc_methods[0].c_str(), &disc_orders[0],
00242                                                  disc_methods[2].c_str(), &disc_orders[2], &fMonotoneTypeID, &fVolumetric,
00243                                                  &fNoConserve, &fValidate, dof_tag_names[0].c_str(), dof_tag_names[2].c_str(),
00244                                                  weights_identifiers[1].size(), disc_methods[0].size(),
00245                                                  disc_methods[2].size(), dof_tag_names[0].size(),
00246                                                  dof_tag_names[2].size() );
00247     CHECKIERR( ierr, "failed to compute remapping projection weights for ATM-LND scalar "
00248                      "non-conservative field" );
00249 
00250     /* Compute the weights to preoject the solution from ATM component to LND compoenent */
00251     ierr = iMOAB_ComputeScalarProjectionWeights( lndatmPID, weights_identifiers[1].c_str(), disc_methods[2].c_str(), &disc_orders[2],
00252                                                  disc_methods[0].c_str(), &disc_orders[0], &fMonotoneTypeID, &fVolumetric,
00253                                                  &fNoConserve, &fValidate, dof_tag_names[2].c_str(), dof_tag_names[0].c_str(),
00254                                                  weights_identifiers[1].size(), disc_methods[2].size(),
00255                                                  disc_methods[0].size(), dof_tag_names[2].size(),
00256                                                  dof_tag_names[0].size() );
00257     CHECKIERR( ierr, "failed to compute remapping projection weights for LND-ATM scalar "
00258                      "non-conservative field" );
00259 
00260     /* We have the mesh intersection now. Let us compute the remapping weights */
00261     fNoConserve = 0;
00262     ierr = iMOAB_ComputeScalarProjectionWeights( atmocnPID, weights_identifiers[2].c_str(), disc_methods[0].c_str(), &disc_orders[0],
00263                                                  disc_methods[1].c_str(), &disc_orders[1], &fMonotoneTypeID, &fVolumetric,
00264                                                  &fNoConserve, &fValidate, dof_tag_names[0].c_str(), dof_tag_names[1].c_str(),
00265                                                  weights_identifiers[1].size(), disc_methods[0].size(),
00266                                                  disc_methods[1].size(), dof_tag_names[0].size(),
00267                                                  dof_tag_names[1].size() );
00268     CHECKIERR( ierr, "failed to compute remapping projection weights for scalar conservative field" );
00269 
00270     /* We have the remapping weights now. Let us apply the weights onto the tag we defined
00271        on the srouce mesh and get the projection on the target mesh */
00272     ierr = iMOAB_ApplyScalarProjectionWeights( atmocnPID, weights_identifiers[0].c_str(), bottomTempField.c_str(),
00273                                                bottomTempProjectedNCField.c_str(), weights_identifiers[0].size(),
00274                                                bottomTempField.size(), bottomTempProjectedNCField.size() );
00275     CHECKIERR( ierr, "failed to apply projection weights for scalar non-conservative field" );
00276 
00277     /* We have the remapping weights now. Let us apply the weights onto the tag we defined
00278        on the srouce mesh and get the projection on the target mesh */
00279     ierr = iMOAB_ApplyScalarProjectionWeights( atmocnPID, weights_identifiers[2].c_str(), bottomTempField.c_str(),
00280                                                bottomTempProjectedField.c_str(), weights_identifiers[1].size(),
00281                                                bottomTempField.size(), bottomTempProjectedField.size() );
00282     CHECKIERR( ierr, "failed to apply projection weights for scalar conservative field" );
00283 
00284     /* We have the remapping weights now. Let us apply the weights onto the tag we defined
00285        on the srouce mesh and get the projection on the target mesh */
00286     ierr = iMOAB_ApplyScalarProjectionWeights( atmlndPID, weights_identifiers[1].c_str(), bottomTempField.c_str(),
00287                                                bottomTempProjectedField.c_str(), weights_identifiers[1].size(),
00288                                                bottomTempField.size(), bottomTempProjectedField.size() );
00289     CHECKIERR( ierr, "failed to apply projection weights for ATM-LND scalar field" );
00290 
00291     /* We have the remapping weights now. Let us apply the weights onto the tag we defined
00292        on the srouce mesh and get the projection on the target mesh */
00293     ierr = iMOAB_ApplyScalarProjectionWeights( lndatmPID, weights_identifiers[1].c_str(), bottomTempProjectedField.c_str(),
00294                                                bottomTempFieldATM.c_str(), weights_identifiers[1].size(),
00295                                                bottomTempField.size(), bottomTempProjectedField.size() );
00296     CHECKIERR( ierr, "failed to apply projection weights for LND-ATM scalar field" );
00297 
00298     /*
00299      * the file can be written in parallel, and it will contain additional tags defined by the user
00300      * we may extend the method to write only desired tags to the file
00301      */
00302     {
00303         // free allocated data
00304         char outputFileAtmTgt[] = "fIntxAtmTarget.h5m";
00305         char outputFileOcnTgt[] = "fIntxOcnTarget.h5m";
00306         char outputFileLndTgt[] = "fIntxLndTarget.h5m";
00307         char writeOptions[]     = "";
00308 
00309         ierr = iMOAB_WriteMesh( ocnPID, outputFileOcnTgt, writeOptions, strlen( outputFileOcnTgt ),
00310                                 strlen( writeOptions ) );
00311 
00312         ierr = iMOAB_WriteMesh( lndPID, outputFileLndTgt, writeOptions, strlen( outputFileLndTgt ),
00313                                 strlen( writeOptions ) );
00314 
00315         ierr = iMOAB_WriteMesh( atmPID, outputFileAtmTgt, writeOptions, strlen( outputFileAtmTgt ),
00316                                 strlen( writeOptions ) );
00317     }
00318 
00319     /*
00320      * deregistering application will delete all mesh entities associated with the application and
00321      * will free allocated tag storage.
00322      */
00323     ierr = iMOAB_DeregisterApplication( lndPID );
00324     CHECKIERR( ierr, "failed to de-register application3" );
00325     ierr = iMOAB_DeregisterApplication( ocnPID );
00326     CHECKIERR( ierr, "failed to de-register application2" );
00327     ierr = iMOAB_DeregisterApplication( atmPID );
00328     CHECKIERR( ierr, "failed to de-register application1" );
00329 
00330     /*
00331      * this method will delete MOAB instance
00332      */
00333     ierr = iMOAB_Finalize();
00334     CHECKIERR( ierr, "failed to finalize MOAB" );
00335 
00336 #ifdef MOAB_HAVE_MPI
00337     MPI_Finalize();
00338 #endif
00339 
00340     return 0;
00341 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines