Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
library/FieldUtils/InputModules/InputNek5000.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: InputNek5000.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Reads a Nek5000 checkpoint file.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 #include <fstream>
38 #include <string>
39 #include <vector>
40 using namespace std;
41 
42 #include <boost/algorithm/string.hpp>
43 
44 #include "InputNek5000.h"
45 
46 namespace Nektar
47 {
48 namespace FieldUtils
49 {
50 
51 ModuleKey InputNek5000::m_className[1] = {
53  ModuleKey(eInputModule, "fld5000"), InputNek5000::create,
54  "Reads Nek5000 field file.")
55 };
56 
57 /**
58  * @brief Swap endian ordering of the input variable.
59  */
60 template <typename T>
61 void swap_endian(T &u)
62 {
63  union
64  {
65  T u;
66  unsigned char u8[sizeof(T)];
67  } source, dest;
68 
69  source.u = u;
70 
71  for (size_t k = 0; k < sizeof(T); k++)
72  {
73  dest.u8[k] = source.u8[sizeof(T) - k - 1];
74  }
75 
76  u = dest.u;
77 }
78 
79 /**
80  * @brief Set up InputNek5000 object.
81  *
82  */
83 InputNek5000::InputNek5000(FieldSharedPtr f) : InputModule(f)
84 {
85  m_allowedFiles.insert("fld5000");
86 }
87 
88 /**
89  *
90  */
92 {
93 }
94 
95 /**
96  * @brief Process Nek5000 input file.
97  *
98  * This routine reads a binary-format Nek5000 field file, loads the data into
99  * memory and populates the field definitions to match the data format. Nek5000
100  * is a classic nodal-Lagrangian spectral element code at a single polynomial
101  * order, meaning that the field data are set up according to this structure.
102  *
103  * This module is adapted from the VisIt visualisation software, which supports
104  * a number of Nek5000 inputs.
105  */
106 void InputNek5000::Process(po::variables_map &vm)
107 {
108  if (m_f->m_verbose)
109  {
110  if (m_f->m_comm->TreatAsRankZero())
111  {
112  cout << "Processing Nek5000 field file" << endl;
113  }
114  }
115 
116  string fldending = "fld5000";
117  ifstream file(m_f->m_inputfiles[fldending][0].c_str(), ios::binary);
118 
119  // Header: 132 bytes for binary.
120  vector<char> data(132);
121  file.read(&data[0], 132);
122 
123  // Check header: should be the four characters #std
124  string check(&data[0], 4);
125  string header(&data[4], 128);
126 
127  ASSERTL0(check == "#std", "Unable to read file");
128 
129  // Determine whether we need to byte-swap data: 4-byte float at byte 80
130  // should be 6.54321.
131  bool byteSwap = false;
132  float test;
133 
134  file.read((char *)(&test), 4);
135  if (test > 6.5 && test < 6.6)
136  {
137  byteSwap = false;
138  }
139  else
140  {
141  swap_endian(test);
142  ASSERTL0(test > 6.5 && test < 6.6,
143  "Unable to determine endian-ness of input file");
144  byteSwap = true;
145  }
146 
147  stringstream ss;
148  ss.str(header);
149 
150  int nBytes, nBlocksXYZ[3], nBlocks, nTotBlocks, dir, nDirs, nCycle, nDim;
151  NekDouble time;
152 
153  // Read header information (this is written as ASCII)
154  string remain;
155  ss >> nBytes >> nBlocksXYZ[0] >> nBlocksXYZ[1] >> nBlocksXYZ[2]
156  >> nBlocks >> nTotBlocks >> time >> nCycle >> dir >> nDirs >> remain;
157  boost::trim(remain);
158 
159  nDim = nBlocksXYZ[2] == 1 ? 2 : 3;
160 
161  // Print some basic information for input if in verbose mode.
162  if (m_f->m_verbose)
163  {
164  cout << "Found header information:" << endl;
165  cout << " -- " << (byteSwap ? "" : "do not ") << "need to swap endian"
166  << endl;
167  cout << " -- Data byte size : " << nBytes << endl;
168  cout << " -- Number of xyz blocks : " << nBlocksXYZ[0] << "x"
169  << nBlocksXYZ[1] << "x" << nBlocksXYZ[2] << endl;
170  cout << " -- Blocks in file/total : " << nBlocks << "/" << nTotBlocks
171  << endl;
172  cout << " -- Simulation time : " << time << endl;
173  cout << " -- Number of cycles : " << nCycle << endl;
174  cout << " -- Number of dirs : " << dir << "/" << nDirs << endl;
175  cout << " -- Remaining header : " << remain << endl;
176  }
177 
178  // Major limitation: we don't read out of multiple directories
179  ASSERTL0(nDirs == 1, "Number of directories must be one");
180 
181  // We also don't read partial files.
182  ASSERTL0(nBlocks == nTotBlocks, "Partial field output not supported");
183 
184  // We don't support non-double data
185  ASSERTL0(nBytes == 8, "Data file must contain double-precision data");
186 
187  // Set up a field definition
189  LibUtilities::FieldDefinitions>::AllocateSharedPtr();
190  fielddef->m_shapeType = LibUtilities::eHexahedron;
191  fielddef->m_numHomogeneousDir = 0;
192  fielddef->m_homoStrips = false;
193  fielddef->m_pointsDef = false;
194  fielddef->m_uniOrder = true;
195  fielddef->m_numPointsDef = false;
196 
197  for (int i = 0; i < nDim; ++i)
198  {
199  fielddef->m_basis.push_back(LibUtilities::eGLL_Lagrange);
200  fielddef->m_numModes.push_back(nBlocksXYZ[i]);
201  }
202 
203  // Read element IDs
204  NekUInt32 maxID = 0, minID = numeric_limits<NekUInt32>::max();
205  for (NekUInt32 i = 0; i < nBlocks; ++i)
206  {
207  NekUInt32 blockNum;
208  file.read((char *)&blockNum, 4);
209  if (byteSwap)
210  {
211  swap_endian(blockNum);
212  }
213  fielddef->m_elementIDs.push_back(blockNum-1);
214 
215  maxID = maxID > blockNum ? maxID : blockNum;
216  minID = minID < blockNum ? minID : blockNum;
217  }
218 
219  // Determine how many fields we have to extract
220  size_t blockSize = nBlocksXYZ[0] * nBlocksXYZ[1] * nBlocksXYZ[2];
221  size_t dataSize = blockSize * nBlocks;
222 
223  for (string::size_type i = 0; i < remain.size(); ++i)
224  {
225  switch (remain[i])
226  {
227  case 'U':
228  fielddef->m_fields.push_back("u");
229  fielddef->m_fields.push_back("v");
230  if (nDim == 3)
231  {
232  fielddef->m_fields.push_back("w");
233  }
234  break;
235  case 'P':
236  fielddef->m_fields.push_back("p");
237  break;
238  case 'T':
239  fielddef->m_fields.push_back("T");
240  break;
241  case '1':
242  case '2':
243  case '3':
244  fielddef->m_fields.push_back(string("scalar") + remain[i]);
245  break;
246  case ' ':
247  continue;
248  default:
249  cerr << "Field contains unknown variable: " << remain[i] << endl;
250  abort();
251  }
252  }
253 
254  m_f->m_data.resize(1);
255  m_f->m_data[0].resize(fielddef->m_fields.size() * dataSize);
256 
257  // Now reprocess since different fields need different logic
258  for (size_t i = 0, cnt = 0; i < remain.size(); ++i)
259  {
260  switch (remain[i])
261  {
262  case 'U':
263  {
264  size_t cntVel[3] = {
265  cnt, cnt + dataSize, cnt + 2*dataSize
266  };
267 
268  for (size_t j = 0; j < nBlocks; ++j)
269  {
270  for (size_t k = 0; k < nDim; ++k)
271  {
272  file.read(
273  (char *)&m_f->m_data[0][cntVel[k]],
274  blockSize * sizeof(NekDouble));
275  cntVel[k] += blockSize;
276  }
277  }
278 
279  cnt += nDim * dataSize;
280  break;
281  }
282  case 'P':
283  {
284  file.read(
285  (char *)&m_f->m_data[0][cnt],
286  dataSize * sizeof(NekDouble));
287  cnt += dataSize;
288  break;
289  }
290  case '1':
291  case '2':
292  case '3':
293  {
294  file.read(
295  (char *)&m_f->m_data[0][cnt],
296  dataSize * sizeof(NekDouble));
297  cnt += dataSize;
298  break;
299  }
300  case ' ':
301  continue;
302  }
303  }
304 
305  m_f->m_fielddef.push_back(fielddef);
306 }
307 }
308 }
void swap_endian(T &u)
Swap endian ordering of the input variable.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
boost::shared_ptr< FieldDefinitions > FieldDefinitionsSharedPtr
Definition: FieldIO.h:181
STL namespace.
pair< ModuleType, string > ModuleKey
virtual void Process(po::variables_map &vm)
Process Nek5000 input file.
Abstract base class for input modules.
Metadata that describes the storage properties of field output.
Definition: FieldIO.h:104
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:767
double NekDouble
#define dest(otri, vertexptr)
boost::uint32_t NekUInt32
Lagrange for SEM basis .
Definition: BasisType.h:53
ModuleFactory & GetModuleFactory()
FieldSharedPtr m_f
Field object.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215