Nektar++
SessionFunction.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: SessionFunction.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2017 Kilian Lackhove
10 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
11 // Department of Aeronautics, Imperial College London (UK), and Scientific
12 // Computing and Imaging Institute, University of Utah (USA).
13 //
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: Session Function
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
38 
40 
41 #include <boost/format.hpp>
42 
43 using namespace std;
44 
45 namespace Nektar
46 {
47 namespace SolverUtils
48 {
49 
50 /**
51  * Representation of a FUNCTION defined in the session xml file.
52  *
53  * @param session The session where the function was defined.
54  * @param field The field the function is defined on.
55  * @param functionName The name of the function.
56  * @param toCache Store the evaluated function for later use.
57  */
58 SessionFunction::SessionFunction(
60  const MultiRegions::ExpListSharedPtr &field, std::string functionName,
61  bool toCache)
62  : m_session(session), m_field(field), m_name(functionName),
63  m_toCache(toCache)
64 {
65  ASSERTL0(m_session->DefinesFunction(m_name),
66  "Function '" + m_name + "' does not exist.");
67 }
68 
69 /**
70  * Evaluates a function defined in the xml session file at each quadrature
71  * point.
72  *
73  * @param pArray The array into which to write the values.
74  * @param pTime The time at which to evaluate the function.
75  * @param domain The domain to evaluate the function in.
76  */
78  const NekDouble pTime, const int domain)
79 {
80  std::vector<std::string> vFieldNames = m_session->GetVariables();
81 
82  for (int i = 0; i < vFieldNames.size(); i++)
83  {
84  Evaluate(vFieldNames[i], pArray[i], pTime, domain);
85  }
86 }
87 
88 /**
89  * Evaluates a function defined in the xml session file at each quadrature
90  * point.
91  *
92  * @param pFieldNames The names of the fields to evaluate the function for.
93  * @param pArray The array into which to write the values.
94  * @param pTime The time at which to evaluate the function.
95  * @param domain The domain to evaluate the function in.
96  */
97 void SessionFunction::Evaluate(std::vector<std::string> pFieldNames,
99  const NekDouble &pTime, const int domain)
100 {
101  ASSERTL1(pFieldNames.size() == pArray.size(),
102  "Function '" + m_name +
103  "' variable list size mismatch with array storage.");
104 
105  for (int i = 0; i < pFieldNames.size(); i++)
106  {
107  Evaluate(pFieldNames[i], pArray[i], pTime, domain);
108  }
109 }
110 
111 /**
112  * Evaluates a function defined in the xml session file at each quadrature
113  * point.
114  *
115  * @param pFieldNames The names of the fields to evaluate the function for.
116  * @param pFields The fields into which to write the values.
117  * @param pTime The time at which to evaluate the function.
118  * @param domain The domain to evaluate the function in.
119  */
121  std::vector<std::string> pFieldNames,
123  const NekDouble &pTime, const int domain)
124 {
125  ASSERTL0(pFieldNames.size() == pFields.size(),
126  "Field list / name list size mismatch.");
127 
128  for (int i = 0; i < pFieldNames.size(); i++)
129  {
130  Evaluate(pFieldNames[i], pFields[i]->UpdatePhys(), pTime, domain);
131  pFields[i]->FwdTransLocalElmt(pFields[i]->GetPhys(),
132  pFields[i]->UpdateCoeffs());
133  }
134 }
135 
136 /**
137  * Evaluates a function defined in the xml session file at each quadrature
138  * point.
139  *
140  * @param pFieldName The name of the field to evaluate the function for.
141  * @param pArray The array into which to write the values.
142  * @param pTime The time at which to evaluate the function.
143  * @param domain The domain to evaluate the function in.
144  */
145 void SessionFunction::Evaluate(std::string pFieldName,
146  Array<OneD, NekDouble> &pArray,
147  const NekDouble &pTime, const int domain)
148 {
150  m_session->GetFunctionType(m_name, pFieldName, domain);
151 
152  unsigned int nq = m_field->GetNpoints();
153 
154  std::pair<std::string, int> key(pFieldName, domain);
155  // sorry
156  if ((m_arrays.find(key) != m_arrays.end()) &&
158  ((m_lastCached.find(key) != m_lastCached.end()) &&
159  (pTime - m_lastCached[key] < NekConstants::kNekZeroTol))))
160  {
161  // found cached field
162  if (pArray.size() < nq)
163  {
164  pArray = Array<OneD, NekDouble>(nq);
165  }
166  Vmath::Vcopy(nq, m_arrays[key], 1, pArray, 1);
167 
168  return;
169  }
170 
172  {
173  EvaluateExp(pFieldName, pArray, pTime, domain);
174  }
175  else if (vType == LibUtilities::eFunctionTypeFile ||
177  {
178  std::string filename =
179  m_session->GetFunctionFilename(m_name, pFieldName, domain);
180 
181  if (boost::filesystem::path(filename).extension() == ".pts" ||
182  boost::filesystem::path(filename).extension() == ".csv")
183  {
184  EvaluatePts(pFieldName, pArray, pTime, domain);
185  }
186  else
187  {
188  EvaluateFld(pFieldName, pArray, pTime, domain);
189  }
190  }
191  else
192  {
193  ASSERTL0(false, "unknown eFunctionType");
194  }
195 
196  if (m_toCache)
197  {
198  m_arrays[key] = Array<OneD, NekDouble>(nq);
199  Vmath::Vcopy(nq, pArray, 1, m_arrays[key], 1);
200 
201  m_lastCached[key] = pTime;
202  }
203 }
204 
205 /**
206  * @brief Provide a description of a function for a given field name.
207  *
208  * @param pFieldName Field name.
209  * @param domain The domain to evaluate the function in.
210  */
211 std::string SessionFunction::Describe(std::string pFieldName, const int domain)
212 {
213  std::string retVal;
214 
216  m_session->GetFunctionType(m_name, pFieldName, domain);
218  {
220  m_session->GetFunction(m_name, pFieldName, domain);
221  retVal = ffunc->GetExpression();
222  }
223  else if (vType == LibUtilities::eFunctionTypeFile ||
225  {
226  std::string filename =
227  m_session->GetFunctionFilename(m_name, pFieldName, domain);
228  retVal = "from file " + filename;
229  }
230  else
231  {
232  ASSERTL0(false, "unknown eFunctionType");
233  }
234 
235  return retVal;
236 }
237 
238 /**
239  * Evaluates a function from expression
240  *
241  * @param pFieldName The name of the field to evaluate the function for.
242  * @param pArray The array into which to write the values.
243  * @param pTime The time at which to evaluate the function.
244  * @param domain The domain to evaluate the function in.
245  */
246 void SessionFunction::EvaluateExp(string pFieldName,
247  Array<OneD, NekDouble> &pArray,
248  const NekDouble &pTime, const int domain)
249 {
250  unsigned int nq = m_field->GetNpoints();
251  if (pArray.size() < nq)
252  {
253  pArray = Array<OneD, NekDouble>(nq);
254  }
255 
256  Array<OneD, NekDouble> x0(nq);
257  Array<OneD, NekDouble> x1(nq);
258  Array<OneD, NekDouble> x2(nq);
259 
260  // Get the coordinates (assuming all fields have the same
261  // discretisation)
262  m_field->GetCoords(x0, x1, x2);
264  m_session->GetFunction(m_name, pFieldName, domain);
265 
266  ffunc->Evaluate(x0, x1, x2, pTime, pArray);
267 }
268 
269 /**
270  * Evaluates a function from fld file
271  *
272  * @param pFieldName The name of the field to evaluate the function for.
273  * @param pArray The array into which to write the values.
274  * @param pTime The time at which to evaluate the function.
275  * @param domain The domain to evaluate the function in.
276  */
277 void SessionFunction::EvaluateFld(string pFieldName,
278  Array<OneD, NekDouble> &pArray,
279  const NekDouble &pTime, const int domain)
280 {
281  unsigned int nq = m_field->GetNpoints();
282  if (pArray.size() < nq)
283  {
284  pArray = Array<OneD, NekDouble>(nq);
285  }
286 
287  std::string filename =
288  m_session->GetFunctionFilename(m_name, pFieldName, domain);
289  std::string fileVar =
290  m_session->GetFunctionFilenameVariable(m_name, pFieldName, domain);
291 
292  if (fileVar.length() == 0)
293  {
294  fileVar = pFieldName;
295  }
296 
297  // In case of eFunctionTypeTransientFile, generate filename from
298  // format string
300  m_session->GetFunctionType(m_name, pFieldName, domain);
302  {
303  try
304  {
305 #if (defined _WIN32 && _MSC_VER < 1900)
306  // We need this to make sure boost::format has always
307  // two digits in the exponents of Scientific notation.
308  unsigned int old_exponent_format;
309  old_exponent_format = _set_output_format(_TWO_DIGIT_EXPONENT);
310  filename = boost::str(boost::format(filename) % pTime);
311  _set_output_format(old_exponent_format);
312 #else
313  filename = boost::str(boost::format(filename) % pTime);
314 #endif
315  }
316  catch (...)
317  {
318  ASSERTL0(false, "Invalid Filename in function \"" + m_name +
319  "\", variable \"" + fileVar + "\"")
320  }
321  }
322 
323  // Define list of global element ids
324  int numexp = m_field->GetExpSize();
325  Array<OneD, int> ElementGIDs(numexp);
326  for (int i = 0; i < numexp; ++i)
327  {
328  ElementGIDs[i] = m_field->GetExp(i)->GetGeom()->GetGlobalID();
329  }
330 
331  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
332  std::vector<std::vector<NekDouble>> FieldData;
335  fldIO->Import(filename, FieldDef, FieldData,
337 
338  int idx = -1;
339  Array<OneD, NekDouble> vCoeffs(m_field->GetNcoeffs(), 0.0);
340  // Loop over all the expansions
341  for (int i = 0; i < FieldDef.size(); ++i)
342  {
343  // Find the index of the required field in the
344  // expansion segment
345  for (int j = 0; j < FieldDef[i]->m_fields.size(); ++j)
346  {
347  if (FieldDef[i]->m_fields[j] == fileVar)
348  {
349  idx = j;
350  }
351  }
352 
353  if (idx >= 0)
354  {
355  m_field->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
356  FieldDef[i]->m_fields[idx], vCoeffs);
357  }
358  else
359  {
360  cout << "Field " + fileVar + " not found." << endl;
361  }
362  }
363 
364  m_field->BwdTrans(vCoeffs, pArray);
365 }
366 
367 /**
368  * Evaluates a function from pts file
369  *
370  * @param pFieldName The name of the field to evaluate the function for.
371  * @param pArray The array into which to write the values.
372  * @param pTime The time at which to evaluate the function.
373  * @param domain The domain to evaluate the function in.
374  */
375 void SessionFunction::EvaluatePts(string pFieldName,
376  Array<OneD, NekDouble> &pArray,
377  const NekDouble &pTime, const int domain)
378 {
379  unsigned int nq = m_field->GetNpoints();
380  if (pArray.size() < nq)
381  {
382  pArray = Array<OneD, NekDouble>(nq);
383  }
384 
385  std::string filename =
386  m_session->GetFunctionFilename(m_name, pFieldName, domain);
387  std::string fileVar =
388  m_session->GetFunctionFilenameVariable(m_name, pFieldName, domain);
389 
390  if (fileVar.length() == 0)
391  {
392  fileVar = pFieldName;
393  }
394 
395  // In case of eFunctionTypeTransientFile, generate filename from
396  // format string
398  m_session->GetFunctionType(m_name, pFieldName, domain);
400  {
401  try
402  {
403 #if (defined _WIN32 && _MSC_VER < 1900)
404  // We need this to make sure boost::format has always
405  // two digits in the exponents of Scientific notation.
406  unsigned int old_exponent_format;
407  old_exponent_format = _set_output_format(_TWO_DIGIT_EXPONENT);
408  filename = boost::str(boost::format(filename) % pTime);
409  _set_output_format(old_exponent_format);
410 #else
411  filename = boost::str(boost::format(filename) % pTime);
412 #endif
413  }
414  catch (...)
415  {
416  ASSERTL0(false, "Invalid Filename in function \"" + m_name +
417  "\", variable \"" + fileVar + "\"")
418  }
419  }
420 
422  // check if we already loaded this file. For transient files,
423  // funcFilename != filename so we can make sure we only keep the
424  // latest pts field per funcFilename.
425  std::string funcFilename =
426  m_session->GetFunctionFilename(m_name, pFieldName, domain);
427 
429  if (boost::filesystem::path(filename).extension() == ".pts")
430  {
431  LibUtilities::PtsIO ptsIO(m_session->GetComm());
432  ptsIO.Import(filename, inPts);
433  }
434  else if (boost::filesystem::path(filename).extension() == ".csv")
435  {
436  LibUtilities::CsvIO csvIO(m_session->GetComm());
437  csvIO.Import(filename, inPts);
438  }
439  else
440  {
441  ASSERTL1(false, "Unsupported file type");
442  }
443 
444  Array<OneD, Array<OneD, NekDouble>> pts(inPts->GetDim() +
445  inPts->GetNFields());
446  for (int i = 0; i < inPts->GetDim() + inPts->GetNFields(); ++i)
447  {
448  pts[i] = Array<OneD, NekDouble>(nq);
449  }
450  if (inPts->GetDim() == 1)
451  {
452  m_field->GetCoords(pts[0]);
453  }
454  else if (inPts->GetDim() == 2)
455  {
456  m_field->GetCoords(pts[0], pts[1]);
457  }
458  else if (inPts->GetDim() == 3)
459  {
460  m_field->GetCoords(pts[0], pts[1], pts[2]);
461  }
463  inPts->GetDim(), inPts->GetFieldNames(), pts);
464 
466  if (m_interpolators.find(funcFilename) != m_interpolators.end())
467  {
468  interp = m_interpolators[funcFilename];
469  }
470  else
471  {
473  if (m_session->GetComm()->GetRank() == 0)
474  {
476  this);
477  }
478  interp.CalcWeights(inPts, outPts);
479  if (m_session->GetComm()->GetRank() == 0)
480  {
481  cout << endl;
482  if (m_session->DefinesCmdLineArgument("verbose"))
483  {
484  interp.PrintStatistics();
485  }
486  }
487  }
488 
489  if (m_toCache)
490  {
491  m_interpolators[funcFilename] = interp;
492  }
493 
494  // TODO: only interpolate the field we actually want
495  interp.Interpolate(inPts, outPts);
496 
497  int fieldInd;
498  vector<string> fieldNames = outPts->GetFieldNames();
499  for (fieldInd = 0; fieldInd < fieldNames.size(); ++fieldInd)
500  {
501  if (outPts->GetFieldName(fieldInd) == fileVar)
502  {
503  break;
504  }
505  }
506  ASSERTL0(fieldInd != fieldNames.size(), "field not found");
507 
508  pArray = outPts->GetPts(fieldInd + outPts->GetDim());
509 }
510 
511 // end of namespaces
512 } // namespace SolverUtils
513 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
A class that contains algorithms for interpolation between pts fields, expansions and different meshe...
FIELD_UTILS_EXPORT void Interpolate(const std::vector< MultiRegions::ExpListSharedPtr > expInField, std::vector< MultiRegions::ExpListSharedPtr > &expOutField, NekDouble def_value=0.0)
Interpolate from an expansion to an expansion.
static std::shared_ptr< FieldIO > CreateForFile(const LibUtilities::SessionReaderSharedPtr session, const std::string &filename)
Construct a FieldIO object for a given input filename.
Definition: FieldIO.cpp:227
void PrintStatistics()
Returns if the weights have already been computed.
void SetProgressCallback(FuncPointerT func, ObjectPointerT obj)
sets a callback funtion which gets called every time the interpolation progresses
void CalcWeights(const LibUtilities::PtsFieldSharedPtr ptsInField, LibUtilities::PtsFieldSharedPtr &ptsOutField, bool reuseTree=false)
Compute interpolation weights without doing any interpolation.
void Import(const std::string &inFile, PtsFieldSharedPtr &ptsField, FieldMetaDataMap &fieldmetadatamap=NullFieldMetaDataMap, DomainRangeShPtr &Range=NullDomainRangeShPtr)
Import a pts field from file.
Definition: PtsIO.cpp:68
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
MultiRegions::ExpListSharedPtr m_field
The expansion we want to evaluate this function for.
bool m_toCache
Store resulting arrays (and interpolators)
std::map< std::pair< std::string, int >, Array< OneD, NekDouble > > m_arrays
Cached result arrays.
SOLVER_UTILS_EXPORT void EvaluateFld(std::string pFieldName, Array< OneD, NekDouble > &pArray, const NekDouble &pTime=0.0, const int domain=0)
SOLVER_UTILS_EXPORT void EvaluatePts(std::string pFieldName, Array< OneD, NekDouble > &pArray, const NekDouble &pTime=0.0, const int domain=0)
Evaluates a function from pts file.
SOLVER_UTILS_EXPORT void Evaluate(Array< OneD, Array< OneD, NekDouble >> &pArray, const NekDouble pTime=0.0, const int domain=0)
Evaluates a function defined in the xml session file at each quadrature point.
SOLVER_UTILS_EXPORT void PrintProgressbar(const int position, const int goal) const
std::map< std::pair< std::string, int >, NekDouble > m_lastCached
Last time the cache for this variable & domain combo was updated.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
std::map< std::string, FieldUtils::Interpolator > m_interpolators
Interpolator for pts file input for a variable & domain combination.
SOLVER_UTILS_EXPORT std::string Describe(std::string pFieldName, const int domain=0)
Provide a description of a function for a given field name.
SOLVER_UTILS_EXPORT void EvaluateExp(std::string pFieldName, Array< OneD, NekDouble > &pArray, const NekDouble &pTime=0.0, const int domain=0)
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:301
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:130
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:53
std::shared_ptr< PtsField > PtsFieldSharedPtr
Definition: PtsField.h:190
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
static const NekDouble kNekZeroTol
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
double NekDouble
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255