MOAB: Mesh Oriented datABase  (version 5.5.0)
read_mpas_nc.cpp
Go to the documentation of this file.
1 #include "TestUtil.hpp"
2 #include "moab/Core.hpp"
3 #include "moab/ReadUtilIface.hpp"
4 #include "MBTagConventions.hpp"
5 
6 using namespace moab;
7 
8 std::string example = TestDir + "unittest/io/mpasx1.642.t.2.nc";
9 
10 #ifdef MOAB_HAVE_MPI
11 #include "moab_mpi.h"
12 #include "moab/ParallelComm.hpp"
13 #endif
14 
15 void test_read_all();
16 void test_read_onevar();
17 void test_read_onetimestep();
18 void test_read_nomesh();
19 void test_read_novars();
20 void test_read_no_mixed_elements(); // Test read option NO_MIXED_ELEMENTS
21 void test_read_no_edges(); // Test read option NO_EDGES
22 void test_gather_onevar(); // Test gather set with one variable
23 
24 void get_options( std::string& opts );
25 
26 const double eps = 1e-20;
27 
28 int main( int argc, char* argv[] )
29 {
30  int result = 0;
31 
32 #ifdef MOAB_HAVE_MPI
33  int fail = MPI_Init( &argc, &argv );
34  if( fail ) return 1;
35 #else
36  argv[0] = argv[argc - argc]; // To remove the warnings in serial mode about unused variables
37 #endif
38 
39  result += RUN_TEST( test_read_all );
40  result += RUN_TEST( test_read_onevar );
41  result += RUN_TEST( test_read_onetimestep );
42  result += RUN_TEST( test_read_nomesh );
43  result += RUN_TEST( test_read_novars );
44  result += RUN_TEST( test_read_no_mixed_elements );
45  result += RUN_TEST( test_read_no_edges );
46  result += RUN_TEST( test_gather_onevar );
47 
48 #ifdef MOAB_HAVE_MPI
49  fail = MPI_Finalize();
50  if( fail ) return 1;
51 #endif
52 
53  return result;
54 }
55 
56 void test_read_all()
57 {
58  Core moab;
59  Interface& mb = moab;
60 
61  std::string opts;
62  get_options( opts );
63 
64  // Read mesh and read all variables at all timesteps
65  ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
66 
67 #ifdef MOAB_HAVE_MPI
68  ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
69  int procs = pcomm->proc_config().proc_size();
70 #else
71  int procs = 1;
72 #endif
73 
74  // Make check runs this test on one processor
75  if( 1 == procs )
76  {
77  // For each tag, check two values
78  double val[2];
79 
80  // Check tags for vertex variable vorticity
81  Tag vorticity_tag0, vorticity_tag1;
82  rval = mb.tag_get_handle( "vorticity0", 1, MB_TYPE_DOUBLE, vorticity_tag0 );CHECK_ERR( rval );
83  rval = mb.tag_get_handle( "vorticity1", 1, MB_TYPE_DOUBLE, vorticity_tag1 );CHECK_ERR( rval );
84 
85  // Get vertices (1280 edges)
86  Range verts;
87  rval = mb.get_entities_by_type( 0, MBVERTEX, verts );CHECK_ERR( rval );
88  CHECK_EQUAL( (size_t)1280, verts.size() );
89  CHECK_EQUAL( (size_t)1, verts.psize() );
90 
91  // Check vorticity tag values on first two vertices
92  EntityHandle vert_ents[] = { verts[0], verts[1] };
93  rval = mb.tag_get_data( vorticity_tag0, vert_ents, 2, val );CHECK_ERR( rval );
94  CHECK_REAL_EQUAL( 1.1, val[0], eps );
95  CHECK_REAL_EQUAL( 1.2, val[1], eps );
96  rval = mb.tag_get_data( vorticity_tag1, vert_ents, 2, val );CHECK_ERR( rval );
97  CHECK_REAL_EQUAL( 2.1, val[0], eps );
98  CHECK_REAL_EQUAL( 2.2, val[1], eps );
99 
100  // Check tags for edge variable u
101  Tag u_tag0, u_tag1;
102  rval = mb.tag_get_handle( "u0", 1, MB_TYPE_DOUBLE, u_tag0 );CHECK_ERR( rval );
103  rval = mb.tag_get_handle( "u1", 1, MB_TYPE_DOUBLE, u_tag1 );CHECK_ERR( rval );
104 
105  // Get edges (1920 edges)
106  Range edges;
107  rval = mb.get_entities_by_type( 0, MBEDGE, edges );CHECK_ERR( rval );
108  CHECK_EQUAL( (size_t)1920, edges.size() );
109  CHECK_EQUAL( (size_t)1, edges.psize() );
110 
111  // Check u tag values on two specified edges
112  EntityHandle edge_ents[] = { edges[5], edges[6] };
113  rval = mb.tag_get_data( u_tag0, edge_ents, 2, val );CHECK_ERR( rval );
114  CHECK_REAL_EQUAL( 1.113138721544778, val[0], eps );
115  CHECK_REAL_EQUAL( -1.113138721930009, val[1], eps );
116  rval = mb.tag_get_data( u_tag1, edge_ents, 2, val );CHECK_ERR( rval );
117  CHECK_REAL_EQUAL( 2.113138721544778, val[0], eps );
118  CHECK_REAL_EQUAL( -2.113138721930009, val[1], eps );
119 
120  // Check tags for cell variable ke
121  Tag ke_tag0, ke_tag1;
122  rval = mb.tag_get_handle( "ke0", 1, MB_TYPE_DOUBLE, ke_tag0 );CHECK_ERR( rval );
123  rval = mb.tag_get_handle( "ke1", 1, MB_TYPE_DOUBLE, ke_tag1 );CHECK_ERR( rval );
124 
125  // Get cells (12 pentagons and 630 hexagons)
126  Range cells;
127  rval = mb.get_entities_by_type( 0, MBPOLYGON, cells );CHECK_ERR( rval );
128  CHECK_EQUAL( (size_t)642, cells.size() );
129  CHECK_EQUAL( (size_t)1, cells.psize() );
130 
131  // Check ke tag values on first pentagon and first hexagon
132  EntityHandle cell_ents[] = { cells[0], cells[12] };
133  rval = mb.tag_get_data( ke_tag0, cell_ents, 2, val );CHECK_ERR( rval );
134  CHECK_REAL_EQUAL( 15.001, val[0], eps );
135  CHECK_REAL_EQUAL( 16.013, val[1], eps );
136  rval = mb.tag_get_data( ke_tag1, cell_ents, 2, val );CHECK_ERR( rval );
137  CHECK_REAL_EQUAL( 25.001, val[0], eps );
138  CHECK_REAL_EQUAL( 26.013, val[1], eps );
139  }
140 }
141 
142 void test_read_onevar()
143 {
144  Core moab;
145  Interface& mb = moab;
146 
147  std::string opts;
148  get_options( opts );
149 
150  // Read mesh and read cell variable ke at all timesteps
151  opts += ";VARIABLE=ke";
152  ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
153 
154 #ifdef MOAB_HAVE_MPI
155  ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
156  int procs = pcomm->proc_config().proc_size();
157 #else
158  int procs = 1;
159 #endif
160 
161  // Make check runs this test on one processor
162  if( 1 == procs )
163  {
164  // Check ke tags
165  Tag ke_tag0, ke_tag1;
166  rval = mb.tag_get_handle( "ke0", 1, MB_TYPE_DOUBLE, ke_tag0 );CHECK_ERR( rval );
167  rval = mb.tag_get_handle( "ke1", 1, MB_TYPE_DOUBLE, ke_tag1 );CHECK_ERR( rval );
168 
169  // Get cells (12 pentagons and 630 hexagons)
170  Range cells;
171  rval = mb.get_entities_by_type( 0, MBPOLYGON, cells );CHECK_ERR( rval );
172  CHECK_EQUAL( (size_t)642, cells.size() );
173  CHECK_EQUAL( (size_t)1, cells.psize() );
174 
175  // Check ke tag values on 4 cells: first pentagon, last pentagon,
176  // first hexagon, and last hexagon
177  EntityHandle cell_ents[] = { cells[0], cells[11], cells[12], cells[641] };
178  double ke0_val[4];
179  rval = mb.tag_get_data( ke_tag0, cell_ents, 4, ke0_val );CHECK_ERR( rval );
180  CHECK_REAL_EQUAL( 15.001, ke0_val[0], eps );
181  CHECK_REAL_EQUAL( 15.012, ke0_val[1], eps );
182  CHECK_REAL_EQUAL( 16.013, ke0_val[2], eps );
183  CHECK_REAL_EQUAL( 16.642, ke0_val[3], eps );
184  double ke1_val[4];
185  rval = mb.tag_get_data( ke_tag1, cell_ents, 4, ke1_val );CHECK_ERR( rval );
186  CHECK_REAL_EQUAL( 25.001, ke1_val[0], eps );
187  CHECK_REAL_EQUAL( 25.012, ke1_val[1], eps );
188  CHECK_REAL_EQUAL( 26.013, ke1_val[2], eps );
189  CHECK_REAL_EQUAL( 26.642, ke1_val[3], eps );
190  }
191 }
192 
193 void test_read_onetimestep()
194 {
195  Core moab;
196  Interface& mb = moab;
197 
198  std::string opts;
199  get_options( opts );
200 
201  // Read mesh and read all variables at 2nd timestep
202  opts += ";TIMESTEP=1";
203  ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
204 
205  // Check vorticity tags
206  Tag vorticity_tag0, vorticity_tag1;
207  rval = mb.tag_get_handle( "vorticity0", 1, MB_TYPE_DOUBLE, vorticity_tag0 );
208  // Tag vorticity0 should not exist
209  CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
210  rval = mb.tag_get_handle( "vorticity1", 1, MB_TYPE_DOUBLE, vorticity_tag1 );CHECK_ERR( rval );
211 }
212 
213 void test_read_nomesh()
214 {
215  Core moab;
216  Interface& mb = moab;
217 
218  // Need a file set for nomesh to work right
219  EntityHandle file_set;
220  ErrorCode rval = mb.create_meshset( MESHSET_SET, file_set );CHECK_ERR( rval );
221 
222  std::string orig, opts;
223  get_options( orig );
224 
225  // Read mesh and read all variables at 1st timestep
226  opts = orig + ";TIMESTEP=0";
227  rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
228 
229  // Check u tags
230  Tag u_tag0, u_tag1;
231  rval = mb.tag_get_handle( "u0", 1, MB_TYPE_DOUBLE, u_tag0 );CHECK_ERR( rval );
232  // Tag u1 should not exist
233  rval = mb.tag_get_handle( "u1", 1, MB_TYPE_DOUBLE, u_tag1 );
234  CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
235 
236  // Read all variables at 2nd timestep 0, no need to read mesh
237  opts = orig + ";TIMESTEP=1;NOMESH";
238  rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
239 
240  // Check tag u1 again
241  rval = mb.tag_get_handle( "u1", 1, MB_TYPE_DOUBLE, u_tag1 );
242  // Tag u1 should exist at this time
243  CHECK_ERR( rval );
244 }
245 
246 void test_read_novars()
247 {
248  Core moab;
249  Interface& mb = moab;
250 
251  // Need a file set for nomesh to work right
252  EntityHandle file_set;
253  ErrorCode rval = mb.create_meshset( MESHSET_SET, file_set );CHECK_ERR( rval );
254 
255  std::string orig, opts;
256  get_options( orig );CHECK_ERR( rval );
257 
258  // Read header info only, no mesh, no variables
259  opts = orig + ";NOMESH;VARIABLE=";
260  rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
261 
262  // Read mesh, but still no variables
263  opts = orig + ";VARIABLE=;TIMESTEP=0";
264  rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
265 
266  // Check ke tags
267  Tag ke_tag0, ke_tag1;
268  rval = mb.tag_get_handle( "ke0", 1, MB_TYPE_DOUBLE, ke_tag0 );
269  // Tag ke0 should not exist
270  CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
271  rval = mb.tag_get_handle( "ke1", 1, MB_TYPE_DOUBLE, ke_tag1 );
272  // Tag ke1 should not exist
273  CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
274 
275  // Read ke at 1st timestep, no need to read mesh
276  opts = orig + ";VARIABLE=ke;TIMESTEP=0;NOMESH";
277  rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
278 
279  // Check ke tags again
280  rval = mb.tag_get_handle( "ke0", 1, MB_TYPE_DOUBLE, ke_tag0 );
281  // Tag ke0 should exist at this time
282  CHECK_ERR( rval );
283  // Tag ke1 should still not exist
284  rval = mb.tag_get_handle( "ke1", 1, MB_TYPE_DOUBLE, ke_tag1 );
285  CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
286 
287  // Read ke at 2nd timestep, no need to read mesh
288  opts = orig + ";VARIABLE=ke;TIMESTEP=1;NOMESH";
289  rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
290 
291  // Check tag ke1 again
292  rval = mb.tag_get_handle( "ke1", 1, MB_TYPE_DOUBLE, ke_tag1 );
293  // Tag ke1 should exist at this time
294  CHECK_ERR( rval );
295 
296 #ifdef MOAB_HAVE_MPI
297  ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
298  int procs = pcomm->proc_config().proc_size();
299 #else
300  int procs = 1;
301 #endif
302 
303  // Make check runs this test on one processor
304  if( 1 == procs )
305  {
306  // Get cells (12 pentagons and 630 hexagons)
307  Range cells;
308  rval = mb.get_entities_by_type( file_set, MBPOLYGON, cells );CHECK_ERR( rval );
309  CHECK_EQUAL( (size_t)642, cells.size() );
310  CHECK_EQUAL( (size_t)1, cells.psize() );
311 
312  // Check ke tag values on 4 cells: first pentagon, last pentagon,
313  // first hexagon, and last hexagon
314  EntityHandle cell_ents[] = { cells[0], cells[11], cells[12], cells[641] };
315  double ke0_val[4];
316  rval = mb.tag_get_data( ke_tag0, cell_ents, 4, ke0_val );CHECK_ERR( rval );
317  CHECK_REAL_EQUAL( 15.001, ke0_val[0], eps );
318  CHECK_REAL_EQUAL( 15.012, ke0_val[1], eps );
319  CHECK_REAL_EQUAL( 16.013, ke0_val[2], eps );
320  CHECK_REAL_EQUAL( 16.642, ke0_val[3], eps );
321  double ke1_val[4];
322  rval = mb.tag_get_data( ke_tag1, cell_ents, 4, ke1_val );CHECK_ERR( rval );
323  CHECK_REAL_EQUAL( 25.001, ke1_val[0], eps );
324  CHECK_REAL_EQUAL( 25.012, ke1_val[1], eps );
325  CHECK_REAL_EQUAL( 26.013, ke1_val[2], eps );
326  CHECK_REAL_EQUAL( 26.642, ke1_val[3], eps );
327  }
328 }
329 
330 void test_read_no_mixed_elements()
331 {
332  Core moab;
333  Interface& mb = moab;
334 
335  std::string opts;
336  get_options( opts );
337 
338  // Read mesh with no mixed elements and read all variables at all timesteps
339  opts += ";NO_MIXED_ELEMENTS";
340  ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
341 
342 #ifdef MOAB_HAVE_MPI
343  ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
344  int procs = pcomm->proc_config().proc_size();
345 #else
346  int procs = 1;
347 #endif
348 
349  // Make check runs this test on one processor
350  if( 1 == procs )
351  {
352  // Check ke tags
353  Tag ke_tag0, ke_tag1;
354  rval = mb.tag_get_handle( "ke0", 1, MB_TYPE_DOUBLE, ke_tag0 );CHECK_ERR( rval );
355  rval = mb.tag_get_handle( "ke1", 1, MB_TYPE_DOUBLE, ke_tag1 );CHECK_ERR( rval );
356 
357  // Get cells (12 pentagons and 630 hexagons)
358  Range cells;
359  rval = mb.get_entities_by_type( 0, MBPOLYGON, cells );CHECK_ERR( rval );
360  CHECK_EQUAL( (size_t)642, cells.size() );
361  // Only one group of cells (pentagons are padded to hexagons,
362  // e.g. connectivity [1 2 3 4 5] => [1 2 3 4 5 5])
363  CHECK_EQUAL( (size_t)1, cells.psize() );
364 
365  // Check ke tag values on 4 cells: first pentagon, last pentagon,
366  // first hexagon, and last hexagon
367  EntityHandle cell_ents[] = { cells[0], cells[11], cells[12], cells[641] };
368  double ke0_val[4];
369  rval = mb.tag_get_data( ke_tag0, cell_ents, 4, ke0_val );CHECK_ERR( rval );
370  CHECK_REAL_EQUAL( 15.001, ke0_val[0], eps );
371  CHECK_REAL_EQUAL( 15.012, ke0_val[1], eps );
372  CHECK_REAL_EQUAL( 16.013, ke0_val[2], eps );
373  CHECK_REAL_EQUAL( 16.642, ke0_val[3], eps );
374  double ke1_val[4];
375  rval = mb.tag_get_data( ke_tag1, cell_ents, 4, ke1_val );CHECK_ERR( rval );
376  CHECK_REAL_EQUAL( 25.001, ke1_val[0], eps );
377  CHECK_REAL_EQUAL( 25.012, ke1_val[1], eps );
378  CHECK_REAL_EQUAL( 26.013, ke1_val[2], eps );
379  CHECK_REAL_EQUAL( 26.642, ke1_val[3], eps );
380  }
381 }
382 
383 void test_read_no_edges()
384 {
385  Core moab;
386  Interface& mb = moab;
387 
388  std::string opts;
389  get_options( opts );
390 
391  opts += ";NO_EDGES;VARIABLE=";
392  ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
393 
394  // Get edges
395  Range edges;
396  rval = mb.get_entities_by_type( 0, MBEDGE, edges );CHECK_ERR( rval );
397  CHECK_EQUAL( (size_t)0, edges.size() );
398 }
399 
400 void test_gather_onevar()
401 {
402  Core moab;
403  Interface& mb = moab;
404 
405  EntityHandle file_set;
406  ErrorCode rval = mb.create_meshset( MESHSET_SET, file_set );CHECK_ERR( rval );
407 
408  std::string opts;
409  get_options( opts );
410 
411  // Read cell variable ke and create gather set on processor 0
412  opts += ";VARIABLE=ke;GATHER_SET=0";
413  rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
414 
415 #ifdef MOAB_HAVE_MPI
416  ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
417  int rank = pcomm->proc_config().proc_rank();
418 
419  Range cells, cells_owned;
420  rval = mb.get_entities_by_type( file_set, MBPOLYGON, cells );CHECK_ERR( rval );
421 
422  // Get local owned cells
423  rval = pcomm->filter_pstatus( cells, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &cells_owned );CHECK_ERR( rval );
424 
425  EntityHandle gather_set = 0;
426  if( 0 == rank )
427  {
428  // Get gather set
429  ReadUtilIface* readUtilIface;
430  mb.query_interface( readUtilIface );
431  rval = readUtilIface->get_gather_set( gather_set );CHECK_ERR( rval );
432  assert( gather_set != 0 );
433  }
434 
435  Tag ke_tag0, gid_tag;
436  rval = mb.tag_get_handle( "ke0", 1, MB_TYPE_DOUBLE, ke_tag0, MB_TAG_DENSE );CHECK_ERR( rval );
437 
438  gid_tag = mb.globalId_tag();
439 
440  pcomm->gather_data( cells_owned, ke_tag0, gid_tag, gather_set, 0 );
441 
442  if( 0 == rank )
443  {
444  // Get gather set cells
445  Range gather_set_cells;
446  rval = mb.get_entities_by_type( gather_set, MBPOLYGON, gather_set_cells );CHECK_ERR( rval );
447  CHECK_EQUAL( (size_t)642, gather_set_cells.size() );
448  CHECK_EQUAL( (size_t)1, gather_set_cells.psize() );
449 
450  // Check ke0 tag values on 4 gather set cells: first pentagon, last pentagon,
451  // first hexagon, and last hexagon
452  double ke0_val[4];
453  EntityHandle cell_ents[] = { gather_set_cells[0], gather_set_cells[11], gather_set_cells[12],
454  gather_set_cells[641] };
455  rval = mb.tag_get_data( ke_tag0, cell_ents, 4, ke0_val );CHECK_ERR( rval );
456  CHECK_REAL_EQUAL( 15.001, ke0_val[0], eps );
457  CHECK_REAL_EQUAL( 15.012, ke0_val[1], eps );
458  CHECK_REAL_EQUAL( 16.013, ke0_val[2], eps );
459  CHECK_REAL_EQUAL( 16.642, ke0_val[3], eps );
460  }
461 #endif
462 }
463 
464 void get_options( std::string& opts )
465 {
466 #ifdef MOAB_HAVE_MPI
467  // Use parallel options
468  opts = std::string( ";;PARALLEL=READ_PART;PARTITION_METHOD=TRIVIAL" );
469 #else
470  opts = std::string( ";;" );
471 #endif
472 }