Commit dfb590bc authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Added forcing section to XML input.

Created forcing tag to instantiate forcing objects.
Functions still defined in conditions for now.
parent 7112ac4a
......@@ -544,7 +544,7 @@ namespace Nektar
m_session->LoadParameter("NumQuadPointsError",
m_NumQuadPointsError, 0);
if (m_session->DefinesFunction("BodyForce"))
/* if (m_session->DefinesFunction("BodyForce"))
{
m_forces = Array<OneD, MultiRegions::ExpListSharedPtr>(v_GetForceDimension());
int nq = m_fields[0]->GetNpoints();
......@@ -632,7 +632,7 @@ namespace Nektar
}
}
}
*/
// If a tangent vector policy is defined then the local tangent
// vectors on each element need to be generated
if (m_session->DefinesGeometricInfo("TANGENTDIR"))
......@@ -642,15 +642,9 @@ namespace Nektar
// Zero all physical fields initially
ZeroPhysFields();
// Forcing term
if(m_session->DefinesFunction("SpongeCoefficient"))
{
m_SpongeForcing = SolverUtils::GetForcingFactory().CreateInstance("SpongeForcing", m_session, m_fields);
}
if (m_session->DefinesFunction("BodyForce"))
{
m_BodyForcing = SolverUtils::GetForcingFactory().CreateInstance("BodyForcing", m_session, m_fields);
}
// Forcing terms
m_forcing = SolverUtils::Forcing::Load(m_session, m_fields);
}
/**
......
......@@ -406,10 +406,8 @@ namespace Nektar
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions;
/// Pointer to graph defining mesh.
SpatialDomains::MeshGraphSharedPtr m_graph;
/// Sponge forcing term
SolverUtils::ForcingSharedPtr m_SpongeForcing;
/// Sponge forcing term
SolverUtils::ForcingSharedPtr m_BodyForcing;
/// Forcing terms
std::vector<ForcingSharedPtr> m_forcing;
/// Filename.
std::string m_filename;
/// Name of the session.
......
......@@ -54,9 +54,10 @@ namespace Nektar
}
void Forcing::InitObject(
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields)
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const TiXmlElement* pForce)
{
v_InitObject(pFields);
v_InitObject(pFields, pForce);
}
void Forcing::Apply(
......@@ -67,6 +68,39 @@ namespace Nektar
v_Apply(fields, inarray, outarray);
}
/**
*
*/
vector<ForcingSharedPtr> Forcing::Load(
const LibUtilities::SessionReaderSharedPtr& pSession,
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields)
{
vector<ForcingSharedPtr> vForceList;
if (!pSession->DefinesElement("Nektar/Forcing"))
{
return vForceList;
}
TiXmlElement* vForcing = pSession->GetElement("Nektar/Forcing");
if (vForcing)
{
TiXmlElement* vForce = vForcing->FirstChildElement("FORCE");
while (vForce)
{
string vType = vForce->Attribute("TYPE");
vForceList.push_back(GetForcingFactory().CreateInstance(
vType, pSession, pFields, vForce));
vForce = vForce->NextSiblingElement("FORCE");
}
}
return vForceList;
}
void Forcing::EvaluateFunction(
Array<OneD, MultiRegions::ExpListSharedPtr> pFields,
LibUtilities::SessionReaderSharedPtr pSession,
......
......@@ -47,6 +47,21 @@ namespace Nektar
{
namespace SolverUtils
{
// Forward declaration
class Forcing;
/// A shared pointer to an EquationSystem object
SOLVER_UTILS_EXPORT typedef boost::shared_ptr<Forcing> ForcingSharedPtr;
/// Declaration of the forcing factory
typedef LibUtilities::NekFactory<std::string, Forcing,
const LibUtilities::SessionReaderSharedPtr&,
const Array<OneD, MultiRegions::ExpListSharedPtr>&,
const TiXmlElement*> ForcingFactory;
/// Declaration of the forcing factory singleton
SOLVER_UTILS_EXPORT ForcingFactory& GetForcingFactory();
/**
* @class Forcing
* @brief Defines a forcing term to be explicitly applied.
......@@ -56,7 +71,8 @@ namespace SolverUtils
public:
/// Initialise the forcing object
SOLVER_UTILS_EXPORT void InitObject(
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields);
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const TiXmlElement* pForce);
/// Apply the forcing
SOLVER_UTILS_EXPORT void Apply(
......@@ -64,6 +80,10 @@ namespace SolverUtils
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
static std::vector<ForcingSharedPtr> Load(
const LibUtilities::SessionReaderSharedPtr& pSession,
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields);
protected:
/// Session reader
LibUtilities::SessionReaderSharedPtr m_session;
......@@ -76,7 +96,8 @@ namespace SolverUtils
Forcing(const LibUtilities::SessionReaderSharedPtr&);
virtual void v_InitObject(
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields) = 0;
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const TiXmlElement* pForce) = 0;
virtual void v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
......@@ -91,18 +112,6 @@ namespace SolverUtils
NekDouble pTime = NekDouble(0));
};
/// A shared pointer to an EquationSystem object
SOLVER_UTILS_EXPORT typedef boost::shared_ptr<Forcing> ForcingSharedPtr;
/// Declaration of the forcing factory
typedef LibUtilities::NekFactory<std::string, Forcing,
const LibUtilities::SessionReaderSharedPtr&,
const Array<OneD, MultiRegions::ExpListSharedPtr>& > ForcingFactory;
/// Declaration of the forcing factory singleton
SOLVER_UTILS_EXPORT ForcingFactory& GetForcingFactory();
}
}
......
......@@ -41,7 +41,7 @@ namespace SolverUtils
{
std::string ForcingBody::className = GetForcingFactory().
RegisterCreatorFunction("BodyForcing",
RegisterCreatorFunction("Body",
ForcingBody::create,
"Body Forcing");
......@@ -52,28 +52,8 @@ namespace SolverUtils
}
void ForcingBody::v_InitObject(
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields)
{
ReadForceInfo(pFields);
}
void ForcingBody::v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
if (m_session->DefinesFunction("BodyForce"))
{
for (int i = 0; i < m_NumVariable; i++)
{
Vmath::Vadd(outarray[i].num_elements(), outarray[i], 1,
m_Forcing[i], 1, outarray[i], 1);
}
}
}
void ForcingBody::ReadForceInfo(
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields)
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const TiXmlElement* pForce)
{
std::string m_SolverInfo = m_session->GetSolverInfo("SolverType");
int nvariables = m_session->GetVariables().size();
......@@ -110,8 +90,23 @@ namespace SolverUtils
m_Forcing[i], "BodyForce");
}
}
}
void ForcingBody::v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
if (m_session->DefinesFunction("BodyForce"))
{
for (int i = 0; i < m_NumVariable; i++)
{
Vmath::Vadd(outarray[i].num_elements(), outarray[i], 1,
m_Forcing[i], 1, outarray[i], 1);
}
}
}
}
}
......@@ -57,11 +57,12 @@ namespace SolverUtils
/// Creates an instance of this class
SOLVER_UTILS_EXPORT static ForcingSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession,
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields)
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const TiXmlElement* pForce)
{
ForcingSharedPtr p = MemoryManager<ForcingBody>::
AllocateSharedPtr(pSession);
p->InitObject(pFields);
p->InitObject(pFields, pForce);
return p;
}
......@@ -70,7 +71,8 @@ namespace SolverUtils
protected:
virtual void v_InitObject(
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields);
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const TiXmlElement* pForce);
virtual void v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
......@@ -80,8 +82,6 @@ namespace SolverUtils
private:
ForcingBody(const LibUtilities::SessionReaderSharedPtr& pSession);
void ReadForceInfo(
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields);
};
}
......
......@@ -41,7 +41,7 @@ namespace SolverUtils
{
std::string ForcingSponge::className = GetForcingFactory().
RegisterCreatorFunction("SpongeForcing",
RegisterCreatorFunction("Sponge",
ForcingSponge::create,
"Forcing Sponge");
......@@ -51,42 +51,8 @@ namespace SolverUtils
}
void ForcingSponge::v_InitObject(
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields)
{
ReadSpongeInfo(pFields);
}
void ForcingSponge::v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
if (m_session->DefinesFunction("RefFields"))
{
for (int i = 0; i < m_NumVariable; i++)
{
Vmath::Vsub(m_Forcing[i].num_elements(), inarray[i], 1,
m_Refflow[i], 1, m_Forcing[i], 1);
Vmath::Vmul(m_Forcing[i].num_elements(), m_Sponge[i], 1,
m_Forcing[i], 1, m_Forcing[i], 1);
Vmath::Vadd(outarray[i].num_elements(), m_Forcing[i], 1,
outarray[i], 1, outarray[i], 1);
}
}
else
{
for (int i = 0; i < m_NumVariable; i++)
{
Vmath::Vmul(m_Forcing[i].num_elements(), m_Sponge[i], 1,
inarray[i], 1, m_Forcing[i], 1);
Vmath::Vadd(outarray[i].num_elements(), m_Forcing[i], 1,
outarray[i], 1, outarray[i], 1);
}
}
}
void ForcingSponge::ReadSpongeInfo(
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields)
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const TiXmlElement* pForce)
{
std::string m_SolverInfo = m_session->GetSolverInfo("SolverType");
int nvariables = m_session->GetVariables().size();
......@@ -100,14 +66,19 @@ namespace SolverUtils
{
m_NumVariable = nvariables; // e.g. (u v w) for 3D case
}
std::string s_FieldStr;
if (m_session->DefinesFunction("SpongeCoefficient"))
{
m_Sponge = 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);
m_Sponge[i] = Array<OneD, NekDouble> (npts, 0.0);
m_Forcing[i] = Array<OneD, NekDouble> (npts, 0.0);
EvaluateFunction(pFields, m_session, s_FieldStr,
m_Sponge[i], "SpongeCoefficient");
}
}
if (m_session->DefinesFunction("RefFields"))
......@@ -115,27 +86,42 @@ namespace SolverUtils
m_Refflow = Array<OneD, Array<OneD, NekDouble> > (m_NumVariable);
for (int i = 0; i < m_NumVariable; ++i)
{
s_FieldStr = m_session->GetVariable(i);
m_Refflow[i] = Array<OneD, NekDouble> (npts, 0.0);
EvaluateFunction(pFields, m_session, s_FieldStr,
m_Refflow[i], "RefFields");
}
}
}
std::string s_FieldStr;
for (int i = 0; i < m_NumVariable; ++i)
void ForcingSponge::v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
if (m_session->DefinesFunction("RefFields"))
{
s_FieldStr = m_session->GetVariable(i);
if (m_session->DefinesFunction("SpongeCoefficient"))
for (int i = 0; i < m_NumVariable; i++)
{
EvaluateFunction(pFields, m_session, s_FieldStr,
m_Sponge[i], "SpongeCoefficient");
Vmath::Vsub(m_Forcing[i].num_elements(), inarray[i], 1,
m_Refflow[i], 1, m_Forcing[i], 1);
Vmath::Vmul(m_Forcing[i].num_elements(), m_Sponge[i], 1,
m_Forcing[i], 1, m_Forcing[i], 1);
Vmath::Vadd(outarray[i].num_elements(), m_Forcing[i], 1,
outarray[i], 1, outarray[i], 1);
}
if (m_session->DefinesFunction("RefFields"))
}
else
{
for (int i = 0; i < m_NumVariable; i++)
{
EvaluateFunction(pFields, m_session, s_FieldStr,
m_Refflow[i], "RefFields");
Vmath::Vmul(m_Forcing[i].num_elements(), m_Sponge[i], 1,
inarray[i], 1, m_Forcing[i], 1);
Vmath::Vadd(outarray[i].num_elements(), m_Forcing[i], 1,
outarray[i], 1, outarray[i], 1);
}
}
}/// The end of reading sponge info.
}
}
}
/// Hui XU 2013 Jul 21
......@@ -56,11 +56,12 @@ namespace SolverUtils
/// Creates an instance of this class
SOLVER_UTILS_EXPORT static ForcingSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession,
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields)
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const TiXmlElement* pForce)
{
ForcingSharedPtr p = MemoryManager<ForcingSponge>::
AllocateSharedPtr(pSession);
p->InitObject(pFields);
p->InitObject(pFields, pForce);
return p;
}
......@@ -72,7 +73,8 @@ namespace SolverUtils
Array<OneD, Array<OneD, NekDouble> > m_Refflow;
virtual void v_InitObject(
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields);
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const TiXmlElement* pForce);
virtual void v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
......@@ -82,8 +84,6 @@ namespace SolverUtils
private:
ForcingSponge(const LibUtilities::SessionReaderSharedPtr& pSession);
void ReadSpongeInfo(
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields);
};
}
......
......@@ -457,10 +457,6 @@ namespace Nektar
timer.Start();
// evaluate convection terms
m_advObject->DoAdvection(m_fields, m_nConvectiveFields, m_velocity,inarray,outarray,m_time);
if(m_session->DefinesFunction("SpongeCoefficient"))
{
m_SpongeForcing->Apply(m_fields,inarray,outarray);
}
timer.Stop();
if(m_showTimings&&IsRoot)
{
......@@ -475,16 +471,17 @@ namespace Nektar
}
}
//add the force
if(m_session->DefinesFunction("BodyForce"))
// apply forcing
timer.Start();
std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
{
timer.Start();
m_BodyForcing->Apply(m_fields,inarray,outarray);
timer.Stop();
if(m_showTimings&&IsRoot)
{
cout << "\t Body ForceTime : "<< timer.TimePerTest(1) << endl;
}
(*x)->Apply(m_fields, inarray, outarray);
}
timer.Stop();
if(m_showTimings&&IsRoot)
{
cout << "\t Body ForceTime : "<< timer.TimePerTest(1) << endl;
}
if(m_pressureHBCs[0].num_elements() > 0)
......
Markdown is supported
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