Commit 0a82e8ac authored by Michael Turner's avatar Michael Turner
Browse files

Merge branch 'feature/SpongeBC' into 'master'

feature: Sponge Layers

See merge request !769
parents e861f516 79fc5473
......@@ -21,6 +21,8 @@ v5.0.0
- Remove unused files from BasicUtils (!841)
- Remove checks for old boost versions which are no longer supported (!841)
- Refactor ParseUtils to be more consistent (!843)
- Added support for using the distance to a specific region (e.g. outlet) in the
function definitions for the Absorption Forcing (!769)
**NekMesh**:
- Add feature to read basic 2D geo files as CAD (!731)
......
......@@ -24,9 +24,17 @@ This force type allows the user to apply an absorption layer (essentially a poro
<COEFF> [FUNCTION NAME] <COEFF/>
<REFFLOW> [FUNCTION NAME] <REFFLOW/>
<REFFLOWTIME> [FUNCTION NAME] <REFFLOWTIME/>
<BOUNDARYREGIONS> 1,4 <BOUNDARYREGIONS/>
</FORCE>
\end{lstlisting}
If a list of \inltt{BOUNDARYREGIONS} is specified, the distance to these regions is available as additional variable \inltt{r} in the definition of the \inltt{COEFF} function:
\begin{lstlisting}[style=XMLStyle]
<FUNCTION NAME="AbsorptionCoefficient">
<E VAR="p" EVARS="r" VALUE="-5000 * exp(-0.5 * (3*r / 0.4)^2)" />
<E VAR="u" EVARS="r" VALUE="-5000 * exp(-0.5 * (3*r / 0.4)^2)" />
<E VAR="v" EVARS="r" VALUE="-5000 * exp(-0.5 * (3*r / 0.4)^2)" />
</FUNCTION>
\end{lstlisting}
\subsection{Body}
This force type specifies the name of a body forcing function expressed in the \inltt{CONDITIONS} section.
......
......@@ -56,24 +56,34 @@ namespace Nektar
public:
LIB_UTILITIES_EXPORT Equation(const Equation &src):
m_vlist (src.m_vlist),
m_expr (src.m_expr),
m_expr_id (src.m_expr_id),
m_evaluator (src.m_evaluator)
{
}
LIB_UTILITIES_EXPORT Equation(const SessionReaderSharedPtr& session, const std::string& expr = ""):
LIB_UTILITIES_EXPORT Equation(const SessionReaderSharedPtr& session,
const std::string& expr = "",
const std::string& vlist = ""):
m_vlist (vlist),
m_expr (expr),
m_expr_id (-1),
m_evaluator (session->GetExpressionEvaluator())
{
boost::algorithm::trim(m_expr);
boost::algorithm::trim(m_vlist);
if (m_vlist.empty())
{
m_vlist = "x y z t";
}
try
{
if (!m_expr.empty())
{
m_expr_id = m_evaluator.DefineFunction("x y z t", m_expr);
m_expr_id = m_evaluator.DefineFunction(m_vlist, m_expr);
}
}
catch (const std::runtime_error& e)
......@@ -171,12 +181,24 @@ namespace Nektar
const Array<OneD, const NekDouble>& z,
const Array<OneD, const NekDouble>& t,
Array<OneD, NekDouble>& result) const
{
std::vector<Array<OneD, const NekDouble> > points;
points.push_back(x);
points.push_back(y);
points.push_back(z);
points.push_back(t);
Evaluate(points, result);
}
LIB_UTILITIES_EXPORT void Evaluate(
const std::vector<Array<OneD, const NekDouble> > points,
Array<OneD, NekDouble>& result) const
{
try
{
if (m_expr_id != -1)
{
m_evaluator.Evaluate(m_expr_id, x,y,z,t, result);
m_evaluator.Evaluate(m_expr_id, points, result);
}
}
catch (const std::runtime_error& e)
......@@ -208,6 +230,11 @@ namespace Nektar
return m_expr;
}
LIB_UTILITIES_EXPORT std::string GetVlist(void) const
{
return m_vlist;
}
/// Returns time spend on expression evaluation at
/// points (it does not include parse/pre-processing time).
LIB_UTILITIES_EXPORT NekDouble GetTime() const
......@@ -216,6 +243,7 @@ namespace Nektar
}
private:
std::string m_vlist;
std::string m_expr;
int m_expr_id;
AnalyticExpressionEvaluator& m_evaluator;
......
......@@ -2380,6 +2380,13 @@ namespace Nektar
domainStr = variable->Attribute("DOMAIN");
}
// If no domain string put to 0
std::string evarsStr = "x y z t";
if (variable->Attribute("EVARS"))
{
evarsStr = evarsStr + std::string(" ") + variable->Attribute("EVARS");
}
// Parse list of variables
std::vector<std::string> varSplit;
std::vector<unsigned int> domainList;
......@@ -2405,7 +2412,7 @@ namespace Nektar
// set expression
funcDef.m_expression = MemoryManager<Equation>
::AllocateSharedPtr(GetSharedThisPtr(),fcnStr);
::AllocateSharedPtr(GetSharedThisPtr(), fcnStr, evarsStr);
}
// Files are denoted by F
......
......@@ -6,6 +6,7 @@
//
// The MIT License
//
// Copyright (c) 2016 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).
......@@ -35,6 +36,9 @@
#include <SolverUtils/Forcing/ForcingAbsorption.h>
#include <LibUtilities/BasicUtils/Equation.h>
#include <LibUtilities/BasicUtils/ParseUtils.h>
using namespace std;
namespace Nektar
......@@ -60,41 +64,26 @@ namespace SolverUtils
const TiXmlElement* pForce)
{
m_NumVariable = pNumForcingFields;
int npts = pFields[0]->GetTotPoints();
const TiXmlElement* funcNameElmt;
funcNameElmt = pForce->FirstChildElement("COEFF");
ASSERTL0(funcNameElmt, "Requires COEFF tag, specifying function "
"name which prescribes absorption layer coefficient.");
int npts = pFields[0]->GetTotPoints();
string funcName = funcNameElmt->GetText();
ASSERTL0(m_session->DefinesFunction(funcName),
"Function '" + funcName + "' not defined.");
CalcAbsorption(pFields, pForce);
std::string s_FieldStr;
m_Absorption = Array<OneD, Array<OneD, NekDouble> > (m_NumVariable);
m_Forcing = Array<OneD, Array<OneD, NekDouble> > (m_NumVariable);
m_Forcing = Array<OneD, Array<OneD, NekDouble> >(m_NumVariable);
for (int i = 0; i < m_NumVariable; ++i)
{
s_FieldStr = m_session->GetVariable(i);
ASSERTL0(m_session->DefinesFunction(funcName, s_FieldStr),
"Variable '" + s_FieldStr + "' not defined.");
m_Absorption[i] = Array<OneD, NekDouble> (npts, 0.0);
m_Forcing[i] = Array<OneD, NekDouble> (npts, 0.0);
GetFunction(pFields, m_session, funcName)->Evaluate(s_FieldStr, m_Absorption[i]);
m_Forcing[i] = Array<OneD, NekDouble>(npts, 0.0);
}
funcNameElmt = pForce->FirstChildElement("REFFLOW");
const TiXmlElement* funcNameElmt = pForce->FirstChildElement("REFFLOW");
if (funcNameElmt)
{
string funcName = funcNameElmt->GetText();
ASSERTL0(m_session->DefinesFunction(funcName),
"Function '" + funcName + "' not defined.");
m_Refflow = Array<OneD, Array<OneD, NekDouble> > (m_NumVariable);
for (int i = 0; i < m_NumVariable; ++i)
{
s_FieldStr = m_session->GetVariable(i);
std::string s_FieldStr = m_session->GetVariable(i);
ASSERTL0(m_session->DefinesFunction(funcName, s_FieldStr),
"Variable '" + s_FieldStr + "' not defined.");
m_Refflow[i] = Array<OneD, NekDouble> (npts, 0.0);
......@@ -102,14 +91,138 @@ namespace SolverUtils
}
m_hasRefFlow = true;
}
funcNameElmt = pForce->FirstChildElement("REFFLOWTIME");
if (funcNameElmt)
{
m_funcNameTime = funcNameElmt->GetText();
m_hasRefFlowTime = true;
}
}
void ForcingAbsorption::CalcAbsorption(
const Nektar::Array<Nektar::OneD, Nektar::MultiRegions::ExpListSharedPtr>
&pFields,
const TiXmlElement *pForce)
{
const TiXmlElement *funcNameElmt = pForce->FirstChildElement("COEFF");
ASSERTL0(funcNameElmt,
"Requires COEFF tag, specifying function "
"name which prescribes absorption layer coefficient.");
string funcName = funcNameElmt->GetText();
ASSERTL0(m_session->DefinesFunction(funcName),
"Function '" + funcName + "' not defined.");
int npts = pFields[0]->GetTotPoints();
m_Absorption = Array<OneD, Array<OneD, NekDouble> >(m_NumVariable);
for (int i = 0; i < m_NumVariable; ++i)
{
m_Absorption[i] = Array<OneD, NekDouble>(npts, 0.0);
}
funcNameElmt = pForce->FirstChildElement("BOUNDARYREGIONS");
if (funcNameElmt)
{
ASSERTL0(ParseUtils::GenerateVector(funcNameElmt->GetText(),
m_bRegions),
"Unable to process list of BOUNDARYREGIONS in Absorption "
"Forcing: " +
std::string(funcNameElmt->GetText()));
// alter m_bRegions so that it contains the boundaryRegions of this rank
std::vector<unsigned int> localBRegions;
SpatialDomains::BoundaryConditions bcs(m_session, pFields[0]->GetGraph());
SpatialDomains::BoundaryRegionCollection regions = bcs.GetBoundaryRegions();
SpatialDomains::BoundaryRegionCollection::iterator it1;
int n = 0;
for (it1 = regions.begin(); it1 != regions.end(); ++it1)
{
if (std::find(m_bRegions.begin(), m_bRegions.end(), it1->first) != m_bRegions.end())
{
localBRegions.push_back(n);
}
n++;
}
m_bRegions = localBRegions;
if (m_bRegions.size() == 0)
{
return;
}
std::vector<Array<OneD, const NekDouble> > points;
Array<OneD, Array<OneD, NekDouble> > x(3);
for (int i = 0; i < 3; i++)
{
x[i] = Array<OneD, NekDouble>(npts, 0.0);
}
pFields[0]->GetCoords(x[0], x[1], x[2]);
for (int i = 0; i < 3; i++)
{
points.push_back(x[i]);
}
Array<OneD, NekDouble> t(npts, 0.0);
points.push_back(t);
Array<OneD, NekDouble> r(npts, 0.0);
std::vector<unsigned int>::iterator it;
std::vector<BPointPair> inPoints;
Array<OneD, Array<OneD, NekDouble> > b(3);
for (it = m_bRegions.begin(); it != m_bRegions.end(); ++it)
{
int bpts = pFields[0]->GetBndCondExpansions()[*it]->GetNpoints();
for (int i = 0; i < 3; i++)
{
b[i] = Array<OneD, NekDouble>(bpts, 0.0);
}
pFields[0]->GetBndCondExpansions()[*it]->GetCoords(
b[0], b[1], b[2]);
for (int i = 0;
i < pFields[0]->GetBndCondExpansions()[*it]->GetNpoints();
++i)
{
inPoints.push_back(
BPointPair(BPoint(b[0][i], b[1][i], b[2][i]), i));
}
}
m_rtree = MemoryManager<BRTree>::AllocateSharedPtr();
m_rtree->insert(inPoints.begin(), inPoints.end());
for (int i = 0; i < npts; ++i)
{
std::vector<BPointPair> result;
BPoint sPoint(x[0][i], x[1][i], x[2][i]);
m_rtree->query(bgi::nearest(sPoint, 1), std::back_inserter(result));
r[i] = bg::distance(sPoint, result[0].first);
}
points.push_back(r);
std::string s_FieldStr;
for (int i = 0; i < m_NumVariable; ++i)
{
s_FieldStr = m_session->GetVariable(i);
ASSERTL0(m_session->DefinesFunction(funcName, s_FieldStr),
"Variable '" + s_FieldStr + "' not defined.");
LibUtilities::EquationSharedPtr ffunc =
m_session->GetFunction(funcName, s_FieldStr);
ASSERTL0(ffunc->GetVlist() == "x y z t r",
"EVARS of " + funcName + " must be 'r'");
ffunc->Evaluate(points, m_Absorption[i]);
}
}
else
{
for (int i = 0; i < m_NumVariable; ++i)
{
std::string s_FieldStr = m_session->GetVariable(i);
GetFunction(pFields, m_session, funcName)->Evaluate(s_FieldStr, m_Absorption[i]);
}
}
}
void ForcingAbsorption::v_Apply(
......
......@@ -38,12 +38,19 @@
#include <string>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/index/rtree.hpp>
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <MultiRegions/ExpList.h>
#include <SolverUtils/SolverUtilsDeclspec.h>
#include <SolverUtils/Forcing/Forcing.h>
namespace bg = boost::geometry;
namespace bgi = boost::geometry::index;
namespace Nektar
{
namespace SolverUtils
......@@ -70,11 +77,18 @@ namespace SolverUtils
static std::string className;
protected:
typedef bg::model::point<NekDouble, 3, bg::cs::cartesian> BPoint;
typedef std::pair<BPoint, unsigned int> BPointPair;
typedef bgi::rtree<BPointPair, bgi::rstar<16> > BRTree;
bool m_hasRefFlow;
bool m_hasRefFlowTime;
Array<OneD, Array<OneD, NekDouble> > m_Absorption;
Array<OneD, Array<OneD, NekDouble> > m_Refflow;
std::string m_funcNameTime;
std::vector<unsigned int> m_bRegions;
std::shared_ptr<BRTree> m_rtree;
SOLVER_UTILS_EXPORT virtual void v_InitObject(
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
......@@ -91,8 +105,10 @@ namespace SolverUtils
ForcingAbsorption(const LibUtilities::SessionReaderSharedPtr& pSession);
virtual ~ForcingAbsorption(void){};
void CalcAbsorption(
const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
const TiXmlElement *pForce);
};
}
}
// Hui XU 21 Jul 2013 Created
......
......@@ -370,7 +370,7 @@ void APE::DoOdeRhs(const Array<OneD, const Array<OneD, NekDouble> >&inarray,
for (auto &x : m_forcing)
{
x->Apply(m_fields, outarray, outarray, m_time);
x->Apply(m_fields, inarray, outarray, m_time);
}
}
......
......@@ -8,14 +8,14 @@
</files>
<metrics>
<metric type="L2" id="1">
<value variable="p" tolerance="1e-12">7.93225</value>
<value variable="u" tolerance="1e-12">0.0189841</value>
<value variable="v" tolerance="1e-12">0.0109285</value>
<value variable="p" tolerance="1e-4">6.7577</value>
<value variable="u" tolerance="1e-7">0.014607</value>
<value variable="v" tolerance="1e-7">0.00571716</value>
</metric>
<metric type="Linf" id="2">
<value variable="p" tolerance="1e-12">18.7643</value>
<value variable="u" tolerance="1e-12">0.0322549</value>
<value variable="v" tolerance="1e-12">0.0460461</value>
<value variable="p" tolerance="1e-4">13.659</value>
<value variable="u" tolerance="1e-7">0.0280936</value>
<value variable="v" tolerance="1e-7">0.0108992</value>
</metric>
</metrics>
</test>
......@@ -337,9 +337,9 @@
</BOUNDARYREGIONS>
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="p" VALUE="0"/>
<D VAR="u" VALUE="0"/>
<D VAR="v" VALUE="0"/>
<D VAR="p" USERDEFINEDTYPE="RiemannInvariantBC" VALUE="0" />
<D VAR="u" USERDEFINEDTYPE="RiemannInvariantBC" VALUE="0" />
<D VAR="v" USERDEFINEDTYPE="RiemannInvariantBC" VALUE="0" />
</REGION>
<REGION REF="1">
<D VAR="p" USERDEFINEDTYPE="Wall" VALUE="0" />
......@@ -358,10 +358,27 @@
<E VAR="u" VALUE="0"/>
<E VAR="v" VALUE="0"/>
</FUNCTION>
<FUNCTION NAME="ExactSolution"> <!-- Not really the exact solution -->
<E VAR="p" VALUE="0"/>
<E VAR="u" VALUE="0"/>
<E VAR="v" VALUE="0"/>
</FUNCTION>
<FUNCTION NAME="InitialConditions">
<E VAR="p" VALUE="100*exp(-32*(x^2+y^2))"/> <!-- Gaussian pulse located at the origin -->
<E VAR="u" VALUE="0"/>
<E VAR="v" VALUE="0"/>
</FUNCTION>
<!-- 0.4 m thick sponge layer around the RiemannInvariantBC BCs -->
<FUNCTION NAME="SpongeCoefficient">
<E VAR="p" EVARS="r" VALUE="-5000 * exp(-0.5 * (3*r / 0.4)^2)" />
<E VAR="u" EVARS="r" VALUE="-5000 * exp(-0.5 * (3*r / 0.4)^2)" />
<E VAR="v" EVARS="r" VALUE="-5000 * exp(-0.5 * (3*r / 0.4)^2)" />
</FUNCTION>
</CONDITIONS>
<FORCING>
<FORCE TYPE="Absorption">
<COEFF> SpongeCoefficient </COEFF>
<BOUNDARYREGIONS> 0 </BOUNDARYREGIONS>
</FORCE>
</FORCING>
</NEKTAR>
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment