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

ADRSolver diffusion working, compressible solver still not working

parent 03586dce
......@@ -112,9 +112,9 @@ namespace Nektar
// Override Riemann solver scalar and vector callbacks.
map<string, RSScalarFuncType>::iterator it1;
//map<string, RSScalarFuncType>::iterator it2;
map<string, RSVecFuncType>::iterator it2;
map<string, RSScalarFuncType> scalars = m_riemann->GetScalars();
//m_vectors = m_riemann->GetVectors();
map<string, RSVecFuncType> vectors = m_riemann->GetVectors();
for (it1 = scalars.begin(); it1 != scalars.end(); ++it1)
{
......@@ -123,14 +123,12 @@ namespace Nektar
m_riemann->SetScalar(it1->first, &HomoRSScalar::Exec, tmp);
}
/*
for (it2 = vectors.begin(); it2 != vectors.end(); ++it2)
{
m_riemann->SetVector(it2->first,
boost::bind(&Advection3DHomogeneous1D::ModifiedRSVector,
this, it2->first));
boost::shared_ptr<HomoRSVector> tmp = MemoryManager<HomoRSVector>
::AllocateSharedPtr(it2->second, m_numPlanes);
m_riemann->SetVector(it2->first, &HomoRSVector::Exec, tmp);
}
*/
m_fluxVecStore = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >(
nConvectiveFields);
......
......@@ -37,6 +37,7 @@
#define NEKTAR_SOLVERUTILS_ADVECTION3DHOMOGENEOUS1D
#include <SolverUtils/Advection/Advection.h>
#include <SolverUtils/Advection/HomogeneousRSScalar.hpp>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
......@@ -44,44 +45,6 @@ namespace Nektar
{
namespace SolverUtils
{
/**
* @brief Wrapper class for Riemann solver scalars.
*/
class HomoRSScalar
{
public:
HomoRSScalar(RSScalarFuncType func,
int nPlanes)
: m_planeNumber(0),
m_numPlanes (nPlanes),
m_func (func),
m_tmp ()
{
}
const Array<OneD, const NekDouble>& Exec()
{
if (m_planeNumber == 0)
{
m_tmp = m_func();
}
const int nPts = m_tmp.num_elements() / m_numPlanes;
const int offset = m_planeNumber * nPts;
Array<OneD, NekDouble> tmp(nPts, m_tmp + offset);
m_planeNumber = (m_planeNumber + 1) % m_numPlanes;
return tmp;
}
private:
int m_planeNumber;
int m_numPlanes;
RSScalarFuncType m_func;
Array<OneD, const NekDouble> m_tmp;
};
class Advection3DHomogeneous1D : public Advection
{
public:
......
......@@ -64,18 +64,64 @@ namespace Nektar
const int nPts = m_tmp.num_elements() / m_numPlanes;
const int offset = m_planeNumber * nPts;
Array<OneD, NekDouble> tmp(nPts, m_tmp + offset);
m_tmp2 = Array<OneD, NekDouble>(nPts, m_tmp + offset);
m_planeNumber = (m_planeNumber + 1) % m_numPlanes;
return tmp;
return m_tmp2;
}
private:
RSScalarFuncType m_func;
int m_planeNumber;
int m_numPlanes;
RSScalarFuncType m_func;
Array<OneD, const NekDouble> m_tmp;
Array<OneD, const NekDouble> m_tmp2;
};
/**
* @brief Wrapper class for Riemann solver scalars.
*/
class HomoRSVector
{
public:
HomoRSVector(RSVecFuncType func,
int nPlanes)
: m_func (func),
m_planeNumber(0),
m_numPlanes (nPlanes),
m_tmp ()
{
}
const Array<OneD, const Array<OneD, NekDouble> >& Exec()
{
if (m_planeNumber == 0)
{
m_tmp = m_func();
}
const int nDim = m_tmp.num_elements();
const int nPts = m_tmp[0].num_elements() / m_numPlanes;
const int offset = m_planeNumber * nPts;
m_tmp2 = Array<OneD, Array<OneD, NekDouble> >(nDim);
for (int i = 0; i < m_tmp.num_elements(); ++i)
{
m_tmp2[i] = Array<OneD, NekDouble>(nPts, m_tmp[i] + offset);
}
m_planeNumber = (m_planeNumber + 1) % m_numPlanes;
return m_tmp2;
}
private:
RSVecFuncType m_func;
int m_planeNumber;
int m_numPlanes;
Array<OneD, Array<OneD, NekDouble> > m_tmp;
Array<OneD, Array<OneD, NekDouble> > m_tmp2;
};
}
}
......@@ -126,7 +126,8 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)=0;
virtual void v_SetHomoDerivs(Array<OneD, Array<OneD, NekDouble> > &deriv)
virtual void v_SetHomoDerivs(
Array<OneD, Array<OneD, NekDouble> > &deriv)
{
}
......
......@@ -241,7 +241,11 @@ namespace Nektar
m_numPointsPlane, tmp2 = outarray[j] + m_planePos[i]);
}
m_planeDiff->SetHomoDerivs(m_homoDerivPlane[i]);
if (m_fluxVectorNS)
{
m_planeDiff->SetHomoDerivs(m_homoDerivPlane[i]);
}
m_planeDiff->Diffuse(nConvectiveFields,
m_fieldsPlane,
m_inarrayPlane,
......
......@@ -150,7 +150,7 @@ namespace Nektar
{
const Array<OneD, const Array<OneD, NekDouble> > &normals =
m_vectors["N"]();
const Array<OneD, NekDouble> &velLoc = m_scalars["velLoc"]();
const Array<OneD, NekDouble> &velLoc = m_auxiliary["velLoc"]();
switch(normals.num_elements())
{
......@@ -249,7 +249,7 @@ namespace Nektar
{
const Array<OneD, const Array<OneD, NekDouble> > &normals =
m_vectors["N"]();
const Array<OneD, NekDouble> &velLoc = m_scalars["velLoc"]();
const Array<OneD, NekDouble> &velLoc = m_auxiliary["velLoc"]();
switch(normals.num_elements())
{
......
......@@ -79,7 +79,7 @@ namespace Nektar
{
m_scalars[name] = fp;
}
template<typename FuncPointerT, typename ObjectPointerT>
void SetVector(std::string name,
FuncPointerT func,
......@@ -105,6 +105,14 @@ namespace Nektar
{
m_params[name] = fp;
}
template<typename FuncPointerT, typename ObjectPointerT>
void SetAuxiliary(std::string name,
FuncPointerT func,
ObjectPointerT obj)
{
m_auxiliary[name] = boost::bind(func, obj);
}
std::map<std::string, RSScalarFuncType> &GetScalars()
{
......@@ -131,6 +139,8 @@ namespace Nektar
std::map<std::string, RSVecFuncType> m_vectors;
/// Map of parameter function types.
std::map<std::string, RSParamFuncType > m_params;
/// Map of scalar function types.
std::map<std::string, RSScalarFuncType> m_auxiliary;
/// Rotation matrices for each trace quadrature point.
Array<OneD, Array<OneD, NekDouble> > m_rotMat;
/// Rotation storage
......
......@@ -52,7 +52,7 @@ namespace Nektar
const LibUtilities::SessionReaderSharedPtr& pSession)
: UnsteadySystem(pSession)
{
m_planeNumber = 0;
}
/**
......@@ -179,7 +179,7 @@ namespace Nektar
// Setting up parameters for advection operator Riemann solver
m_riemannSolver->SetParam (
"gamma", &CompressibleFlowSystem::GetGamma, this);
m_riemannSolver->SetScalar(
m_riemannSolver->SetAuxiliary(
"velLoc", &CompressibleFlowSystem::GetVelLoc, this);
m_riemannSolver->SetVector(
"N", &CompressibleFlowSystem::GetNormals, this);
......@@ -187,7 +187,7 @@ namespace Nektar
// Setting up parameters for diffusion operator Riemann solver
m_riemannSolverLDG->SetParam (
"gamma", &CompressibleFlowSystem::GetGamma, this);
m_riemannSolverLDG->SetScalar(
m_riemannSolverLDG->SetAuxiliary(
"velLoc", &CompressibleFlowSystem::GetVelLoc, this);
m_riemannSolverLDG->SetVector(
"N", &CompressibleFlowSystem::GetNormals, this);
......@@ -1020,9 +1020,41 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > fields(nVariables);
// Reorder storage to list time-integrated fields first
for (i = 0; i < nVariables; ++i)
// For 3DHomogenoeus1D
if (m_expdim == 2 && m_HomogeneousType == eHomogeneous1D)
{
int nSolutionPts = m_fields[0]->GetTotPoints();
int nSolutionPtsPlane = m_fields[0]->GetPlane(0)->GetTotPoints();
int nPlanes = nSolutionPts/nSolutionPtsPlane;
Array<OneD, Array<OneD, NekDouble> > fields_temp(nVariables);
for (i = 0; i < nVariables; ++i)
{
fields[i] = Array<OneD, NekDouble>(nPts, 0.0);
fields_temp[i] = m_fields[i]->UpdatePhys();
Vmath::Vcopy(nPts,
&fields_temp[i][m_planeNumber*nSolutionPtsPlane], 1,
&fields[i][0], 1);
}
if(m_planeNumber == nPlanes - 1)
{
m_planeNumber = 0;
}
else
{
m_planeNumber = m_planeNumber + 1;
}
}
else // For general case
{
fields[i] = m_fields[i]->UpdatePhys();
for (i = 0; i < nVariables; ++i)
{
fields[i] = m_fields[i]->UpdatePhys();
}
}
// Thermodynamic related quantities
......@@ -1392,17 +1424,61 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > fields_interp(nVariables);
// Reorder storage to list time-integrated fields first
for (i = 0; i < nVariables; ++i)
// For 3DHomogenoeus1D
if (m_expdim == 2 && m_HomogeneousType == eHomogeneous1D)
{
fields[i] = m_fields[i]->UpdatePhys();
fields_interp[i] = Array<OneD, NekDouble> (nPts);
int nSolutionPts = m_fields[0]->GetTotPoints();
int nSolutionPtsPlane = m_fields[0]->GetPlane(0)->GetTotPoints();
int nPlanes = nSolutionPts/nSolutionPtsPlane;
Array<OneD, Array<OneD, NekDouble> > fields_temp(nVariables);
for (i = 0; i < nVariables; ++i)
{
fields[i] = Array<OneD, NekDouble>(nPts, 0.0);
fields_temp[i] = m_fields[i]->UpdatePhys();
Vmath::Vcopy(nPts,
&fields_temp[i][m_planeNumber*nSolutionPtsPlane], 1,
&fields[i][0], 1);
fields_interp[i] = Array<OneD, NekDouble> (nPts);
}
if(m_planeNumber == nPlanes - 1)
{
m_planeNumber = 0;
}
else
{
m_planeNumber = m_planeNumber + 1;
}
}
else // For general case
{
for (i = 0; i < nVariables; ++i)
{
fields[i] = m_fields[i]->UpdatePhys();
fields_interp[i] = Array<OneD, NekDouble> (nPts);
}
}
for (i = 0; i < nVariables; ++i)
{
// Interpolation to higher space
m_fields[0]->PhysInterp1DScaled(OneDptscale,fields[i],
// For 3DHomogenoeus1D
if (m_expdim == 2 && m_HomogeneousType == eHomogeneous1D)
{
m_fields[0]->GetPlane(0)->PhysInterp1DScaled(
OneDptscale,fields[i], fields_interp[i]);
}
else // For general case
{
m_fields[0]->PhysInterp1DScaled(OneDptscale,fields[i],
fields_interp[i] );
}
}
for (i = 0; i < variables_phys; ++i)
......
......@@ -106,6 +106,10 @@ namespace Nektar
NekDouble m_Cp;
NekDouble m_Prandtl;
// Plane (used only for Discontinous projection
// with 3DHomogenoeus1D expansion)
int m_planeNumber;
CompressibleFlowSystem(
const LibUtilities::SessionReaderSharedPtr& pSession);
......
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