Nektar++
OutputFileBase.cpp
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////////////
2//
3// File: OutputFileBase.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// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Base class for outputting to a file
32//
33////////////////////////////////////////////////////////////////////////////////
34
35#include <set>
36#include <string>
37using namespace std;
38
39#include "OutputFileBase.h"
41#include <boost/format.hpp>
42#include <iomanip>
43
44namespace Nektar::FieldUtils
45{
46
48{
49 m_requireEquiSpaced = false;
50 m_prohibitWrite = false;
51 m_config["writemultiplefiles"] = ConfigOption(
52 true, "0",
53 "Write multiple files in parallel or when using nparts option");
54}
55
57{
58}
59
60void OutputFileBase::v_Process(po::variables_map &vm)
61{
62 m_f->SetUpExp(vm);
63
64 string filename = m_config["outfile"].as<string>();
65
66 if (filename == "")
67 {
68 m_prohibitWrite = true;
69 }
70
71 if (m_f->m_fieldPts != LibUtilities::NullPtsField)
72 {
73 ASSERTL0(!m_f->m_writeBndFld, "Boundary can't be obtained from pts.");
74 if (WriteFile(filename, vm))
75 {
77
78 if (vm.count("error"))
79 {
81 }
82 }
83 }
84 else if (m_f->m_exp.size())
85 {
86 // reset expansion definition to use equispaced points if required.
87 if (m_requireEquiSpaced && (vm.count("no-equispaced") == 0) &&
88 m_f->m_exp[0]->GetNumElmts() != 0 && !m_equispacedSetup)
89 {
91 }
92 if (m_f->m_writeBndFld)
93 {
94 if (m_f->m_verbose &&
95 m_f->m_comm->GetSpaceComm()->TreatAsRankZero())
96 {
97 cout << "\t" << GetModuleName()
98 << ": Writing boundary file(s): ";
99 for (int i = 0; i < m_f->m_bndRegionsToWrite.size(); ++i)
100 {
101 cout << m_f->m_bndRegionsToWrite[i];
102 if (i < m_f->m_bndRegionsToWrite.size() - 1)
103 {
104 cout << ", ";
105 }
106 }
107 cout << endl;
108 }
109
110 int nfields = m_f->m_variables.size();
111 int normdim = m_f->m_graph->GetMeshDimension();
112
113 // Prepare for normals output
114 if (m_f->m_addNormals)
115 {
116 // Prepare for creating expansions for normals
117 m_f->m_exp.resize(nfields + normdim);
118
119 // Include normal name in m_variables
120 string normstr[3] = {"Norm_x", "Norm_y", "Norm_z"};
121 for (int j = 0; j < normdim; ++j)
122 {
123 m_f->m_exp[nfields + j] =
124 m_f->AppendExpList(m_f->m_numHomogeneousDir);
125 m_f->m_variables.push_back(normstr[j]);
126 }
127 }
128
129 // Move m_exp to a new expansion vector
130 vector<MultiRegions::ExpListSharedPtr> exp(m_f->m_exp.size());
131 exp.swap(m_f->m_exp);
132
134 BndExp(exp.size());
135 for (int i = 0; i < exp.size(); ++i)
136 {
137 BndExp[i] = exp[i]->GetBndCondExpansions();
138 }
139
140 // get hold of partition boundary regions so we can match it to
141 // desired region extraction
143 exp[0]->GetGraph());
145 bcs.GetBoundaryRegions();
146 map<int, int> BndRegionMap;
147 map<int, LibUtilities::CommSharedPtr> BndRegionComm;
148 int cnt = 0;
149 for (auto &breg_it : bregions)
150 {
151 BndRegionMap[breg_it.first] = cnt++;
152 BndRegionComm[breg_it.first] =
153 bcs.GetBoundaryCommunicators()[breg_it.first];
154 }
155
156 // find ending of output file and insert _b1, _b2
157 int dot = filename.find_last_of('.') + 1;
158 string ext = filename.substr(dot, filename.length() - dot);
159 string name = filename.substr(0, dot - 1);
160
161 // Store temporary communicator
162 LibUtilities::CommSharedPtr tmpComm = m_f->m_comm;
163
164 for (int i = 0; i < m_f->m_bndRegionsToWrite.size(); ++i)
165 {
166 string outname = name + "_b" +
167 std::to_string(m_f->m_bndRegionsToWrite[i]) +
168 "." + ext;
169
170 if (!WriteFile(outname, vm))
171 {
172 continue;
173 }
174 RegisterConfig("outfile", outname);
175
176 if (BndRegionMap.count(m_f->m_bndRegionsToWrite[i]) == 1)
177 {
178 m_f->m_comm = BndRegionComm[m_f->m_bndRegionsToWrite[i]];
179
180 int Border = BndRegionMap[m_f->m_bndRegionsToWrite[i]];
181
182 // set up m_exp to point to boundary expansion
183 for (int j = 0; j < exp.size(); ++j)
184 {
185 m_f->m_exp[j] = BndExp[j][Border];
186 }
187
188 for (int j = 0; j < exp.size(); ++j)
189 {
190 m_f->m_exp[j] = BndExp[j][Border];
191 m_f->m_exp[j]->BwdTrans(m_f->m_exp[j]->GetCoeffs(),
192 m_f->m_exp[j]->UpdatePhys());
193 }
194
195 if (m_f->m_addNormals)
196 {
197 // Get normals
199 exp[0]->GetBoundaryNormals(Border, NormPhys);
200
201 // add normal coefficients to expansions
202 for (int j = 0; j < normdim; ++j)
203 {
204 m_f->m_exp[nfields + j] =
205 BndExp[nfields + j][Border];
207 m_f->m_exp[nfields + j]->GetTotPoints(),
208 NormPhys[j], 1,
209 m_f->m_exp[nfields + j]->UpdatePhys(), 1);
210 m_f->m_exp[nfields + j]->FwdTransLocalElmt(
211 m_f->m_exp[nfields + j]->GetPhys(),
212 m_f->m_exp[nfields + j]->UpdateCoeffs());
213 }
214 }
215 v_OutputFromExp(vm);
216 // output error for regression checking.
217 if (vm.count("error"))
218 {
220 }
221
222 // Reset communicator
223 m_f->m_comm = tmpComm;
224 }
225
226 // put outfile back to filename in case of nparts option
227 RegisterConfig("outfile", filename);
228 }
229 // Restore m_exp
230 exp.swap(m_f->m_exp);
231 }
232 else
233 {
234 if (WriteFile(filename, vm))
235 {
236 v_OutputFromExp(vm);
237 // output error for regression checking.
238 if (vm.count("error"))
239 {
241 }
242 }
243 }
244 }
245 else if (m_f->m_data.size())
246 {
247 ASSERTL0(!m_f->m_writeBndFld, "Boundary extraction requires xml file.");
248 if (WriteFile(filename, vm))
249 {
251 }
252 }
253}
254
255bool OutputFileBase::WriteFile(std::string &filename, po::variables_map &vm)
256{
257 if (m_prohibitWrite)
258 {
259 return true;
260 }
261
262 // Get path to file. If procid was defined, get the full name to avoid
263 // checking files from other partitions
264 fs::path outFile;
265 if (vm.count("nparts"))
266 {
267 outFile = GetFullOutName(filename, vm);
268 }
269 else
270 {
271 outFile = GetPath(filename, vm);
272 }
273
275 if (m_f->m_comm)
276 {
277 comm = m_f->m_comm;
278 }
279 else
280 {
281 comm = LibUtilities::GetCommFactory().CreateInstance("Serial", 0, 0);
282 }
283
284 int count = fs::exists(outFile) ? 1 : 0;
285 comm->GetSpaceComm()->AllReduce(count, LibUtilities::ReduceSum);
286
287 int writeFile = 1;
288 if (count && (vm.count("force-output") == 0))
289 {
290 if (vm.count("nparts") == 0) // do not do check if --nparts is enabled.
291 {
292 writeFile = 0; // set to zero for reduce all to be correct.
293
294 if (comm->GetSpaceComm()->TreatAsRankZero())
295 {
296 string answer;
297 cout << "Did you wish to overwrite " << outFile << " (y/n)? ";
298 getline(cin, answer);
299 if (answer.compare("y") == 0)
300 {
301 writeFile = 1;
302 }
303 else
304 {
305 cout << "Not writing file '" << filename
306 << "' because it already exists" << endl;
307 }
308 }
309 comm->GetSpaceComm()->AllReduce(writeFile, LibUtilities::ReduceSum);
310 }
311 }
312 return (writeFile == 0) ? false : true;
313}
314
315void OutputFileBase::ConvertExpToEquispaced(po::variables_map &vm)
316{
317 // Information to create new expansion
318 int numFields = m_f->m_exp.size();
319 m_f->m_fielddef = m_f->m_exp[0]->GetFieldDefinitions();
320
321 // Set points to equispaced
322 int nPointsNew = 0;
323 if (vm.count("output-points"))
324 {
325 nPointsNew = vm["output-points"].as<int>();
326 }
327 m_f->m_graph->SetExpansionInfoToEvenlySpacedPoints(nPointsNew);
328
329 // Save original expansion
330 vector<MultiRegions::ExpListSharedPtr> expOld = m_f->m_exp;
331
332 // Create new expansion
333 m_f->m_exp[0] = m_f->SetUpFirstExpList(m_f->m_numHomogeneousDir, true);
334 for (int i = 1; i < numFields; ++i)
335 {
336 m_f->m_exp[i] = m_f->AppendExpList(m_f->m_numHomogeneousDir);
337 }
338
339 // Extract result to new expansion
340 for (int i = 0; i < numFields; ++i)
341 {
342 m_f->m_exp[i]->ExtractCoeffsToCoeffs(expOld[i], expOld[i]->GetCoeffs(),
343 m_f->m_exp[i]->UpdateCoeffs());
344 m_f->m_exp[i]->BwdTrans(m_f->m_exp[i]->GetCoeffs(),
345 m_f->m_exp[i]->UpdatePhys());
346 }
347
348 // Extract boundary expansion if needed
349 if (m_f->m_writeBndFld)
350 {
353
354 for (int i = 0; i < numFields; ++i)
355 {
356 BndExpOld = expOld[i]->GetBndCondExpansions();
357 for (int j = 0; j < BndExpOld.size(); ++j)
358 {
359 BndExp = m_f->m_exp[i]->UpdateBndCondExpansion(j);
360 BndExp->ExtractCoeffsToCoeffs(BndExpOld[j],
361 BndExpOld[j]->GetCoeffs(),
362 BndExp->UpdateCoeffs());
363 }
364 }
365 }
366
367 m_f->m_fielddef = std::vector<LibUtilities::FieldDefinitionsSharedPtr>();
368
369 // Make a note so that we know that the fielddefs have now been setup.
370 m_equispacedSetup = true;
371}
372
374{
375 int coordim = m_f->m_fieldPts->GetDim();
376 std::string coordVars[] = {"x", "y", "z"};
377
378 vector<string> variables = m_f->m_variables;
379 variables.insert(variables.begin(), coordVars, coordVars + coordim);
380 // Get fields and coordinates
381 Array<OneD, Array<OneD, NekDouble>> fields(variables.size());
382
383 // We can just grab everything from points. This should be a
384 // reference, not a copy.
385 m_f->m_fieldPts->GetPts(fields);
386 for (int i = 0; i < fields.size(); ++i)
387 {
388 // calculate L2 and Linf value
389 int npts = fields[i].size();
390
391 NekDouble l2err = 0.0;
392 NekDouble linferr = 0.0;
393 for (int j = 0; j < npts; ++j)
394 {
395 l2err += fields[i][j] * fields[i][j];
396 linferr = max(linferr, fabs(fields[i][j]));
397 }
398
399 m_f->m_comm->GetSpaceComm()->AllReduce(l2err, LibUtilities::ReduceSum);
400 m_f->m_comm->GetSpaceComm()->AllReduce(npts, LibUtilities::ReduceSum);
401 m_f->m_comm->GetSpaceComm()->AllReduce(linferr,
403
404 l2err /= npts;
405 l2err = sqrt(l2err);
406
407 if (m_f->m_comm->GetSpaceComm()->TreatAsRankZero())
408 {
409 cout << "L 2 error (variable " << variables[i] << ") : " << l2err
410 << endl;
411
412 cout << "L inf error (variable " << variables[i]
413 << ") : " << linferr << endl;
414 }
415 }
416}
417
419{
420 int coordim =
421 m_f->m_exp[0]->GetExp(0)->GetCoordim() + m_f->m_numHomogeneousDir;
422 int totpoints = m_f->m_exp[0]->GetTotPoints();
423 std::string coordVars[] = {"x", "y", "z"};
424
425 // Set up storage for coordinates
426 Array<OneD, Array<OneD, NekDouble>> coords(coordim);
427 for (int i = 0; i < coordim; ++i)
428 {
429 coords[i] = Array<OneD, NekDouble>(totpoints);
430 }
431
432 // Get coordinates
433 if (coordim == 1)
434 {
435 m_f->m_exp[0]->GetCoords(coords[0]);
436 }
437 else if (coordim == 2)
438 {
439 m_f->m_exp[0]->GetCoords(coords[0], coords[1]);
440 }
441 else
442 {
443 m_f->m_exp[0]->GetCoords(coords[0], coords[1], coords[2]);
444 }
445
446 for (int j = 0; j < coordim; ++j)
447 {
448 NekDouble l2err = m_f->m_exp[0]->L2(coords[j]);
449 NekDouble linferr = m_f->m_exp[0]->Linf(coords[j]);
450
451 if (m_f->m_comm->GetSpaceComm()->TreatAsRankZero())
452 {
453 cout << "L 2 error (variable " << coordVars[j] << ") : " << l2err
454 << endl;
455
456 cout << "L inf error (variable " << coordVars[j]
457 << ") : " << linferr << endl;
458 }
459 }
460
461 for (int j = 0; j < m_f->m_exp.size(); ++j)
462 {
463 NekDouble l2err = m_f->m_exp[j]->L2(m_f->m_exp[j]->GetPhys());
464 NekDouble linferr = m_f->m_exp[j]->Linf(m_f->m_exp[j]->GetPhys());
465
466 if (m_f->m_comm->GetSpaceComm()->TreatAsRankZero() &&
467 m_f->m_variables.size() > 0)
468 {
469 cout << "L 2 error (variable " << m_f->m_variables[j]
470 << ") : " << l2err << endl;
471
472 cout << "L inf error (variable " << m_f->m_variables[j]
473 << ") : " << linferr << endl;
474 }
475 }
476}
477
478} // namespace Nektar::FieldUtils
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
FIELD_UTILS_EXPORT void RegisterConfig(std::string key, std::string value="")
Register a configuration option with a module.
Definition: Module.cpp:100
std::string GetModuleName()
Definition: Module.h:200
FieldSharedPtr m_f
Field object.
Definition: Module.h:239
std::map< std::string, ConfigOption > m_config
List of configuration values.
Definition: Module.h:272
virtual void v_OutputFromExp(po::variables_map &vm)=0
Write from m_exp to output file.
bool WriteFile(std::string &filename, po::variables_map &vm)
fs::path GetFullOutName(std::string &filename, po::variables_map &vm)
void ConvertExpToEquispaced(po::variables_map &vm)
virtual void v_OutputFromData(po::variables_map &vm)=0
Write from data to output file.
virtual void v_OutputFromPts(po::variables_map &vm)=0
Write from pts to output file.
fs::path GetPath(std::string &filename, po::variables_map &vm)
void v_Process(po::variables_map &vm) override
Write fld to output file.
Abstract base class for output modules.
Definition: Module.h:316
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
std::map< int, LibUtilities::CommSharedPtr > GetBoundaryCommunicators() const
Definition: Conditions.h:261
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:235
std::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:1026
static PtsFieldSharedPtr NullPtsField
Definition: PtsField.h:185
CommFactory & GetCommFactory()
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:211
double NekDouble
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
STL namespace.
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285
Represents a command-line configuration option.
Definition: Module.h:129