Commit c38f8b88 authored by Dave Moxey's avatar Dave Moxey
Browse files

Fixed some bugs, ADRSolver regression tests pass.

Added some initial code to wrap RiemannSolver calls.
parent 93fc0e14
......@@ -68,15 +68,15 @@ namespace Nektar
const LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
spaceDim = pFields[0]->GetCoordim(0);
m_spaceDim = pFields[0]->GetCoordim(0);
if (pSession->->DefinesSolverInfo("HOMOGENEOUS"))
if (pSession->DefinesSolverInfo("HOMOGENEOUS"))
{
std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
std::string HomoStr = pSession->GetSolverInfo("HOMOGENEOUS");
if (HomoStr == "HOMOGENEOUS1D" || HomoStr == "Homogeneous1D" ||
HomoStr == "1D" || HomoStr == "Homo1D")
{
spaceDim++;
m_spaceDim++;
}
else
{
......
......@@ -114,6 +114,27 @@ namespace Nektar
&Advection3DHomogeneous1D::ModifiedFluxVector, this);
m_planeCounter = 0;
// Override Riemann solver scalar and vector callbacks.
/*
map<string, RSScalarFuncType>::iterator it1;
map<string, RSScalarFuncType>::iterator it2;
m_scalars = m_riemann->GetScalars();
m_vectors = m_riemann->GetVectors();
for (it1 = scalars.begin(); it1 != scalars.end(); ++it1)
{
m_riemann->AddScalar(
boost::bind(&Advection3DHomogeneous1D::ModifiedRSScalar,
this, it1->first));
}
for (it2 = vectors.begin(); it2 != vectors.end(); ++it2)
{
m_riemann->AddVector(
boost::bind(&Advection3DHomogeneous1D::ModifiedRSVector,
this, it2->first));
}
*/
m_fluxVecStore = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >(
nConvectiveFields);
......@@ -233,5 +254,16 @@ namespace Nektar
// Increment the plane counter.
m_planeCounter = (m_planeCounter + 1) % m_numPlanes;
}
/*
const Array<OneD, const NekDouble>
&Advection3DHomogeneous1D::ModifiedRSScalar(string name)
{
if (m_rsPlaneNumber == 0)
{
}
}
*/
}
}
......@@ -57,29 +57,28 @@ namespace Nektar
protected:
Advection3DHomogeneous1D(std::string advType);
std::string m_advType;
SolverUtils::AdvectionSharedPtr m_planeAdv;
int m_numPoints;
int m_numPointsPlane;
int m_numPlanes;
int m_planeCounter;
Array<OneD, unsigned int> m_planes;
Array<OneD, unsigned int> m_planePos;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > m_fluxVecStore;
Array<OneD, Array<OneD, NekDouble> > m_inarrayPlane;
Array<OneD, Array<OneD, NekDouble> > m_outarrayPlane;
Array<OneD, MultiRegions::ExpListSharedPtr> m_fieldsPlane;
Array<OneD, Array<OneD, NekDouble> > m_advVelPlane;
std::string m_advType;
SolverUtils::AdvectionSharedPtr m_planeAdv;
int m_numPoints;
int m_numPointsPlane;
int m_numPlanes;
int m_planeCounter;
Array<OneD, unsigned int> m_planes;
Array<OneD, unsigned int> m_planePos;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > m_fluxVecStore;
Array<OneD, Array<OneD, NekDouble> > m_inarrayPlane;
Array<OneD, Array<OneD, NekDouble> > m_outarrayPlane;
Array<OneD, MultiRegions::ExpListSharedPtr> m_fieldsPlane;
Array<OneD, Array<OneD, NekDouble> > m_advVelPlane;
Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble> > > >
m_fluxVecPlane;
m_fluxVecPlane;
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_Advect(
const int nConvectiveFields,
const int nConvField,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
......@@ -89,6 +88,12 @@ namespace Nektar
void ModifiedFluxVector(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux);
/*
const Array<OneD, const NekDouble> &ModifiedRSScalar(
string name);
const Array<OneD, const Array<OneD, NekDouble> > &ModifiedRSVector(
string name);
*/
};
}
}
......
......@@ -87,8 +87,9 @@ namespace Nektar
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
v_SetupMetrics(pSession, pFields);
v_SetupCFunctions(pSession, pFields);
Advection::v_InitObject(pSession, pFields);
v_SetupMetrics (pSession, pFields);
v_SetupCFunctions (pSession, pFields);
}
/**
......@@ -815,7 +816,7 @@ namespace Nektar
fluxvector[i] =
Array<OneD, Array<OneD, NekDouble> >(m_spaceDim);
for (j = 0; j < nVel; ++j)
for (j = 0; j < m_spaceDim; ++j)
{
fluxvector[i][j] = Array<OneD, NekDouble>(nSolutionPts);
}
......
......@@ -48,6 +48,20 @@ namespace Nektar
}
/**
* @brief Initialise AdvectionNonConservative objects and store them
* before starting the time-stepping.
*
* @param pSession Pointer to session reader.
* @param pFields Pointer to fields.
*/
void AdvectionNonConservative::v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
Advection::v_InitObject(pSession, pFields);
}
void AdvectionNonConservative::v_Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
......
......@@ -54,9 +54,13 @@ namespace Nektar
protected:
AdvectionNonConservative();
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_Advect(
const int nConvectiveFields,
const int nConvective,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
......
......@@ -59,7 +59,7 @@ namespace Nektar
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
Advection::v_InitObject(pSession, pFields);
}
/**
......@@ -80,8 +80,7 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
int nVel = advVel.num_elements();
int nExpDim = fields[0]->GetCoordim(0);
int nDim = fields[0]->GetCoordim(0);
int nPointsTot = fields[0]->GetTotPoints();
int nCoeffs = fields[0]->GetNcoeffs();
int nTracePointsTot = fields[0]->GetTrace()->GetTotPoints();
......@@ -91,11 +90,12 @@ namespace Nektar
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > fluxvector(
nConvectiveFields);
// Allocate storage for flux vector F(u).
for (i = 0; i < nConvectiveFields; ++i)
{
fluxvector[i] =
Array<OneD, Array<OneD, NekDouble> >(nVel);
for (j = 0; j < nVel; ++j)
Array<OneD, Array<OneD, NekDouble> >(m_spaceDim);
for (j = 0; j < m_spaceDim; ++j)
{
fluxvector[i][j] = Array<OneD, NekDouble>(nPointsTot);
}
......@@ -111,7 +111,7 @@ namespace Nektar
{
tmp[i] = Array<OneD, NekDouble>(nCoeffs, 0.0);
for (j = 0; j < nExpDim; ++j)
for (j = 0; j < nDim; ++j)
{
fields[i]->IProductWRTDerivBase(j, fluxvector[i][j],
outarray[i]);
......
......@@ -56,11 +56,11 @@ namespace Nektar
AdvectionWeakDG();
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_Advect(
const int nConvective,
const int nConvective,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
......
......@@ -47,7 +47,6 @@
namespace Nektar
{
template <typename Dim, typename DataType>
class Array;
......@@ -59,7 +58,7 @@ namespace Nektar
const Array<OneD, const Array<OneD, NekDouble> >& ()> RSVecFuncType;
typedef boost::function<
NekDouble ()> RSParamFuncType;
class RiemannSolver
{
public:
......@@ -67,31 +66,60 @@ namespace Nektar
const Array<OneD, const Array<OneD, NekDouble> > &Fwd,
const Array<OneD, const Array<OneD, NekDouble> > &Bwd,
Array<OneD, Array<OneD, NekDouble> > &flux);
template<typename FuncPointerT, typename ObjectPointerT>
void AddScalar(std::string name,
FuncPointerT func,
template<typename FuncPointerT, typename ObjectPointerT>
void SetScalar(std::string name,
FuncPointerT func,
ObjectPointerT obj)
{
m_scalars[name] = boost::bind(func, obj);
}
void SetScalar(std::string name, RSScalarFuncType fp)
{
m_scalars[name] = fp;
}
template<typename FuncPointerT, typename ObjectPointerT>
void AddVector(std::string name,
FuncPointerT func,
template<typename FuncPointerT, typename ObjectPointerT>
void SetVector(std::string name,
FuncPointerT func,
ObjectPointerT obj)
{
m_vectors[name] = boost::bind(func, obj);
}
template<typename FuncPointerT, typename ObjectPointerT>
void AddParam(std::string name,
FuncPointerT func,
void SetVector(std::string name, RSVecFuncType fp)
{
m_vectors[name] = fp;
}
template<typename FuncPointerT, typename ObjectPointerT>
void SetParam(std::string name,
FuncPointerT func,
ObjectPointerT obj)
{
m_params[name] = boost::bind(func, obj);
}
void SetParam(std::string name, RSParamFuncType fp)
{
m_params[name] = fp;
}
std::map<std::string, RSScalarFuncType> &GetScalars()
{
return m_scalars;
}
std::map<std::string, RSVecFuncType> &GetVectors()
{
return m_vectors;
}
std::map<std::string, RSParamFuncType> &GetParams()
{
return m_params;
}
protected:
/// Indicates whether the Riemann solver requires a rotation to be
......@@ -129,8 +157,8 @@ namespace Nektar
bool CheckScalars(std::string name);
bool CheckVectors(std::string name);
bool CheckParams (std::string name);
};
};
/// A shared pointer to an EquationSystem object
typedef boost::shared_ptr<RiemannSolver> RiemannSolverSharedPtr;
/// Datatype of the NekFactory used to instantiate classes derived
......
......@@ -115,19 +115,17 @@ namespace Nektar
{
m_advection->SetFluxVector(
&UnsteadyAdvection::GetFluxVectorDeAlias, this);
}
else
{
m_advection->SetFluxVector(
&UnsteadyAdvection::GetFluxVector, this);
}
m_session->LoadSolverInfo(
"UpwindType", riemName, "Upwind");
m_riemannSolver = SolverUtils::
GetRiemannSolverFactory().CreateInstance(riemName);
m_riemannSolver->AddScalar(
m_riemannSolver->SetScalar(
"Vn", &UnsteadyAdvection::GetNormalVelocity, this);
m_advection->SetRiemannSolver(m_riemannSolver);
......
......@@ -88,7 +88,7 @@ namespace Nektar
m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
m_riemannSolver = SolverUtils::GetRiemannSolverFactory().
CreateInstance(riemName);
m_riemannSolver->AddScalar("Vn", &UnsteadyAdvectionDiffusion::
m_riemannSolver->SetScalar("Vn", &UnsteadyAdvectionDiffusion::
GetNormalVelocity, this);
m_advection->SetRiemannSolver(m_riemannSolver);
m_advection->InitObject (m_session, m_fields);
......@@ -195,16 +195,8 @@ namespace Nektar
m_traceVn, 1,
m_traceVn, 1);
}
if(m_planeNumber == n_planes - 1)
{
m_planeNumber = 0;
}
else
{
m_planeNumber = m_planeNumber + 1;
}
m_planeNumber = (m_planeNumber + 1) % n_planes;
}
else // For general case
{
......@@ -224,7 +216,8 @@ namespace Nektar
}
/* @brief Compute the right-hand side for the unsteady linear advection
/**
* @brief Compute the right-hand side for the unsteady linear advection
* diffusion problem.
*
* @param inarray Given fields.
......@@ -386,46 +379,15 @@ namespace Nektar
{
ASSERTL1(flux[0].num_elements() == m_velocity.num_elements(),
"Dimension of flux array and velocity array do not match");
int i , j;
int nq = physfield[0].num_elements();
if (m_expdim == 2 && m_HomogeneousType == eHomogeneous1D)
{
Array<OneD, Array<OneD, NekDouble> >
advVel_plane(m_velocity.num_elements());
int nPointsTot = m_fields[0]->GetTotPoints();
int nPointsTot_plane = m_fields[0]->GetPlane(0)->GetTotPoints();
int n_planes = nPointsTot/nPointsTot_plane;
for (int i = 0; i < m_velocity.num_elements(); ++i)
{
advVel_plane[i] = Array<OneD, NekDouble>(nPointsTot_plane, 0.0);
Vmath::Vcopy(nPointsTot_plane,
&m_velocity[i][m_planeNumber*nPointsTot_plane], 1,
&advVel_plane[i][0], 1);
}
for (int i = 0; i < flux.num_elements(); ++i)
{
for (j = 0; j < flux[0].num_elements(); ++j)
{
Vmath::Vmul(nq, physfield[i], 1, advVel_plane[j], 1,
flux[i][j], 1);
}
}
}
else
const int nq = m_fields[0]->GetNpoints();
for (int i = 0; i < flux.num_elements(); ++i)
{
for (i = 0; i < flux.num_elements(); ++i)
for (int j = 0; j < flux[0].num_elements(); ++j)
{
for (j = 0; j < flux[0].num_elements(); ++j)
{
Vmath::Vmul(nq, physfield[i], 1, m_velocity[j], 1,
flux[i][j], 1);
}
Vmath::Vmul(nq, physfield[i], 1, m_velocity[j], 1,
flux[i][j], 1);
}
}
}
......
......@@ -83,7 +83,7 @@ namespace Nektar
m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
m_riemannSolver = SolverUtils::GetRiemannSolverFactory().CreateInstance(riemName);
m_riemannSolver->AddScalar("Vn", &UnsteadyInviscidBurger::GetNormalVelocity, this);
m_riemannSolver->SetScalar("Vn", &UnsteadyInviscidBurger::GetNormalVelocity, this);
m_advection->SetRiemannSolver(m_riemannSolver);
m_advection->InitObject (m_session, m_fields);
......
......@@ -179,24 +179,24 @@ namespace Nektar
.CreateInstance("UpwindLDG");
// Setting up parameters for advection operator Riemann solver
m_riemannSolver->AddParam (
m_riemannSolver->SetParam (
"gamma",
&CompressibleFlowSystem::GetGamma, this);
m_riemannSolver->AddScalar(
m_riemannSolver->SetScalar(
"velLoc",
&CompressibleFlowSystem::GetVelLoc, this);
m_riemannSolver->AddVector(
m_riemannSolver->SetVector(
"N",
&CompressibleFlowSystem::GetNormals, this);
// Setting up parameters for diffusion operator Riemann solver
m_riemannSolverLDG->AddParam (
m_riemannSolverLDG->SetParam (
"gamma",
&CompressibleFlowSystem::GetGamma, this);
m_riemannSolverLDG->AddScalar(
m_riemannSolverLDG->SetScalar(
"velLoc",
&CompressibleFlowSystem::GetVelLoc, this);
m_riemannSolverLDG->AddVector(
m_riemannSolverLDG->SetVector(
"N",
&CompressibleFlowSystem::GetNormals, this);
......
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