Commit 65a01278 authored by Yumnah Mohamied's avatar Yumnah Mohamied
Browse files

high order outflow condition and sponge layer adjustments

parent a576c8b9
......@@ -53,7 +53,7 @@
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#ifndef MAX_PARAM
#define MAX_PARAM 4 // default maximum number of parameters to support
#define MAX_PARAM 5 // default maximum number of parameters to support
#endif
namespace Nektar
......
......@@ -1255,8 +1255,8 @@ namespace Nektar
// Find nearest element
if (!elmtIdDist.empty())
{
NekDouble min_d = elmtIdDist[0].second;
int min_id = elmtIdDist[0].first;
NekDouble min_d = elmtIdDist[0].first;
int min_id = elmtIdDist[0].second;
for (int i = 1; i < elmtIdDist.size(); ++i)
{
......@@ -1266,10 +1266,6 @@ namespace Nektar
}
}
// retrieve local coordinates of chosen point
(*m_exp)[min_id]->GetGeom()->ContainsPoint(gloCoords,
locCoords,
tol, resid);
return min_id;
}
else {
......@@ -1313,24 +1309,6 @@ namespace Nektar
}
/**
* The operation is evaluated locally by the elemental
* function StdRegions#StdExpansion#GetSurfaceNormal.
*/
void ExpList::GetSurfaceNormal(Array<OneD, NekDouble> &SurfaceNormal,
const int k)
{
int i;
Array<OneD, NekDouble> normals;
for(i = 0; i < (*m_exp).size(); ++i)
{
//e_SN = SurfaceNormal + m_phys_offset[i];
normals = (*m_exp)[i]->GetSurfaceNormal()[k];
Vmath::Vcopy(normals.num_elements(), &normals[0], 1, &SurfaceNormal[0] + m_phys_offset[i], 1);
}
}
/**
* Configures geometric info, such as tangent direction, on each
* expansion.
......
......@@ -347,12 +347,6 @@ namespace Nektar
Array<OneD, NekDouble> &outarray,
int BndID);
/// This function calculates Surface Normal vector of a smooth
/// manifold.
MULTI_REGIONS_EXPORT void GetSurfaceNormal(
Array<OneD,NekDouble> &SurfaceNormal,
const int k);
/// Apply geometry information to each expansion.
MULTI_REGIONS_EXPORT void ApplyGeomInfo();
......
......@@ -70,9 +70,10 @@ namespace Nektar
void Forcing::Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr>& fields,
const Array<OneD, Array<OneD, NekDouble> >& inarray,
Array<OneD, Array<OneD, NekDouble> >& outarray)
Array<OneD, Array<OneD, NekDouble> >& outarray,
const NekDouble& time)
{
v_Apply(fields, inarray, outarray);
v_Apply(fields, inarray, outarray, time);
}
......@@ -114,6 +115,26 @@ namespace Nektar
return vForceList;
}
void Forcing::EvaluateTimeFunction(
LibUtilities::SessionReaderSharedPtr pSession,
std::string pFieldName,
Array<OneD, NekDouble>& pArray,
const std::string& pFunctionName,
NekDouble pTime)
{
ASSERTL0(pSession->DefinesFunction(pFunctionName),
"Function '" + pFunctionName + "' does not exist.");
LibUtilities::EquationSharedPtr ffunc =
pSession->GetFunction(pFunctionName, pFieldName);
Array<OneD, NekDouble> x0(1,0.0);
Array<OneD, NekDouble> x1(1,0.0);
Array<OneD, NekDouble> x2(1,0.0);
ffunc->Evaluate(x0, x1, x2, pTime, pArray);
}
void Forcing::EvaluateFunction(
Array<OneD, MultiRegions::ExpListSharedPtr> pFields,
......@@ -139,11 +160,11 @@ namespace Nektar
Array<OneD, NekDouble> x0(nq);
Array<OneD, NekDouble> x1(nq);
Array<OneD, NekDouble> x2(nq);
pFields[0]->GetCoords(x0, x1, x2);
LibUtilities::EquationSharedPtr ffunc =
pSession->GetFunction(pFunctionName, pFieldName);
pSession->GetFunction(pFunctionName, pFieldName);
ffunc->Evaluate(x0, x1, x2, pTime, pArray);
}
else if (vType == LibUtilities::eFunctionTypeFile)
......
......@@ -80,7 +80,8 @@ namespace SolverUtils
SOLVER_UTILS_EXPORT void Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time);
SOLVER_UTILS_EXPORT static std::vector<ForcingSharedPtr> Load(
const LibUtilities::SessionReaderSharedPtr& pSession,
......@@ -107,15 +108,24 @@ namespace SolverUtils
SOLVER_UTILS_EXPORT virtual void v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)=0;
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)=0;
void EvaluateFunction(
Array<OneD, MultiRegions::ExpListSharedPtr> pFields,
LibUtilities::SessionReaderSharedPtr pSession,
std::string pFieldName, Array<OneD, NekDouble>& pArray,
LibUtilities::SessionReaderSharedPtr pSession,
std::string pFieldName,
Array<OneD, NekDouble>& pArray,
const std::string& pFunctionName,
NekDouble pTime = NekDouble(0));
void EvaluateTimeFunction(
LibUtilities::SessionReaderSharedPtr pSession,
std::string pFieldName,
Array<OneD, NekDouble>& pArray,
const std::string& pFunctionName,
NekDouble pTime = NekDouble(0));
};
}
}
......
......@@ -110,8 +110,9 @@ namespace SolverUtils
void ForcingBody::v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)
{
for (int i = 0; i < m_NumVariable; i++)
{
Vmath::Vadd(outarray[i].num_elements(), outarray[i], 1,
......
......@@ -79,7 +79,8 @@ namespace SolverUtils
SOLVER_UTILS_EXPORT virtual void v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time);
private:
ForcingBody(const LibUtilities::SessionReaderSharedPtr& pSession);
......
......@@ -80,7 +80,8 @@ namespace SolverUtils
void ForcingNoise::v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)
{
for (int i = 0; i < m_NumVariable; i++)
{
......
......@@ -79,7 +79,8 @@ namespace SolverUtils
SOLVER_UTILS_EXPORT virtual void v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time);
private:
ForcingNoise(const LibUtilities::SessionReaderSharedPtr& pSession);
......
......@@ -74,7 +74,8 @@ namespace SolverUtils
void ForcingProgrammatic::v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)
{
for (int i = 0; i < m_NumVariable; i++)
{
......
......@@ -88,7 +88,8 @@ namespace SolverUtils
SOLVER_UTILS_EXPORT virtual void v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time);
private:
ForcingProgrammatic(const LibUtilities::SessionReaderSharedPtr& pSession);
......
......@@ -101,20 +101,49 @@ namespace SolverUtils
}
m_hasRefFlow = true;
}
funcNameElmt = pForce->FirstChildElement("REFFLOWTIME");
if (funcNameElmt)
{
m_funcNameTime = funcNameElmt->GetText();
ASSERTL0(m_session->DefinesFunction(funcName),
"Function '" + funcName + "' not defined.");
m_hasRefFlowTime = true;
}
}
void ForcingSponge::v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)
{
int nq = m_Forcing[0].num_elements();
std::string s_FieldStr;
Array<OneD, NekDouble> TimeScale(1);
Array<OneD, Array<OneD, NekDouble> > RefflowScaled(m_NumVariable);
if (m_hasRefFlow)
{
for (int i = 0; i < m_NumVariable; i++)
{
RefflowScaled[i] = Array<OneD, NekDouble> (nq);
if (m_hasRefFlowTime)
{
s_FieldStr = m_session->GetVariable(i);
EvaluateTimeFunction(m_session, s_FieldStr, TimeScale, m_funcNameTime, time);
Vmath::Smul(nq, TimeScale[0], m_Refflow[i],1,RefflowScaled[i],1);
}
else
{
Vmath::Vcopy(nq, m_Refflow[i],1, RefflowScaled[i],1);
}
Vmath::Vsub(nq, inarray[i], 1,
m_Refflow[i], 1, m_Forcing[i], 1);
RefflowScaled[i], 1, m_Forcing[i], 1);
Vmath::Vmul(nq, m_Sponge[i], 1,
m_Forcing[i], 1, m_Forcing[i], 1);
Vmath::Vadd(nq, m_Forcing[i], 1,
......@@ -132,6 +161,6 @@ namespace SolverUtils
}
}
}
}
}
......@@ -71,9 +71,11 @@ namespace SolverUtils
protected:
bool m_hasRefFlow;
bool m_hasRefFlowTime;
Array<OneD, Array<OneD, NekDouble> > m_Sponge;
Array<OneD, Array<OneD, NekDouble> > m_Refflow;
std::string m_funcNameTime;
SOLVER_UTILS_EXPORT virtual void v_InitObject(
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const unsigned int& pNumForcingFields,
......@@ -82,7 +84,8 @@ namespace SolverUtils
SOLVER_UTILS_EXPORT virtual void v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time);
private:
ForcingSponge(const LibUtilities::SessionReaderSharedPtr& pSession);
......
......@@ -966,6 +966,11 @@ namespace Nektar
v_AddFaceNormBoundaryInt(face,FaceExp,Fn,outarray);
}
int StdExpansion::v_GetElmtId(void)
{
return m_elmt_id;
}
const Array<OneD, const NekDouble>& StdExpansion::v_GetPhysNormals(void)
{
NEKERROR(ErrorUtil::efatal, "This function is not valid for this class");
......@@ -1686,7 +1691,7 @@ namespace Nektar
return result;
}
const NormalVector & StdExpansion::v_GetSurfaceNormal() const
const NormalVector & StdExpansion::v_GetSurfaceNormal(const int id) const
{
ASSERTL0(false, "Cannot get face normals for this expansion.");
static NormalVector result;
......
......@@ -1128,6 +1128,10 @@ namespace Nektar
return v_GetMetricInfo();
}
/// \brief Get the element id of this expansion when used
/// in a list by returning value of #m_elmt_id
STD_REGIONS_EXPORT virtual int v_GetElmtId();
STD_REGIONS_EXPORT virtual const Array<OneD, const NekDouble>& v_GetPhysNormals(void);
STD_REGIONS_EXPORT virtual void v_SetPhysNormals(Array<OneD, const NekDouble> &normal);
......@@ -1251,10 +1255,9 @@ namespace Nektar
return v_GetVertexNormal(vertex);
}
const NormalVector & GetSurfaceNormal() const
const NormalVector & GetSurfaceNormal(const int id) const
{
// @TODO Implement this
return v_GetSurfaceNormal();
return v_GetSurfaceNormal(id);
}
const LibUtilities::PointsKeyVector GetPointsKeys() const
......@@ -1721,7 +1724,7 @@ namespace Nektar
STD_REGIONS_EXPORT virtual void v_ComputeVertexNormal(const int vertex);
STD_REGIONS_EXPORT virtual const NormalVector & v_GetFaceNormal(const int face) const;
STD_REGIONS_EXPORT virtual const NormalVector & v_GetSurfaceNormal() const;
STD_REGIONS_EXPORT virtual const NormalVector & v_GetSurfaceNormal(const int id) const;
STD_REGIONS_EXPORT virtual Array<OneD, unsigned int>
v_GetEdgeInverseBoundaryMap(int eid);
......
......@@ -307,9 +307,9 @@ namespace Nektar
}
const NormalVector & StdExpansion3D::v_GetSurfaceNormal() const
const NormalVector & StdExpansion3D::v_GetSurfaceNormal(const int id) const
{
return m_surfaceNormal;
return v_GetFaceNormal(id);
}
const NormalVector & StdExpansion3D::v_GetFaceNormal(const int face) const
......
......@@ -192,7 +192,6 @@ namespace Nektar
STD_REGIONS_EXPORT virtual void v_NegateFaceNormal(
const int face);
NormalVector m_surfaceNormal;
std::map<int, NormalVector> m_faceNormals;
std::map<int, bool> m_negatedNormals;
......@@ -208,7 +207,7 @@ namespace Nektar
{
return 3;
}
STD_REGIONS_EXPORT const NormalVector & v_GetSurfaceNormal() const;
STD_REGIONS_EXPORT const NormalVector & v_GetSurfaceNormal(const int id) const;
STD_REGIONS_EXPORT const NormalVector & v_GetFaceNormal(const int face) const;
};
} //end of namespace
......
......@@ -63,7 +63,7 @@ namespace Nektar
void CoupledLinearNS::v_InitObject()
{
EquationSystem::v_InitObject();
UnsteadySystem::v_InitObject();
IncNavierStokes::v_InitObject();
int i;
......@@ -1337,7 +1337,7 @@ namespace Nektar
std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
{
(*x)->Apply(m_fields, outarray, outarray);
(*x)->Apply(m_fields, outarray, outarray,time);
}
}
......@@ -1484,7 +1484,8 @@ namespace Nektar
std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
{
(*x)->Apply(m_fields, forcing_phys, forcing_phys);
const NekDouble time=0;
(*x)->Apply(m_fields, forcing_phys, forcing_phys, time);
}
for (unsigned int i = 0; i < ncmpt; ++i)
{
......
......@@ -74,6 +74,7 @@ namespace Nektar
typedef LibUtilities::NekFactory< std::string, Extrapolate,
const LibUtilities::SessionReaderSharedPtr& ,
Array<OneD, MultiRegions::ExpListSharedPtr>& ,
MultiRegions::ExpListSharedPtr& ,
const Array<OneD, int>& ,
const AdvectionTermSharedPtr& > ExtrapolateFactory;
......@@ -85,6 +86,7 @@ namespace Nektar
Extrapolate(
const LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields,
MultiRegions::ExpListSharedPtr pPressure,
const Array<OneD, int> pVel,
const AdvectionTermSharedPtr advObject);
......@@ -147,11 +149,16 @@ namespace Nektar
Array<OneD, NekDouble> &Q,
Array<OneD, const NekDouble> &Advection)=0;
void CalcPressureBCs(
void CalcNeumannPressureBCs(
const Array<OneD, const Array<OneD, NekDouble> > &fields,
const Array<OneD, const Array<OneD, NekDouble> > &N,
NekDouble kinvis);
void CalcOutflowBCs(
const Array<OneD, const Array<OneD, NekDouble> > &fields,
const Array<OneD, const Array<OneD, NekDouble> > &N,
NekDouble kinvis);
void RollOver(
Array<OneD, Array<OneD, NekDouble> > &input);
......@@ -166,6 +173,11 @@ namespace Nektar
Array<OneD, MultiRegions::ExpListSharedPtr> m_fields;
/// Pointer to field holding pressure field
MultiRegions::ExpListSharedPtr m_pressure;
/// int which identifies which components of m_fields contains the
/// velocity (u,v,w);
Array<OneD, int> m_velocity;
AdvectionTermSharedPtr m_advObject;
......@@ -225,14 +237,15 @@ namespace Nektar
/// data structure to old all the information regarding High order pressure BCs
Array<OneD, HBCInfo > m_HBCdata;
/// general standard element used to deal with HOPBC calculations
StdRegions::StdExpansionSharedPtr m_elmt;
/// wave number 2 pi k /Lz
Array<OneD, NekDouble> m_wavenumber;
/// minus Square of wavenumber
Array<OneD, NekDouble> m_negWavenumberSq;
/// Storage for current and previous velocity fields at the otuflow for high order outflow BCs
Array<OneD, Array<OneD, Array<OneD, NekDouble > > > m_outflowVel;
private:
static std::string def;
......
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