Commit 18546a0f authored by Kilian Lackhove's avatar Kilian Lackhove

RiemannSolver: Change the "velLoc" 1D-array to a 2D "vecLocs" array, that...

RiemannSolver: Change the "velLoc" 1D-array to a 2D "vecLocs" array, that holds the component positions for each directed quantity that needs to be rotated.

This also replaces the m_auxiliary member with m_auxScal and m_auxVec, so now auxiliary scalars as well as aux vectors can be stored.
parent 4642564b
......@@ -89,7 +89,8 @@ namespace Nektar
* velocity field is rotated to the normal direction to perform
* dimensional splitting, and the resulting fluxes are rotated back to
* the Cartesian directions before being returned. For the Rotation to
* work, the normal vectors "N" must be set via the SetVector() method.
* work, the normal vectors "N" and the location of the vector
* components in Fwd "vecLocs"must be set via the SetAuxVec() method.
*
* @param Fwd Forwards trace space.
* @param Bwd Backwards trace space.
......@@ -102,7 +103,12 @@ namespace Nektar
{
if (m_requiresRotation)
{
ASSERTL0(CheckVectors("N"), "N not defined.");
ASSERTL0(CheckVectors("N"), "N not defined.");
ASSERTL0(CheckAuxVec("vecLocs"), "vecLocs not defined.");
const Array<OneD, const Array<OneD, NekDouble> > normals =
m_vectors["N"]();
const Array<OneD, const Array<OneD, NekDouble> > vecLocs =
m_auxVec["vecLocs"]();
int nFields = Fwd .num_elements();
int nPts = Fwd[0].num_elements();
......@@ -120,15 +126,11 @@ namespace Nektar
}
}
}
const Array<OneD, const Array<OneD, NekDouble> > normals =
m_vectors["N"]();
rotateToNormal (Fwd, normals, m_rotStorage[0]);
rotateToNormal (Bwd, normals, m_rotStorage[1]);
v_Solve (m_rotStorage[0], m_rotStorage[1],
m_rotStorage[2]);
rotateFromNormal(m_rotStorage[2], normals, flux);
rotateToNormal(Fwd, normals, vecLocs, m_rotStorage[0]);
rotateToNormal(Bwd, normals, vecLocs, m_rotStorage[1]);
v_Solve(m_rotStorage[0], m_rotStorage[1], m_rotStorage[2]);
rotateFromNormal(m_rotStorage[2], normals, vecLocs, flux);
}
else
{
......@@ -142,12 +144,10 @@ namespace Nektar
* This function performs a rotation of a vector so that the first
* component aligns with the trace normal direction.
*
* The vectors components are stored in inarray. Their locations need
* to be specified in the "velLoc" auxiliary array. If the inarray holds
* multiple vectors, their locations must be concatenated, e.g.
* "0,1,4,5" for two 2D vectors with the first vectors component locations
* stored in inarray[0] and inarray[1] and the second vectors component
* locations in inarray[4] and inarray[5].
* The vectors components are stored in inarray. Their locations must
* be specified in the "vecLocs" array. vecLocs[0] contains the locations
* of the first vectors components, vecLocs[1] those of the second and
* so on.
*
* In 2D, this is accomplished through the transform:
*
......@@ -160,26 +160,21 @@ namespace Nektar
void RiemannSolver::rotateToNormal(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
const Array<OneD, const Array<OneD, NekDouble> > &normals,
const Array<OneD, const Array<OneD, NekDouble> > &vecLocs,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
ASSERTL0(CheckAuxiliary("velLoc"), "velLoc not defined.");
const Array<OneD, NekDouble> &velLoc = m_auxiliary["velLoc"]();
int vecCount = velLoc.num_elements() / normals.num_elements();
int vecSize = normals.num_elements();
ASSERTL0(velLoc.num_elements() % normals.num_elements() == 0,
"velLoc and normals mismatch");
for (int i = 0; i < inarray.num_elements(); ++i)
{
Vmath::Vcopy(inarray[i].num_elements(), inarray[i], 1,
outarray[i], 1);
}
for (int i = 0; i < vecCount; i++)
for (int i = 0; i < vecLocs.num_elements(); i++)
{
switch(vecSize)
ASSERTL0(vecLocs[i].num_elements() == normals.num_elements(),
"vecLocs[i] element count mismatch");
switch (normals.num_elements())
{
case 1:
// do nothing
......@@ -187,9 +182,9 @@ namespace Nektar
case 2:
{
const int vx = (int)velLoc[i*vecSize+0];
const int vy = (int)velLoc[i*vecSize+1];
const int nq = inarray[0].num_elements();
const int vx = (int)vecLocs[i][0];
const int vy = (int)vecLocs[i][1];
Vmath::Vmul (nq, inarray [vx], 1, normals [0], 1,
outarray[vx], 1);
......@@ -204,10 +199,10 @@ namespace Nektar
case 3:
{
const int vx = (int)velLoc[i*vecSize+0];
const int vy = (int)velLoc[i*vecSize+1];
const int vz = (int)velLoc[i*vecSize+2];
const int nq = inarray[0].num_elements();
const int vx = (int)vecLocs[i][0];
const int vy = (int)vecLocs[i][1];
const int vz = (int)vecLocs[i][2];
// Generate matrices if they don't already exist.
if (m_rotMat.num_elements() == 0)
......@@ -233,6 +228,7 @@ namespace Nektar
outarray[vz], 1, outarray[vz], 1);
break;
}
default:
ASSERTL0(false, "Invalid space dimension.");
break;
......@@ -251,25 +247,21 @@ namespace Nektar
void RiemannSolver::rotateFromNormal(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
const Array<OneD, const Array<OneD, NekDouble> > &normals,
const Array<OneD, const Array<OneD, NekDouble> > &vecLocs,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
ASSERTL0(CheckAuxiliary("velLoc"), "velLoc not defined.");
const Array<OneD, NekDouble> &velLoc = m_auxiliary["velLoc"]();
int vecCount = velLoc.num_elements() / normals.num_elements();
int vecSize = normals.num_elements();
ASSERTL0(velLoc.num_elements() % normals.num_elements() == 0,
"velLoc and normals mismatch");
for (int i = 0; i < inarray.num_elements(); ++i)
{
Vmath::Vcopy(inarray[i].num_elements(), inarray[i], 1,
outarray[i], 1);
}
for (int i = 0; i < vecCount; i++)
for (int i = 0; i < vecLocs.num_elements(); i++)
{
switch(vecSize)
ASSERTL0(vecLocs[i].num_elements() == normals.num_elements(),
"vecLocs[i] element count mismatch");
switch (normals.num_elements())
{
case 1:
// do nothing
......@@ -277,9 +269,9 @@ namespace Nektar
case 2:
{
const int vx = (int)velLoc[i*vecSize+0];
const int vy = (int)velLoc[i*vecSize+1];
const int nq = normals[0].num_elements();
const int vx = (int)vecLocs[i][0];
const int vy = (int)vecLocs[i][1];
Vmath::Vmul (nq, inarray [vy], 1, normals [1], 1,
outarray[vx], 1);
......@@ -291,12 +283,13 @@ namespace Nektar
outarray[vy], 1, outarray[vy], 1);
break;
}
case 3:
{
const int vx = (int)velLoc[i*vecSize+0];
const int vy = (int)velLoc[i*vecSize+1];
const int vz = (int)velLoc[i*vecSize+2];
const int nq = normals[0].num_elements();
const int vx = (int)vecLocs[i][0];
const int vy = (int)vecLocs[i][1];
const int vz = (int)vecLocs[i][2];
Vmath::Vvtvvtp(nq, inarray [vx], 1, m_rotMat[0], 1,
inarray [vy], 1, m_rotMat[3], 1,
......@@ -315,11 +308,11 @@ namespace Nektar
outarray[vz], 1, outarray[vz], 1);
break;
}
default:
ASSERTL0(false, "Invalid space dimension.");
break;
}
}
}
......@@ -363,16 +356,29 @@ namespace Nektar
}
/**
* @brief Determine whether a scalar has been defined in #m_auxiliary.
* @brief Determine whether a scalar has been defined in #m_auxScal.
*
* @param name Scalar name.
*/
bool RiemannSolver::CheckAuxiliary(std::string name)
bool RiemannSolver::CheckAuxScal(std::string name)
{
std::map<std::string, RSScalarFuncType>::iterator it =
m_auxiliary.find(name);
m_auxScal.find(name);
return it != m_auxScal.end();
}
/**
* @brief Determine whether a vector has been defined in #m_auxVec.
*
* @param name Vector name.
*/
bool RiemannSolver::CheckAuxVec(std::string name)
{
std::map<std::string, RSVecFuncType>::iterator it =
m_auxVec.find(name);
return it != m_auxiliary.end();
return it != m_auxVec.end();
}
/**
......@@ -395,7 +401,6 @@ namespace Nektar
{
m_rotMat[i] = Array<OneD, NekDouble>(nq);
}
for (i = 0; i < normals[0].num_elements(); ++i)
{
// Generate matrix which takes us from (1,0,0) vector to trace
......
......@@ -107,11 +107,19 @@ namespace Nektar
}
template<typename FuncPointerT, typename ObjectPointerT>
void SetAuxiliary(std::string name,
void SetAuxScal(std::string name,
FuncPointerT func,
ObjectPointerT obj)
{
m_auxiliary[name] = boost::bind(func, obj);
m_auxScal[name] = boost::bind(func, obj);
}
template<typename FuncPointerT, typename ObjectPointerT>
void SetAuxVec(std::string name,
FuncPointerT func,
ObjectPointerT obj)
{
m_auxVec[name] = boost::bind(func, obj);
}
std::map<std::string, RSScalarFuncType> &GetScalars()
......@@ -139,8 +147,10 @@ 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;
/// Map of auxiliary scalar function types.
std::map<std::string, RSScalarFuncType> m_auxScal;
/// Map of auxiliary vector function types.
std::map<std::string, RSVecFuncType> m_auxVec;
/// Rotation matrices for each trace quadrature point.
Array<OneD, Array<OneD, NekDouble> > m_rotMat;
/// Rotation storage
......@@ -162,15 +172,18 @@ namespace Nektar
void rotateToNormal (
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
const Array<OneD, const Array<OneD, NekDouble> > &normals,
const Array<OneD, const Array<OneD, NekDouble> > &vecLocs,
Array<OneD, Array<OneD, NekDouble> > &outarray);
void rotateFromNormal(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
const Array<OneD, const Array<OneD, NekDouble> > &normals,
const Array<OneD, const Array<OneD, NekDouble> > &vecLocs,
Array<OneD, Array<OneD, NekDouble> > &outarray);
bool CheckScalars (std::string name);
bool CheckVectors (std::string name);
bool CheckParams (std::string name);
bool CheckAuxiliary (std::string name);
bool CheckScalars (std::string name);
bool CheckVectors (std::string name);
bool CheckParams (std::string name);
bool CheckAuxScal (std::string name);
bool CheckAuxVec (std::string name);
};
/// A shared pointer to an EquationSystem object
......
......@@ -68,10 +68,11 @@ namespace Nektar
m_homoInitialFwd = false;
// Set up locations of velocity vector.
m_velLoc = Array<OneD, NekDouble>(m_spacedim);
m_vecLocs = Array<OneD, Array<OneD, NekDouble> >(1);
m_vecLocs[0] = Array<OneD, NekDouble>(m_spacedim);
for (int i = 0; i < m_spacedim; ++i)
{
m_velLoc[i] = i+1;
m_vecLocs[0][i] = 1 + i;
}
// Get gamma parameter from session file.
......@@ -180,19 +181,19 @@ namespace Nektar
// Setting up parameters for advection operator Riemann solver
m_riemannSolver->SetParam (
"gamma", &CompressibleFlowSystem::GetGamma, this);
m_riemannSolver->SetAuxiliary(
"velLoc", &CompressibleFlowSystem::GetVelLoc, this);
"gamma", &CompressibleFlowSystem::GetGamma, this);
m_riemannSolver->SetAuxVec(
"vecLocs", &CompressibleFlowSystem::GetVecLocs, this);
m_riemannSolver->SetVector(
"N", &CompressibleFlowSystem::GetNormals, this);
"N", &CompressibleFlowSystem::GetNormals, this);
// Setting up parameters for diffusion operator Riemann solver
m_riemannSolverLDG->SetParam (
"gamma", &CompressibleFlowSystem::GetGamma, this);
m_riemannSolverLDG->SetAuxiliary(
"velLoc", &CompressibleFlowSystem::GetVelLoc, this);
"gamma", &CompressibleFlowSystem::GetGamma, this);
m_riemannSolverLDG->SetVector(
"N", &CompressibleFlowSystem::GetNormals, this);
"vecLocs", &CompressibleFlowSystem::GetVecLocs, this);
m_riemannSolverLDG->SetVector(
"N", &CompressibleFlowSystem::GetNormals, this);
// Concluding initialisation of advection / diffusion operators
m_advection->SetRiemannSolver (m_riemannSolver);
......
......@@ -91,7 +91,7 @@ namespace Nektar
SolverUtils::RiemannSolverSharedPtr m_riemannSolverLDG;
SolverUtils::AdvectionSharedPtr m_advection;
SolverUtils::DiffusionSharedPtr m_diffusion;
Array<OneD, NekDouble> m_velLoc;
Array<OneD, Array<OneD, NekDouble> >m_vecLocs;
NekDouble m_gamma;
NekDouble m_pInf;
NekDouble m_rhoInf;
......@@ -196,9 +196,9 @@ namespace Nektar
return m_gamma;
}
const Array<OneD, NekDouble> &GetVelLoc()
const Array<OneD, const Array<OneD, NekDouble> > &GetVecLocs()
{
return m_velLoc;
return m_vecLocs;
}
const Array<OneD, const Array<OneD, NekDouble> > &GetNormals()
......
......@@ -115,9 +115,9 @@ namespace Nektar
m_riemannSolver->SetParam (
"gravity",
&LinearSWE::GetGravity, this);
m_riemannSolver->SetAuxiliary(
"velLoc",
&LinearSWE::GetVelLoc, this);
m_riemannSolver->SetAuxVec(
"vecLocs",
&LinearSWE::GetVecLocs, this);
m_riemannSolver->SetVector(
"N",
&LinearSWE::GetNormals, this);
......@@ -139,9 +139,9 @@ namespace Nektar
// m_riemannSolverLDG->AddParam (
// "gravity",
// &NonlinearSWE::GetGravity, this);
// m_riemannSolverLDG->AddAuxiliary(
// "velLoc",
// &NonlinearSWE::GetVelLoc, this);
// m_riemannSolverLDG->SetAuxVec(
// "vecLocs",
// &NonlinearSWE::GetVecLocs, this);
// m_riemannSolverLDG->AddVector(
// "N",
// &NonlinearSWE::GetNormals, this);
......
......@@ -111,9 +111,9 @@ namespace Nektar
m_riemannSolver->SetParam (
"gravity",
&NonlinearSWE::GetGravity, this);
m_riemannSolver->SetAuxiliary(
"velLoc",
&NonlinearSWE::GetVelLoc, this);
m_riemannSolver->SetAuxVec(
"vecLocs",
&NonlinearSWE::GetVecLocs, this);
m_riemannSolver->SetVector(
"N",
&NonlinearSWE::GetNormals, this);
......@@ -125,9 +125,9 @@ namespace Nektar
// m_riemannSolverLDG->AddParam (
// "gravity",
// &NonlinearSWE::GetGravity, this);
// m_riemannSolverLDG->AddAuxiliary(
// "velLoc",
// &NonlinearSWE::GetVelLoc, this);
// m_riemannSolverLDG->SetAuxVec(
// "vecLocs",
// &NonlinearSWE::GetVecLocs, this);
// m_riemannSolverLDG->AddVector(
// "N",
// &NonlinearSWE::GetNormals, this);
......
......@@ -82,11 +82,12 @@ namespace Nektar
}
// Set up locations of velocity vector.
m_velLoc = Array<OneD, NekDouble>(m_spacedim);
for (int i = 0; i < m_spacedim; ++i)
{
m_velLoc[i] = i+1;
}
m_vecLocs = Array<OneD, Array<OneD, NekDouble> >(1);
m_vecLocs[0] = Array<OneD, NekDouble>(m_spacedim);
for (int i = 0; i < m_spacedim; ++i)
{
m_vecLocs[0][i] = 1 + i;
}
// Load generic input parameters
m_session->LoadParameter("IO_InfoSteps", m_infosteps, 0);
......
......@@ -82,7 +82,7 @@ namespace Nektar
/// Coriolis force
Array<OneD, NekDouble> m_coriolis;
// Location of velocity vector.
Array<OneD, NekDouble> m_velLoc;
Array<OneD, Array<OneD, NekDouble> > m_vecLocs;
/// Initialises UnsteadySystem class members.
ShallowWaterSystem(const LibUtilities::SessionReaderSharedPtr& pSession);
......@@ -109,9 +109,9 @@ namespace Nektar
return m_g;
}
const Array<OneD, NekDouble> &GetVelLoc()
const Array<OneD, const Array<OneD, NekDouble> > &GetVecLocs()
{
return m_velLoc;
return m_vecLocs;
}
const Array<OneD, const Array<OneD, NekDouble> > &GetNormals()
......
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