Commit ea381797 authored by Chris Cantwell's avatar Chris Cantwell

Updated IncNavierStokes advection classes to inherit SolverUtils::Advection.

All tests pass, apart from substepping.
parent 8dbf5969
......@@ -39,7 +39,9 @@
namespace Nektar
{
{
namespace SolverUtils
{
AdvectionTermFactory& GetAdvectionTermFactory()
{
typedef Loki::SingletonHolder<AdvectionTermFactory,
......@@ -133,7 +135,7 @@ namespace Nektar
//////////////////////////////////////////////////////////////////////////////////////////////
void AdvectionTerm::DoAdvection(Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
void AdvectionTerm::DoAdvection(const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
const int nConvectiveFields,
const Array<OneD, int> &vel_loc,
const Array<OneD, const Array<OneD, NekDouble> > &pInarray,
......@@ -166,7 +168,7 @@ namespace Nektar
DoAdvection(pFields,velocity,pInarray,pOutarray,time,pWk);
}
void AdvectionTerm::DoAdvection(Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
void AdvectionTerm::DoAdvection(const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
const Array<OneD, const Array<OneD, NekDouble> > &velocity,
const Array<OneD, const Array<OneD, NekDouble> > &pInarray,
Array<OneD, Array<OneD, NekDouble> > &pOutarray,
......@@ -199,5 +201,5 @@ namespace Nektar
Vmath::Neg(nqtot,pOutarray[i],1);
}
}
}
} //end of namespace
......@@ -41,9 +41,11 @@
#include <LibUtilities/FFT/NektarFFT.h> // for NektarFFTSharedPtr
#include <SpatialDomains/MeshGraph.h> // for MeshGraphSharedPtr
#include <MultiRegions/ExpList.h> // for ExpListSharedPtr
#include <SolverUtils/Advection/Advection.h>
namespace Nektar
{
namespace SolverUtils
{
class AdvectionTerm;
......@@ -59,7 +61,7 @@ namespace Nektar
AdvectionTermFactory& GetAdvectionTermFactory();
/// Base class for the development of solvers.
class AdvectionTerm
class AdvectionTerm : public SolverUtils::Advection
{
public:
/// Destructor
......@@ -68,7 +70,7 @@ namespace Nektar
inline void InitObject();
/// Compute advection term
void DoAdvection(Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
void DoAdvection(const Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
const int nConvectiveFields,
const Array<OneD, int> &vel_loc,
const Array<OneD, const Array<OneD, NekDouble> > &pInarray,
......@@ -76,7 +78,7 @@ namespace Nektar
NekDouble m_time,
Array<OneD, NekDouble> &pWk = NullNekDouble1DArray);
void DoAdvection(Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
void DoAdvection(const Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
const Array<OneD, const Array<OneD, NekDouble> > &Velocity,
const Array<OneD, const Array<OneD, NekDouble> > &pInarray,
Array<OneD, Array<OneD, NekDouble> > &pOutarray,
......@@ -139,7 +141,20 @@ namespace Nektar
virtual void v_InitObject();
virtual void v_ComputeAdvectionTerm(Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
virtual void v_Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)
{
Array<OneD, int> vel_loc(fields.num_elements() - 1);
for (int i = 0; i < fields.num_elements(); ++i) vel_loc[i] = i;
DoAdvection(fields, nConvectiveFields, vel_loc, inarray, outarray, time);
}
virtual void v_ComputeAdvectionTerm(const Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
const Array<OneD, Array<OneD, NekDouble> > &pV,
const Array<OneD, const NekDouble> &pU,
Array<OneD, NekDouble> &pOutarray,
......@@ -155,7 +170,8 @@ namespace Nektar
{
v_InitObject();
}
} //end of namespace
}
} //endof namespace
#endif //NEKTAR_SOLVERS_AUXILIARY_ADRBASE_H
......@@ -14,7 +14,6 @@ namespace Nektar {
AdvectionSystem::AdvectionSystem(const LibUtilities::SessionReaderSharedPtr& pSession)
: UnsteadySystem(pSession)
{
}
AdvectionSystem::~AdvectionSystem()
......@@ -22,6 +21,11 @@ namespace Nektar {
}
void AdvectionSystem::v_InitObject()
{
UnsteadySystem::v_InitObject();
}
//////////////////////////////////////////////////////////////////////////////////////////////
void AdvectionSystem::UpdateBaseFlow(const Array<OneD, Array<OneD, NekDouble> > &inarray)
{
......@@ -30,4 +34,4 @@ namespace Nektar {
//////////////////////////////////////////////////////////////////////////////////////////////
}
}
\ No newline at end of file
}
......@@ -50,6 +50,8 @@ namespace Nektar {
AdvectionSystem(const LibUtilities::SessionReaderSharedPtr &pSession);
virtual ~AdvectionSystem();
virtual void v_InitObject();
AdvectionSharedPtr GetAdvObject()
{
return m_advObject;
......@@ -80,4 +82,4 @@ namespace Nektar {
}
}
#endif
\ No newline at end of file
#endif
......@@ -49,7 +49,7 @@
namespace Nektar
{
string AdjointAdvection::className = GetAdvectionTermFactory().RegisterCreatorFunction("Adjoint", AdjointAdvection::create);
string AdjointAdvection::className = SolverUtils::GetAdvectionFactory().RegisterCreatorFunction("Adjoint", AdjointAdvection::create);
/**
* Constructor. Creates ...
......@@ -58,19 +58,22 @@ namespace Nektar
* \param
*/
AdjointAdvection::AdjointAdvection(
const LibUtilities::SessionReaderSharedPtr& pSession,
const SpatialDomains::MeshGraphSharedPtr& pGraph):
AdvectionTerm(pSession, pGraph)
AdjointAdvection::AdjointAdvection():
Advection()
{
}
void AdjointAdvection::v_InitObject()
{
AdvectionTerm::v_InitObject();
m_boundaryConditions = MemoryManager<SpatialDomains::BoundaryConditions>
::AllocateSharedPtr(m_session, m_graph);
void AdjointAdvection::v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
Advection::v_InitObject(pSession, pFields);
m_session = pSession;
m_boundaryConditions = MemoryManager<SpatialDomains::BoundaryConditions>
::AllocateSharedPtr(m_session, pFields[0]->GetGraph());
//Setting parameters for homogeneous problems
m_HomoDirec = 0;
......@@ -79,7 +82,10 @@ namespace Nektar
m_SingleMode =false;
m_HalfMode =false;
m_MultipleModes =false;
m_spacedim = pFields[0]->GetGraph()->GetSpaceDimension();
m_expdim = pFields[0]->GetGraph()->GetMeshDimension();
m_CoeffState = MultiRegions::eLocal;
if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
{
std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
......@@ -92,6 +98,8 @@ namespace Nektar
m_LhomZ = m_session->GetParameter("LZ");
m_HomoDirec = 1;
m_homogen_dealiasing = pSession->DefinesSolverInfo("DEALIASING");
if(m_session->DefinesSolverInfo("ModeType"))
{
m_session->MatchSolverInfo("ModeType","SingleMode",m_SingleMode,false);
......@@ -189,7 +197,7 @@ namespace Nektar
m_projectionType = MultiRegions::eGalerkin;
}
SetUpBaseFields(m_graph);
SetUpBaseFields(pFields[0]->GetGraph());
ASSERTL0(m_session->DefinesFunction("BaseFlow"),
"Base flow must be defined for linearised forms.");
string file = m_session->GetFunctionFilename("BaseFlow", 0);
......@@ -217,7 +225,7 @@ namespace Nektar
if (m_session->GetFunctionType("BaseFlow", m_session->GetVariable(0))
== LibUtilities::eFunctionTypeFile)
{
ImportFldBase(file,m_graph,1);
ImportFldBase(file,pFields[0]->GetGraph(),1);
}
//analytic base flow
......@@ -254,7 +262,57 @@ namespace Nektar
}
void AdjointAdvection::SetUpBaseFields(SpatialDomains::MeshGraphSharedPtr &mesh)
void AdjointAdvection::v_Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)
{
int i;
int nqtot = fields[0]->GetTotPoints();
Array<OneD, Array<OneD, NekDouble> > velocity(nConvectiveFields);
ASSERTL1(nConvectiveFields == inarray.num_elements(),"Number of convective fields and Inarray are not compatible");
for(i = 0; i < nConvectiveFields; ++i)
{
if(fields[i]->GetWaveSpace() && !m_SingleMode && !m_HalfMode)
{
velocity[i] = Array<OneD, NekDouble>(nqtot,0.0);
fields[i]->HomogeneousBwdTrans(inarray[i],velocity[i]);
}
else
{
velocity[i] = inarray[i];
}
}
Array<OneD, NekDouble > Deriv;
// // Set up Derivative work space;
// if(pWk.num_elements())
// {
// ASSERTL0(pWk.num_elements() >= nqtot*VelDim,"Workspace is not sufficient");
// Deriv = pWk;
// }
// else
// {
Deriv = Array<OneD, NekDouble> (nqtot*nConvectiveFields);
// }
for(i=0; i< nConvectiveFields; ++i)
{
v_ComputeAdvectionTerm(fields,velocity,inarray[i],outarray[i],i,time,Deriv);
Vmath::Neg(nqtot,outarray[i],1);
}
}
void AdjointAdvection::SetUpBaseFields(SpatialDomains::MeshGraphSharedPtr mesh)
{
int nvariables = m_session->GetVariables().size();
int i;
......@@ -275,7 +333,7 @@ namespace Nektar
for(i = 0 ; i < m_base.num_elements(); i++)
{
m_base[i] = MemoryManager<MultiRegions::ContField3DHomogeneous2D>
::AllocateSharedPtr(m_session,BkeyY,BkeyZ,m_LhomY,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
::AllocateSharedPtr(m_session,BkeyY,BkeyZ,m_LhomY,m_LhomZ,m_useFFT,m_homogen_dealiasing,mesh,m_session->GetVariable(i));
}
}
......@@ -303,7 +361,7 @@ namespace Nektar
for(i = 0 ; i < m_base.num_elements(); i++)
{
m_base[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,mesh,m_session->GetVariable(i));
}
}
......@@ -316,7 +374,7 @@ namespace Nektar
for(i = 0 ; i < m_base.num_elements(); i++)
{
m_base[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,mesh,m_session->GetVariable(i));
}
}
......@@ -329,7 +387,7 @@ namespace Nektar
for(i = 0 ; i < m_base.num_elements(); i++)
{
m_base[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,mesh,m_session->GetVariable(i));
m_base[i]->SetWaveSpace(false);
}
......@@ -397,7 +455,7 @@ namespace Nektar
for(i = 0 ; i < m_base.num_elements(); i++)
{
m_base[i] = MemoryManager<MultiRegions::DisContField3DHomogeneous2D>
::AllocateSharedPtr(m_session,BkeyY,BkeyZ,m_LhomY,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
::AllocateSharedPtr(m_session,BkeyY,BkeyZ,m_LhomY,m_LhomZ,m_useFFT,m_homogen_dealiasing,mesh,m_session->GetVariable(i));
}
}
else
......@@ -422,7 +480,7 @@ namespace Nektar
for(i = 0 ; i < m_base.num_elements(); i++)
{
m_base[i] = MemoryManager<MultiRegions::DisContField3DHomogeneous1D>
::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,mesh,m_session->GetVariable(i));
}
......@@ -552,9 +610,9 @@ namespace Nektar
if(m_session->DefinesParameter("N_slices"))
{
m_nConvectiveFields = m_base.num_elements()-1;
int n = m_base.num_elements()-1;
for(int i=0; i<m_nConvectiveFields;++i)
for(int i=0; i<n;++i)
{
Vmath::Vcopy(nqtot, &m_base[i]->GetPhys()[0], 1, &m_interp[i][cnt*nqtot], 1);
......@@ -566,7 +624,7 @@ namespace Nektar
//Evaluation of the advective terms
void AdjointAdvection::v_ComputeAdvectionTerm(
Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
const Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
const Array<OneD, Array<OneD, NekDouble> > &pVelocity,
const Array<OneD, const NekDouble> &pU,
Array<OneD, NekDouble> &pOutarray,
......@@ -574,7 +632,7 @@ namespace Nektar
NekDouble m_time,
Array<OneD, NekDouble> &pWk)
{
int ndim = m_nConvectiveFields;
int ndim = pVelocity.num_elements();
int nPointsTot = pFields[0]->GetNpoints();
Array<OneD, NekDouble> grad0,grad1,grad2;
......@@ -601,7 +659,7 @@ namespace Nektar
if (m_session->GetFunctionType("BaseFlow", 0)
== LibUtilities::eFunctionTypeFile)
{
for(int i=0; i<m_nConvectiveFields;++i)
for(int i=0; i<ndim;++i)
{
UpdateBase(m_slices,m_interp[i],m_base[i]->UpdatePhys(),m_time,m_period);
}
......@@ -851,7 +909,7 @@ namespace Nektar
{
char chkout[16] = "";
sprintf(chkout, "%d", i);
ImportFldBase(file+"_"+chkout+".bse",m_graph,i);
ImportFldBase(file+"_"+chkout+".bse",m_base[0]->GetGraph(),i);
}
......
......@@ -36,8 +36,8 @@
#ifndef NEKTAR_SOLVERS_ADJOINTADVECTION_H
#define NEKTAR_SOLVERS_ADJOINTADVECTION_H
///#include <IncNavierStokesSolver/AdvectionTerms/AdvectionTerm.h>
#include <SolverUtils/Advection/AdvectionTerm.h>
#include <SolverUtils/Advection/Advection.h>
#include <LibUtilities/FFT/NektarFFT.h>
//#define TIMING
......@@ -51,7 +51,7 @@ namespace Nektar
{
class AdjointAdvection: public AdvectionTerm
class AdjointAdvection: public SolverUtils::Advection
{
enum FloquetMatType
{
......@@ -69,17 +69,20 @@ namespace Nektar
friend class MemoryManager<AdjointAdvection>;
/// Creates an instance of this class
static AdvectionTermSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession,
const SpatialDomains::MeshGraphSharedPtr& pGraph) {
AdvectionTermSharedPtr p = MemoryManager<AdjointAdvection>::AllocateSharedPtr(pSession, pGraph);
p->InitObject();
static SolverUtils::AdvectionSharedPtr create(std::string) {
SolverUtils::AdvectionSharedPtr p = MemoryManager<AdjointAdvection>::AllocateSharedPtr();
return p;
}
/// Name of class
static std::string className;
protected:
LibUtilities::SessionReaderSharedPtr m_session;
MultiRegions::ProjectionType m_projectionType;
int m_spacedim;
int m_expdim;
//Storage of the base flow
Array<OneD, MultiRegions::ExpListSharedPtr> m_base;
......@@ -97,20 +100,30 @@ namespace Nektar
bool m_SingleMode; ///< flag to determine if use single mode or not
bool m_HalfMode; ///< flag to determine if use half mode or not
bool m_MultipleModes; ///< flag to determine if use multiple mode or not
bool m_homogen_dealiasing;
MultiRegions::CoeffState m_CoeffState;
DNekBlkMatSharedPtr GetFloquetBlockMatrix(FloquetMatType mattype, bool UseContCoeffs = false) const;
DNekBlkMatSharedPtr GenFloquetBlockMatrix(FloquetMatType mattype, bool UseContCoeffs = false) const;
FloquetBlockMatrixMapShPtr m_FloquetBlockMat;
AdjointAdvection(
const LibUtilities::SessionReaderSharedPtr& pSession,
const SpatialDomains::MeshGraphSharedPtr& pGraph);
AdjointAdvection();
virtual ~AdjointAdvection();
virtual void v_InitObject();
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time);
void SetUpBaseFields(SpatialDomains::MeshGraphSharedPtr &mesh);
void SetUpBaseFields(SpatialDomains::MeshGraphSharedPtr mesh);
void UpdateBase(const NekDouble m_slices,
Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
......@@ -140,7 +153,7 @@ namespace Nektar
private:
//Function for the evaluation of the Adjoint advective terms
virtual void v_ComputeAdvectionTerm(
Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
const Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
const Array<OneD, Array<OneD, NekDouble> > &pV,
const Array<OneD, const NekDouble> &pU,
Array<OneD, NekDouble> &pOutarray,
......
......@@ -36,8 +36,8 @@
#ifndef NEKTAR_SOLVERS_LINEARISEDADVECTION_H
#define NEKTAR_SOLVERS_LINEARISEDADVECTION_H
///#include <IncNavierStokesSolver/AdvectionTerms/AdvectionTerm.h>
#include <SolverUtils/Advection/AdvectionTerm.h>
#include <SolverUtils/Advection/Advection.h>
#include <LibUtilities/FFT/NektarFFT.h>
//#define TIMING
......@@ -51,7 +51,7 @@ namespace Nektar
{
class LinearisedAdvection: public AdvectionTerm
class LinearisedAdvection: public SolverUtils::Advection
{
enum FloquetMatType
{
......@@ -69,17 +69,20 @@ namespace Nektar
friend class MemoryManager<LinearisedAdvection>;
/// Creates an instance of this class
static AdvectionTermSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession,
const SpatialDomains::MeshGraphSharedPtr& pGraph) {
AdvectionTermSharedPtr p = MemoryManager<LinearisedAdvection>::AllocateSharedPtr(pSession, pGraph);
p->InitObject();
static SolverUtils::AdvectionSharedPtr create(std::string) {
SolverUtils::AdvectionSharedPtr p = MemoryManager<LinearisedAdvection>::AllocateSharedPtr();
return p;
}
/// Name of class
static std::string className;
protected:
LibUtilities::SessionReaderSharedPtr m_session;
MultiRegions::ProjectionType m_projectionType;
int m_spacedim;
int m_expdim;
//Storage of the base flow
Array<OneD, MultiRegions::ExpListSharedPtr> m_base;
......@@ -99,20 +102,31 @@ namespace Nektar
bool m_SingleMode; ///< flag to determine if use single mode or not
bool m_HalfMode; ///< flag to determine if use half mode or not
bool m_MultipleModes; ///< flag to determine if use multiple mode or not
bool m_homogen_dealiasing;
MultiRegions::CoeffState m_CoeffState;
DNekBlkMatSharedPtr GetFloquetBlockMatrix(FloquetMatType mattype, bool UseContCoeffs = false) const;
DNekBlkMatSharedPtr GenFloquetBlockMatrix(FloquetMatType mattype, bool UseContCoeffs = false) const;
FloquetBlockMatrixMapShPtr m_FloquetBlockMat;
LinearisedAdvection(const LibUtilities::SessionReaderSharedPtr& pSession,
const SpatialDomains::MeshGraphSharedPtr& pGraph);
LinearisedAdvection();
virtual ~LinearisedAdvection();
virtual void v_InitObject();
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time);
void SetUpBaseFields(SpatialDomains::MeshGraphSharedPtr &mesh);
void SetUpBaseFields(SpatialDomains::MeshGraphSharedPtr mesh);
void UpdateBase(const NekDouble m_slices,
Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
......@@ -137,9 +151,10 @@ namespace Nektar
Array<OneD, std::string> &variables);
private:
//Function for the evaluation of the linearised advective terms
virtual void v_ComputeAdvectionTerm(
Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
const Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
const Array<OneD, Array<OneD, NekDouble> > &pV,
const Array<OneD, const NekDouble> &pU,
Array<OneD, NekDouble> &pOutarray,
......
......@@ -37,8 +37,8 @@
namespace Nektar
{
string NavierStokesAdvection::className = GetAdvectionTermFactory().RegisterCreatorFunction("Convective", NavierStokesAdvection::create);
string NavierStokesAdvection::className2 = GetAdvectionTermFactory().RegisterCreatorFunction("NonConservative", NavierStokesAdvection::create);
string NavierStokesAdvection::className = SolverUtils::GetAdvectionFactory().RegisterCreatorFunction("Convective", NavierStokesAdvection::create);
string NavierStokesAdvection::className2 = SolverUtils::GetAdvectionFactory().RegisterCreatorFunction("NonConservative", NavierStokesAdvection::create);
/**
* Constructor. Creates ...
......@@ -47,10 +47,8 @@ namespace Nektar
* \param
*/
NavierStokesAdvection::NavierStokesAdvection(
const LibUtilities::SessionReaderSharedPtr& pSession,
const SpatialDomains::MeshGraphSharedPtr& pGraph):
AdvectionTerm(pSession, pGraph)
NavierStokesAdvection::NavierStokesAdvection():
Advection()
{
......@@ -60,12 +58,79 @@ namespace Nektar
{
}
void NavierStokesAdvection::v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
m_CoeffState = MultiRegions::eLocal;
m_homogen_dealiasing = pSession->DefinesSolverInfo("dealiasing");
pSession->MatchSolverInfo("SPECTRALHPDEALIASING","True",m_specHP_dealiasing,false);
if(m_specHP_dealiasing == false)
{
pSession->MatchSolverInfo("SPECTRALHPDEALIASING","On",m_specHP_dealiasing,false);
}
pSession->MatchSolverInfo("ModeType","SingleMode",m_SingleMode,false);
pSession->MatchSolverInfo("ModeType","HalfMode",m_HalfMode,false);
Advection::v_InitObject(pSession, pFields);
}
//Advection function
void NavierStokesAdvection::v_Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)
{
int i;
int nqtot = fields[0]->GetTotPoints();
Array<OneD, Array<OneD, NekDouble> > velocity(nConvectiveFields);
ASSERTL1(nConvectiveFields == inarray.num_elements(),"Number of convective fields and Inarray are not compatible");
for(i = 0; i < nConvectiveFields; ++i)
{
if(fields[i]->GetWaveSpace() && !m_SingleMode && !m_HalfMode)
{
velocity[i] = Array<OneD, NekDouble>(nqtot,0.0);
fields[i]->HomogeneousBwdTrans(inarray[i],velocity[i]);
}
else
{
velocity[i] = inarray[i];
}
}
Array<OneD, NekDouble > Deriv;
// // Set up Derivative work space;
// if(pWk.num_elements())
// {
// ASSERTL0(pWk.num_elements() >= nqtot*VelDim,"Workspace is not sufficient");
// Deriv = pWk;
// }
// else
// {
Deriv = Array<OneD, NekDouble> (nqtot*nConvectiveFields);
// }
for(i=0; i< nConvectiveFields; ++i)
{
v_ComputeAdvectionTerm(fields,velocity,inarray[i],outarray[i],i,time,Deriv);
Vmath::Neg(nqtot,outarray[i],1);
}
}
//Evaluation of the advective terms
void NavierStokesAdvection::v_ComputeAdvectionTerm(
Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
const Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
const Array<OneD, Array<OneD, NekDouble> > &pV,
const Array<OneD, const NekDouble> &pU,
Array<OneD, NekDouble> &pOutarray,
......