|
MOAB: Mesh Oriented datABase
(version 5.4.0)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.0)
|
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 }
1.7.6.1 
