Commit 25790739 authored by Dave Moxey's avatar Dave Moxey

Merge branch 'tidy/SessionFunction' into 'master'

Tidy: Move EquationSystem::EvaluateFunction and Forcing::EvaluateFunction into separate class

See merge request !759
parents ec641e68 fe571186
......@@ -10,6 +10,8 @@ v5.0.0
- Fix ThridpartyCCM options (!802)
- Fix Windows CRLF tokens in GEO reader and improve comment handling (!805)
- Use chrono in Timer (!807)
- Fix caching of FUNCTION tags that read from file and provide the same
functionality in FUNCTIONs defined for forcings (!759)
**NekMesh**:
- Add feature to read basic 2D geo files as CAD (!731)
......
......@@ -147,7 +147,7 @@ public:
/// Returns the output field
FIELD_UTILS_EXPORT LibUtilities::PtsFieldSharedPtr GetOutField() const;
/// Print statics of the interpolation weights
/// Returns if the weights have already been computed
FIELD_UTILS_EXPORT void PrintStatistics();
/// sets a callback funtion which gets called every time the interpolation
......
SET(SOLVER_UTILS_SOURCES
Core/Misc.cpp
Core/SessionFunction.cpp
AdvectionSystem.cpp
Advection/Advection.cpp
Advection/Advection3DHomogeneous1D.cpp
......@@ -44,6 +45,7 @@ SET(SOLVER_UTILS_SOURCES
SET(SOLVER_UTILS_HEADERS
Core/Misc.h
Core/SessionFunction.h
AdvectionSystem.h
Advection/Advection.h
Advection/AdvectionFR.h
......
////////////////////////////////////////////////////////////////////////////////
//
// File: SessionFunction.cpp
//
// For more information, please see: http://www.nektar.info/
//
// The MIT License
//
// Copyright (c) 2017 Kilian Lackhove
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Session Function
//
////////////////////////////////////////////////////////////////////////////////
#include <SolverUtils/Core/SessionFunction.h>
#include <LibUtilities/BasicUtils/VmathArray.hpp>
#include <boost/format.hpp>
#include <boost/function.hpp>
using namespace std;
namespace Nektar
{
namespace SolverUtils
{
/**
* Representation of a FUNCTION defined in the session xml file.
*
* @param session The session where the function was defined.
* @param field The field the function is defined on.
* @param functionName The name of the function.
* @param toCache Store the evaluated function for later use.
*/
SessionFunction::SessionFunction(const LibUtilities::SessionReaderSharedPtr &session,
const MultiRegions::ExpListSharedPtr &field,
std::string functionName,
bool toCache)
: m_session(session), m_field(field), m_name(functionName), m_toCache(toCache)
{
ASSERTL0(m_session->DefinesFunction(m_name),
"Function '" + m_name + "' does not exist.");
}
/**
* Evaluates a function defined in the xml session file at each quadrature point.
*
* @param pArray The array into which to write the values.
* @param pTime The time at which to evaluate the function.
* @param domain The domain to evaluate the function in.
*/
void SessionFunction::Evaluate(Array<OneD, Array<OneD, NekDouble> > &pArray,
const NekDouble pTime,
const int domain)
{
std::vector<std::string> vFieldNames = m_session->GetVariables();
for (int i = 0; i < vFieldNames.size(); i++)
{
Evaluate(vFieldNames[i], pArray[i], pTime, domain);
}
}
/**
* Evaluates a function defined in the xml session file at each quadrature point.
*
* @param pFieldNames The names of the fields to evaluate the function for.
* @param pArray The array into which to write the values.
* @param pTime The time at which to evaluate the function.
* @param domain The domain to evaluate the function in.
*/
void SessionFunction::Evaluate(std::vector<std::string> pFieldNames,
Array<OneD, Array<OneD, NekDouble> > &pArray,
const NekDouble &pTime,
const int domain)
{
ASSERTL1(pFieldNames.size() == pArray.num_elements(),
"Function '" + m_name +
"' variable list size mismatch with array storage.");
for (int i = 0; i < pFieldNames.size(); i++)
{
Evaluate(pFieldNames[i], pArray[i], pTime, domain);
}
}
/**
* Evaluates a function defined in the xml session file at each quadrature point.
*
* @param pFieldNames The names of the fields to evaluate the function for.
* @param pFields The fields into which to write the values.
* @param pTime The time at which to evaluate the function.
* @param domain The domain to evaluate the function in.
*/
void SessionFunction::Evaluate(
std::vector<std::string> pFieldNames,
Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
const NekDouble &pTime,
const int domain)
{
ASSERTL0(pFieldNames.size() == pFields.num_elements(),
"Field list / name list size mismatch.");
for (int i = 0; i < pFieldNames.size(); i++)
{
Evaluate(pFieldNames[i], pFields[i]->UpdatePhys(), pTime, domain);
pFields[i]->FwdTrans_IterPerExp(pFields[i]->GetPhys(),
pFields[i]->UpdateCoeffs());
}
}
/**
* Evaluates a function defined in the xml session file at each quadrature point.
*
* @param pFieldName The name of the field to evaluate the function for.
* @param pArray The array into which to write the values.
* @param pTime The time at which to evaluate the function.
* @param domain The domain to evaluate the function in.
*/
void SessionFunction::Evaluate(std::string pFieldName,
Array<OneD, NekDouble> &pArray,
const NekDouble &pTime,
const int domain)
{
LibUtilities::FunctionType vType =
m_session->GetFunctionType(m_name, pFieldName, domain);
unsigned int nq = m_field->GetNpoints();
std::pair<std::string, int> key(pFieldName, domain);
// sorry
if ((m_arrays.find(key) != m_arrays.end()) &&
(vType == LibUtilities::eFunctionTypeFile ||
((m_lastCached.find(key) != m_lastCached.end()) &&
(pTime - m_lastCached[key] < NekConstants::kNekZeroTol))))
{
// found cached field
if (pArray.num_elements() < nq)
{
pArray = Array<OneD, NekDouble>(nq);
}
Vmath::Vcopy(nq, m_arrays[key], 1, pArray, 1);
return;
}
if (vType == LibUtilities::eFunctionTypeExpression)
{
EvaluateExp(pFieldName, pArray, pTime, domain);
}
else if (vType == LibUtilities::eFunctionTypeFile ||
vType == LibUtilities::eFunctionTypeTransientFile)
{
std::string filename =
m_session->GetFunctionFilename(m_name, pFieldName, domain);
if (boost::filesystem::path(filename).extension() == ".pts")
{
EvaluatePts(pFieldName, pArray, pTime, domain);
}
else
{
EvaluateFld(pFieldName, pArray, pTime, domain);
}
}
if (m_toCache)
{
m_arrays[key] = Array<OneD, NekDouble>(nq);
Vmath::Vcopy(nq, pArray, 1, m_arrays[key], 1);
m_lastCached[key] = pTime;
}
}
/**
* @brief Provide a description of a function for a given field name.
*
* @param pFieldName Field name.
* @param domain The domain to evaluate the function in.
*/
std::string SessionFunction::Describe(std::string pFieldName, const int domain)
{
std::string retVal;
LibUtilities::FunctionType vType =
m_session->GetFunctionType(m_name, pFieldName, domain);
if (vType == LibUtilities::eFunctionTypeExpression)
{
LibUtilities::EquationSharedPtr ffunc =
m_session->GetFunction(m_name, pFieldName, domain);
retVal = ffunc->GetExpression();
}
else if (vType == LibUtilities::eFunctionTypeFile ||
LibUtilities::eFunctionTypeTransientFile)
{
std::string filename =
m_session->GetFunctionFilename(m_name, pFieldName, domain);
retVal = "from file " + filename;
}
return retVal;
}
/**
* Evaluates a function from expression
*
* @param pFieldName The name of the field to evaluate the function for.
* @param pArray The array into which to write the values.
* @param pTime The time at which to evaluate the function.
* @param domain The domain to evaluate the function in.
*/
void SessionFunction::EvaluateExp(string pFieldName,
Array<OneD, NekDouble> &pArray,
const NekDouble &pTime,
const int domain)
{
unsigned int nq = m_field->GetNpoints();
if (pArray.num_elements() < nq)
{
pArray = Array<OneD, NekDouble>(nq);
}
Array<OneD, NekDouble> x0(nq);
Array<OneD, NekDouble> x1(nq);
Array<OneD, NekDouble> x2(nq);
// Get the coordinates (assuming all fields have the same
// discretisation)
m_field->GetCoords(x0, x1, x2);
LibUtilities::EquationSharedPtr ffunc =
m_session->GetFunction(m_name, pFieldName, domain);
ffunc->Evaluate(x0, x1, x2, pTime, pArray);
}
/**
* Evaluates a function from fld file
*
* @param pFieldName The name of the field to evaluate the function for.
* @param pArray The array into which to write the values.
* @param pTime The time at which to evaluate the function.
* @param domain The domain to evaluate the function in.
*/
void SessionFunction::EvaluateFld(string pFieldName,
Array<OneD, NekDouble> &pArray,
const NekDouble &pTime,
const int domain)
{
unsigned int nq = m_field->GetNpoints();
if (pArray.num_elements() < nq)
{
pArray = Array<OneD, NekDouble>(nq);
}
std::string filename =
m_session->GetFunctionFilename(m_name, pFieldName, domain);
std::string fileVar =
m_session->GetFunctionFilenameVariable(m_name, pFieldName, domain);
if (fileVar.length() == 0)
{
fileVar = pFieldName;
}
// In case of eFunctionTypeTransientFile, generate filename from
// format string
LibUtilities::FunctionType vType =
m_session->GetFunctionType(m_name, pFieldName, domain);
if (vType == LibUtilities::eFunctionTypeTransientFile)
{
try
{
#if (defined _WIN32 && _MSC_VER < 1900)
// We need this to make sure boost::format has always
// two digits in the exponents of Scientific notation.
unsigned int old_exponent_format;
old_exponent_format = _set_output_format(_TWO_DIGIT_EXPONENT);
filename = boost::str(boost::format(filename) % pTime);
_set_output_format(old_exponent_format);
#else
filename = boost::str(boost::format(filename) % pTime);
#endif
}
catch (...)
{
ASSERTL0(false,
"Invalid Filename in function \"" + m_name +
"\", variable \"" + fileVar + "\"")
}
}
// Define list of global element ids
int numexp = m_field->GetExpSize();
Array<OneD, int> ElementGIDs(numexp);
for (int i = 0; i < numexp; ++i)
{
ElementGIDs[i] = m_field->GetExp(i)->GetGeom()->GetGlobalID();
}
std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
std::vector<std::vector<NekDouble> > FieldData;
LibUtilities::FieldIOSharedPtr fldIO =
LibUtilities::FieldIO::CreateForFile(m_session, filename);
fldIO->Import(filename,
FieldDef,
FieldData,
LibUtilities::NullFieldMetaDataMap,
ElementGIDs);
int idx = -1;
Array<OneD, NekDouble> vCoeffs(m_field->GetNcoeffs(), 0.0);
// Loop over all the expansions
for (int i = 0; i < FieldDef.size(); ++i)
{
// Find the index of the required field in the
// expansion segment
for (int j = 0; j < FieldDef[i]->m_fields.size(); ++j)
{
if (FieldDef[i]->m_fields[j] == fileVar)
{
idx = j;
}
}
if (idx >= 0)
{
m_field->ExtractDataToCoeffs(
FieldDef[i], FieldData[i], FieldDef[i]->m_fields[idx], vCoeffs);
}
else
{
cout << "Field " + fileVar + " not found." << endl;
}
}
m_field->BwdTrans_IterPerExp(vCoeffs, pArray);
}
/**
* Evaluates a function from pts file
*
* @param pFieldName The name of the field to evaluate the function for.
* @param pArray The array into which to write the values.
* @param pTime The time at which to evaluate the function.
* @param domain The domain to evaluate the function in.
*/
void SessionFunction::EvaluatePts(string pFieldName,
Array<OneD, NekDouble> &pArray,
const NekDouble &pTime,
const int domain)
{
unsigned int nq = m_field->GetNpoints();
if (pArray.num_elements() < nq)
{
pArray = Array<OneD, NekDouble>(nq);
}
std::string filename =
m_session->GetFunctionFilename(m_name, pFieldName, domain);
std::string fileVar =
m_session->GetFunctionFilenameVariable(m_name, pFieldName, domain);
if (fileVar.length() == 0)
{
fileVar = pFieldName;
}
// In case of eFunctionTypeTransientFile, generate filename from
// format string
LibUtilities::FunctionType vType =
m_session->GetFunctionType(m_name, pFieldName, domain);
if (vType == LibUtilities::eFunctionTypeTransientFile)
{
try
{
#if (defined _WIN32 && _MSC_VER < 1900)
// We need this to make sure boost::format has always
// two digits in the exponents of Scientific notation.
unsigned int old_exponent_format;
old_exponent_format = _set_output_format(_TWO_DIGIT_EXPONENT);
filename = boost::str(boost::format(filename) % pTime);
_set_output_format(old_exponent_format);
#else
filename = boost::str(boost::format(filename) % pTime);
#endif
}
catch (...)
{
ASSERTL0(false,
"Invalid Filename in function \"" + m_name +
"\", variable \"" + fileVar + "\"")
}
}
LibUtilities::PtsFieldSharedPtr outPts;
// check if we already loaded this file. For transient files,
// funcFilename != filename so we can make sure we only keep the
// latest pts field per funcFilename.
std::string funcFilename =
m_session->GetFunctionFilename(m_name, pFieldName, domain);
LibUtilities::PtsFieldSharedPtr inPts;
LibUtilities::PtsIO ptsIO(m_session->GetComm());
ptsIO.Import(filename, inPts);
Array<OneD, Array<OneD, NekDouble> > pts(inPts->GetDim() +
inPts->GetNFields());
for (int i = 0; i < inPts->GetDim() + inPts->GetNFields(); ++i)
{
pts[i] = Array<OneD, NekDouble>(nq);
}
if (inPts->GetDim() == 1)
{
m_field->GetCoords(pts[0]);
}
else if (inPts->GetDim() == 2)
{
m_field->GetCoords(pts[0], pts[1]);
}
else if (inPts->GetDim() == 3)
{
m_field->GetCoords(pts[0], pts[1], pts[2]);
}
outPts = MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(
inPts->GetDim(), inPts->GetFieldNames(), pts);
FieldUtils::Interpolator interp;
if (m_interpolators.find(funcFilename) != m_interpolators.end())
{
interp = m_interpolators[funcFilename];
}
else
{
interp = FieldUtils::Interpolator(Nektar::FieldUtils::eShepard);
if (m_session->GetComm()->GetRank() == 0)
{
interp.SetProgressCallback(&SessionFunction::PrintProgressbar,
this);
}
interp.CalcWeights(inPts, outPts);
if (m_session->GetComm()->GetRank() == 0)
{
cout << endl;
if (m_session->DefinesCmdLineArgument("verbose"))
{
interp.PrintStatistics();
}
}
}
if (m_toCache)
{
m_interpolators[funcFilename] = interp;
}
// TODO: only interpolate the field we actually want
interp.Interpolate(inPts, outPts);
int fieldInd;
vector<string> fieldNames = outPts->GetFieldNames();
for (fieldInd = 0; fieldInd < fieldNames.size(); ++fieldInd)
{
if (outPts->GetFieldName(fieldInd) == pFieldName)
{
break;
}
}
ASSERTL0(fieldInd != fieldNames.size(), "field not found");
pArray = outPts->GetPts(fieldInd + outPts->GetDim());
}
// end of namespaces
}
}
///////////////////////////////////////////////////////////////////////////////
//
// File SessionFunction.h
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2017 Kilian Lackhove
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Session Function
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_SOLVERUTILS_SESSIONFUNCTION_H
#define NEKTAR_SOLVERUTILS_SESSIONFUNCTION_H
#include <FieldUtils/Interpolator.h>
#include <LibUtilities/BasicUtils/FieldIO.h>
#include <LibUtilities/BasicUtils/FileSystem.h>
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/Progressbar.hpp>
#include <LibUtilities/BasicUtils/PtsField.h>
#include <LibUtilities/BasicUtils/PtsIO.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <MultiRegions/ExpList.h>
#include <SolverUtils/SolverUtilsDeclspec.h>
namespace Nektar
{
namespace SolverUtils
{
class SessionFunction
{
public:
/// Representation of a FUNCTION defined in the session xml file.
SOLVER_UTILS_EXPORT SessionFunction(
const LibUtilities::SessionReaderSharedPtr &session,
const MultiRegions::ExpListSharedPtr &field,
std::string functionName,
bool toCache = false);
/// Evaluates a function defined in the xml session file at each quadrature point.
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 Evaluate(
std::vector<std::string> pFieldNames,
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 Evaluate(
std::vector<std::string> pFieldNames,
Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,