Commit 316d0aac authored by Douglas Serson's avatar Douglas Serson
Browse files

Allow passing Fwd and Bwd arrays to advection and diffusion

parent 21511805
......@@ -79,9 +79,12 @@ void Advection::Advect(
const Array<OneD, Array<OneD, NekDouble> > &pAdvVel,
const Array<OneD, Array<OneD, NekDouble> > &pInarray,
Array<OneD, Array<OneD, NekDouble> > &pOutarray,
const NekDouble &pTime)
const NekDouble &pTime,
const Array<OneD, Array<OneD, NekDouble> > &pFwd,
const Array<OneD, Array<OneD, NekDouble> > &pBwd)
{
v_Advect(nConvectiveFields, pFields, pAdvVel, pInarray, pOutarray, pTime);
v_Advect(nConvectiveFields, pFields, pAdvVel, pInarray,
pOutarray, pTime, pFwd, pBwd);
}
......
......@@ -81,7 +81,9 @@ public:
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time);
const NekDouble &time,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
/**
* @brief Set the flux vector callback function.
......@@ -147,7 +149,9 @@ protected:
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)=0;
const NekDouble &time,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray)=0;
/// Overrides the base flow used during linearised advection
SOLVER_UTILS_EXPORT virtual void v_SetBaseFlow(
......
......@@ -194,7 +194,9 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)
const NekDouble &time,
const Array<OneD, Array<OneD, NekDouble> > &pFwd,
const Array<OneD, Array<OneD, NekDouble> > &pBwd)
{
Array<OneD, NekDouble> tmp(m_numPoints), tmp2;
int nVel = advVel.num_elements();
......
......@@ -84,7 +84,9 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time);
const NekDouble &time,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
private:
void ModifiedFluxVector(
......
......@@ -804,7 +804,9 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)
const NekDouble &time,
const Array<OneD, Array<OneD, NekDouble> > &pFwd,
const Array<OneD, Array<OneD, NekDouble> > &pBwd)
{
int i, j, n;
int phys_offset;
......
......@@ -88,7 +88,9 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time);
const NekDouble &time,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
virtual void v_SetupMetrics(
LibUtilities::SessionReaderSharedPtr pSession,
......
......@@ -68,7 +68,9 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)
const NekDouble &time,
const Array<OneD, Array<OneD, NekDouble> > &pFwd,
const Array<OneD, Array<OneD, NekDouble> > &pBwd)
{
int nDim = advVel.num_elements();
int nPointsTot = fields[0]->GetNpoints();
......
......@@ -65,7 +65,9 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time);
const NekDouble &time,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
};
}
}
......
......@@ -79,7 +79,9 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)
const NekDouble &time,
const Array<OneD, Array<OneD, NekDouble> > &pFwd,
const Array<OneD, Array<OneD, NekDouble> > &pBwd)
{
int nPointsTot = fields[0]->GetTotPoints();
int nCoeffs = fields[0]->GetNcoeffs();
......@@ -119,12 +121,24 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > Bwd (nConvectiveFields);
Array<OneD, Array<OneD, NekDouble> > numflux(nConvectiveFields);
for(i = 0; i < nConvectiveFields; ++i)
if (pFwd == NullNekDoubleArrayofArray ||
pBwd == NullNekDoubleArrayofArray)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
Bwd[i] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
numflux[i] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
for(i = 0; i < nConvectiveFields; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
Bwd[i] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
numflux[i] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
}
}
else
{
for(i = 0; i < nConvectiveFields; ++i)
{
Fwd[i] = pFwd[i];
Bwd[i] = pBwd[i];
}
}
m_riemann->Solve(m_spaceDim, Fwd, Bwd, numflux);
......
......@@ -65,7 +65,9 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time);
const NekDouble &time,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
};
}
}
......
......@@ -59,9 +59,11 @@ namespace Nektar
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd,
const Array<OneD, Array<OneD, NekDouble> > &pBwd)
{
v_Diffuse(nConvectiveFields, fields, inarray, outarray);
v_Diffuse(nConvectiveFields, fields, inarray, outarray, pFwd, pBwd);
}
}
}
......@@ -79,7 +79,9 @@ namespace Nektar
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
SOLVER_UTILS_EXPORT void FluxVec(
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
......@@ -145,7 +147,9 @@ namespace Nektar
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)=0;
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray)=0;
virtual void v_SetHomoDerivs(
Array<OneD, Array<OneD, NekDouble> > &deriv)
......
......@@ -190,7 +190,9 @@ namespace Nektar
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd,
const Array<OneD, Array<OneD, NekDouble> > &pBwd)
{
Array<OneD, NekDouble> tmp(m_numPoints), tmp2;
......
......@@ -80,7 +80,9 @@ namespace Nektar
const int nConvective,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
};
}
}
......
......@@ -74,7 +74,9 @@ namespace Nektar
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd,
const Array<OneD, Array<OneD, NekDouble> > &pBwd)
{
int nBndEdgePts, i, j, k, e;
int nDim = fields[0]->GetCoordim(0);
......@@ -115,7 +117,7 @@ namespace Nektar
// Compute q_{\eta} and q_{\xi}
// Obtain numerical fluxes
v_NumFluxforScalar(fields, inarray, flux);
v_NumFluxforScalar(fields, inarray, flux, pFwd, pBwd);
for (j = 0; j < nDim; ++j)
{
......@@ -219,7 +221,9 @@ namespace Nektar
void DiffusionLDG::v_NumFluxforScalar(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &ufield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &uflux)
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &uflux,
const Array<OneD, Array<OneD, NekDouble> > &pFwd,
const Array<OneD, Array<OneD, NekDouble> > &pBwd)
{
int i, j;
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
......@@ -241,38 +245,45 @@ namespace Nektar
// Get the sign of (v \cdot n), v = an arbitrary vector
// Evaluate upwind flux:
// uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
for (j = 0; j < nDim; ++j)
for (i = 0; i < nvariables ; ++i)
{
for (i = 0; i < nvariables ; ++i)
// Compute Fwd and Bwd value of ufield of i direction
if (pFwd == NullNekDoubleArrayofArray ||
pBwd == NullNekDoubleArrayofArray)
{
// Compute Fwd and Bwd value of ufield of i direction
fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
}
else
{
Fwd = pFwd[i];
Bwd = pBwd[i];
}
// if Vn >= 0, flux = uFwd, i.e.,
// edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
// edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
// if Vn >= 0, flux = uFwd, i.e.,
// edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
// edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
// else if Vn < 0, flux = uBwd, i.e.,
// edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
// edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
fields[i]->GetTrace()->Upwind(/*m_traceNormals[j]*/Vn,
// else if Vn < 0, flux = uBwd, i.e.,
// edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
// edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
fields[i]->GetTrace()->Upwind(/*m_traceNormals[j]*/Vn,
Fwd, Bwd, fluxtemp);
// Imposing weak boundary condition with flux
// if Vn >= 0, uflux = uBwd at Neumann, i.e.,
// edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
// edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
// if Vn >= 0, uflux = uFwd at Neumann, i.e.,
// edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
// edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
// Imposing weak boundary condition with flux
// if Vn >= 0, uflux = uBwd at Neumann, i.e.,
// edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
// edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
if(fields[0]->GetBndCondExpansions().num_elements())
{
v_WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
}
// if Vn >= 0, uflux = uFwd at Neumann, i.e.,
// edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
// edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
if(fields[0]->GetBndCondExpansions().num_elements())
{
v_WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
}
for (j = 0; j < nDim; ++j)
{
// if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}),
// i.e,
// edge::eForward, uFwd \(\tan_{\xi}^Fwd \cdot \vec{n})
......
......@@ -68,12 +68,16 @@ namespace Nektar
const int nConvective,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
virtual void v_NumFluxforScalar(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &ufield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&uflux);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&uflux,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
virtual void v_WeakPenaltyforScalar(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
......
......@@ -106,7 +106,9 @@ namespace Nektar
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd,
const Array<OneD, Array<OneD, NekDouble> > &pBwd)
{
int i, j;
int nDim = fields[0]->GetCoordim(0);
......@@ -139,7 +141,7 @@ namespace Nektar
}
// Compute the numerical fluxes for the first order derivatives
v_NumericalFluxO1(fields, inarray, numericalFluxO1);
v_NumericalFluxO1(fields, inarray, numericalFluxO1, pFwd, pBwd);
for (j = 0; j < nDim; ++j)
{
......@@ -217,7 +219,9 @@ namespace Nektar
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
&numericalFluxO1)
&numericalFluxO1,
const Array<OneD, Array<OneD, NekDouble> > &pFwd,
const Array<OneD, Array<OneD, NekDouble> > &pBwd)
{
int i, j;
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
......@@ -234,17 +238,26 @@ namespace Nektar
}
// Store forwards/backwards space along trace space
Array<OneD, Array<OneD, NekDouble> > Fwd (nScalars);
Array<OneD, Array<OneD, NekDouble> > Bwd (nScalars);
Array<OneD, NekDouble> Fwd;
Array<OneD, NekDouble> Bwd;
Array<OneD, Array<OneD, NekDouble> > numflux(nScalars);
for (i = 0; i < nScalars; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
Bwd[i] = Array<OneD, NekDouble>(nTracePts);
if (pFwd == NullNekDoubleArrayofArray ||
pBwd == NullNekDoubleArrayofArray)
{
Fwd = Array<OneD, NekDouble>(nTracePts);
Bwd = Array<OneD, NekDouble>(nTracePts);
fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd, Bwd);
}
else
{
Fwd = pFwd[i];
Bwd = pBwd[i];
}
numflux[i] = Array<OneD, NekDouble>(nTracePts);
fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
fields[0]->GetTrace()->Upwind(Vn, Fwd[i], Bwd[i], numflux[i]);
fields[0]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux[i]);
}
// Extract internal values of the scalar variables for Neumann bcs
......
......@@ -82,13 +82,17 @@ namespace Nektar
const int nConvective,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
virtual void v_NumericalFluxO1(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
&numericalFluxO1);
&numericalFluxO1,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
virtual void v_WeakPenaltyO1(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
......
......@@ -850,7 +850,9 @@ namespace Nektar
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd,
const Array<OneD, Array<OneD, NekDouble> > &pBwd)
{
int i, j, n;
int phys_offset;
......
......@@ -107,7 +107,9 @@ namespace Nektar
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
virtual void v_NumFluxforScalar(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
......
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