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