Commit 93fc0e14 authored by Dave Moxey's avatar Dave Moxey
Browse files

More changes to allow advection velocity on a per-plane basis.

parent d5596d85
......@@ -63,5 +63,26 @@ namespace Nektar
{
v_Advect(nConvectiveFields, fields, advVel, inarray, outarray);
}
void Advection::v_InitObject(
const LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
spaceDim = pFields[0]->GetCoordim(0);
if (pSession->->DefinesSolverInfo("HOMOGENEOUS"))
{
std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
if (HomoStr == "HOMOGENEOUS1D" || HomoStr == "Homogeneous1D" ||
HomoStr == "1D" || HomoStr == "Homo1D")
{
spaceDim++;
}
else
{
ASSERTL0(false, "Only 1D homogeneous dimension supported.");
}
}
}
}
}
......@@ -115,10 +115,7 @@ namespace Nektar
protected:
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
// Overridden by subclasses if necessary.
}
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_Advect(
const int nConvectiveFields,
......@@ -132,6 +129,8 @@ namespace Nektar
AdvectionFluxVecCB m_fluxVector;
/// Riemann solver for DG-type schemes.
RiemannSolverSharedPtr m_riemann;
/// Storage for space dimension. Used for homogeneous extension.
int m_spaceDim;
};
/// A shared pointer to an Advection object.
......
......@@ -35,16 +35,9 @@
#include <SolverUtils/Advection/Advection3DHomogeneous1D.h>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <LibUtilities/Foundations/Basis.h>
#include <LibUtilities/Foundations/Points.h>
#include <LibUtilities/Polylib/Polylib.h>
#include <StdRegions/StdSegExp.h>
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
#include <boost/math/special_functions/gamma.hpp>
#include <iostream>
#include <iomanip>
namespace Nektar
{
namespace SolverUtils
......@@ -79,8 +72,7 @@ namespace Nektar
{
// Strip trailing string "3DHomogeneous1D" to determine 2D advection
// type, and create an advection object for the plane.
string advName = advType.substr(
0, advType.length()-15);
string advName = advType.substr(0, advType.length()-15);
m_planeAdv = GetAdvectionFactory().CreateInstance(advName, advName);
}
......@@ -199,23 +191,21 @@ namespace Nektar
// this plane.
for (int j = 0; j < nConvectiveFields; ++j)
{
m_fieldsPlane[j] = fields[j]->GetPlane(i);
m_inarrayPlane[j] = Array<OneD, NekDouble>(
m_numPointsPlane, tmp2 = inarray[j] + m_planePos[i]);
m_fieldsPlane [j] = fields[j]->GetPlane(i);
m_inarrayPlane [j] = Array<OneD, NekDouble>(
m_numPointsPlane, tmp2 = inarray [j] + m_planePos[i]);
m_outarrayPlane[j] = Array<OneD, NekDouble>(
m_numPointsPlane, tmp2 = outarray[j] + m_planePos[i]);
}
/*
for (int j = 0; j < nVel; ++j)
{
if (advVel[j].num_elements() != 0)
{
m_advVelPlane[j] = Array<OneD, NekDouble>(
nPointsTot_plane, tmp2 = advVel[j] + i*nPointsTot_plane);
m_numPointsPlane, tmp2 = advVel[j] + m_planePos[i]);
}
}
*/
// Compute advection term for this plane.
m_planeAdv->Advect(nConvectiveFields, m_fieldsPlane,
......
......@@ -89,32 +89,6 @@ namespace Nektar
{
v_SetupMetrics(pSession, pFields);
v_SetupCFunctions(pSession, pFields);
// Initialising the fluxvector
int i, j;
int nConvectiveFields = pFields.num_elements();
int nDimensions = pFields[0]->GetCoordim(0);
int nSolutionPts = pFields[0]->GetTotPoints();
int spaceDim;
spaceDim = nDimensions;
if (pSession->DefinesSolverInfo("HOMOGENEOUS"))
{
spaceDim = 3;
}
m_fluxvector = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >(
nConvectiveFields);
for (i = 0; i < nConvectiveFields; ++i)
{
m_fluxvector[i] = Array<OneD, Array<OneD, NekDouble> >(
spaceDim);
for (j = 0; j < spaceDim; ++j)
{
m_fluxvector[i][j] = Array<OneD, NekDouble>(nSolutionPts);
}
}
}
/**
......@@ -810,7 +784,7 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
int i, n;
int i, j, n;
int phys_offset;
Array<OneD, NekDouble> auxArray1, auxArray2;
......@@ -824,14 +798,29 @@ namespace Nektar
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
int nCoeffs = fields[0]->GetNcoeffs();
// Storage for flux vector.
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > fluxvector
(nConvectiveFields);
// Outarray for Galerkin projection in case of primitive dealising
Array<OneD, Array<OneD, NekDouble> > outarrayCoeff(nConvectiveFields);
Array<OneD, Array<OneD, NekDouble> > outarrayCoeff
(nConvectiveFields);
// Store forwards/backwards space along trace space
Array<OneD, Array<OneD, NekDouble> > Fwd (nConvectiveFields);
Array<OneD, Array<OneD, NekDouble> > Bwd (nConvectiveFields);
Array<OneD, Array<OneD, NekDouble> > numflux(nConvectiveFields);
// Set up storage for flux vector.
for (i = 0; i < nConvectiveFields; ++i)
{
fluxvector[i] =
Array<OneD, Array<OneD, NekDouble> >(m_spaceDim);
for (j = 0; j < nVel; ++j)
{
fluxvector[i][j] = Array<OneD, NekDouble>(nSolutionPts);
}
}
for (i = 0; i < nConvectiveFields; ++i)
{
outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
......@@ -854,7 +843,7 @@ namespace Nektar
Array<OneD, NekDouble> divFC (nSolutionPts);
// Get the advection flux vector (solver specific)
m_fluxVector(inarray, m_fluxvector);
m_fluxVector(inarray, fluxvector);
// Get the discontinuous flux divergence
for(i = 0; i < nConvectiveFields; ++i)
......@@ -864,13 +853,13 @@ namespace Nektar
phys_offset = fields[0]->GetPhys_Offset(n);
fields[i]->GetExp(n)->PhysDeriv(
0, m_fluxvector[i][0] + phys_offset,
0, fluxvector[i][0] + phys_offset,
auxArray2 = DfluxvectorX1 + phys_offset);
}
// Get the correction flux divergence
v_DivCFlux_1D(nConvectiveFields, fields,
m_fluxvector[i][0], numflux[i], divFC);
fluxvector[i][0], numflux[i], divFC);
// Back to the physical space using global operations
Vmath::Vdiv(nSolutionPts, &divFC[0], 1, &m_jac[0], 1,
......@@ -898,7 +887,7 @@ namespace Nektar
Array<OneD, NekDouble> divFC(nSolutionPts);
// Get the advection flux vector (solver specific)
m_fluxVector(inarray, m_fluxvector);
m_fluxVector(inarray, fluxvector);
for (i = 0; i < nConvectiveFields; ++i)
{
......@@ -908,9 +897,9 @@ namespace Nektar
Vmath::Vvtvvtp(nSolutionPts,
&m_gmat[0][0], 1,
&m_fluxvector[i][0][0], 1,
&fluxvector[i][0][0], 1,
&m_gmat[2][0], 1,
&m_fluxvector[i][1][0], 1,
&fluxvector[i][1][0], 1,
&f_hat[0], 1);
Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &f_hat[0], 1,
......@@ -918,9 +907,9 @@ namespace Nektar
Vmath::Vvtvvtp(nSolutionPts,
&m_gmat[1][0], 1,
&m_fluxvector[i][0][0], 1,
&fluxvector[i][0][0], 1,
&m_gmat[3][0], 1,
&m_fluxvector[i][1][0], 1,
&fluxvector[i][1][0], 1,
&g_hat[0], 1);
Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &g_hat[0], 1,
......@@ -956,7 +945,7 @@ namespace Nektar
{
v_DivCFlux_2D(
nConvectiveFields, fields,
m_fluxvector[i][0], m_fluxvector[i][1],
fluxvector[i][0], fluxvector[i][1],
numflux[i], divFC);
}
......@@ -1004,7 +993,7 @@ namespace Nektar
const Array<OneD, const NekDouble> &numericalFlux,
Array<OneD, NekDouble> &divCFlux)
{
int i, n;
int n;
int nLocalSolutionPts, phys_offset;
LibUtilities::BasisSharedPtr Basis;
......@@ -1271,8 +1260,7 @@ namespace Nektar
{
int n, e, i, j, cnt;
int nElements = fields[0]->GetExpSize();
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
int nElements = fields[0]->GetExpSize();
int nLocalSolutionPts;
int nEdgePts;
......
......@@ -75,7 +75,6 @@ namespace Nektar
AdvectionFR(std::string advType);
Array<OneD, Array<OneD, NekDouble> > m_traceNormals;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > m_fluxvector;
std::string m_advType;
......@@ -129,12 +128,6 @@ namespace Nektar
const Array<OneD, const NekDouble> &fluxX3,
const Array<OneD, const NekDouble> &numericalFlux,
Array<OneD, NekDouble> &divCFlux);
virtual void v_SetFluxVector(
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &fluxvector)
{
fluxvector = m_fluxvector;
}
};
}
}
......
......@@ -49,8 +49,8 @@ namespace Nektar
}
/**
* @brief Initiliase AdvectionWeakDG objects and store them before starting
* the time-stepping.
* @brief Initialise AdvectionWeakDG objects and store them before
* starting the time-stepping.
*
* @param pSession Pointer to session reader.
* @param pFields Pointer to fields.
......@@ -97,8 +97,7 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> >(nVel);
for (j = 0; j < nVel; ++j)
{
fluxvector[i][j] =
Array<OneD, NekDouble>(nPointsTot);
fluxvector[i][j] = Array<OneD, NekDouble>(nPointsTot);
}
}
......
......@@ -331,47 +331,14 @@ namespace Nektar
int i , j;
int nq = physfield[0].num_elements();
// For 3DHomogenoeus1D
/*
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 // For general case
for (i = 0; i < flux.num_elements(); ++i)
{
*/
for (i = 0; i < flux.num_elements(); ++i)
for (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,
Vmath::Vmul(nq, physfield[i], 1, m_velocity[j], 1,
flux[i][j], 1);
}
}
//}
}
}
/**
......
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