Skip to content

GitLab

Nektar
Nektar
Commit 18ebba85 authored by Zhenguo Yan's avatar Zhenguo Yan
Browse files
Options

TMP

parent 34e8eb90
Hide whitespace changes
Inline Side-by-side
Showing with 1046 additions and 1318 deletions
library/SolverUtils/Diffusion/DiffusionLDG.cpp
View file @ 18ebba85
......@@ -32,611 +32,388 @@
//
///////////////////////////////////////////////////////////////////////////////
#include <SolverUtils/Diffusion/DiffusionLDG.h>
#include <iostream>
#include <iomanip>
#include <iostream>
#include <boost/algorithm/string/predicate.hpp>
#include <SolverUtils/Diffusion/DiffusionLDG.h>
namespace Nektar
{
namespace SolverUtils
{
std::string DiffusionLDG::type = GetDiffusionFactory().
RegisterCreatorFunction("LDG", DiffusionLDG::create);
namespace SolverUtils
{
std::string DiffusionLDG::type =
GetDiffusionFactory().RegisterCreatorFunction("LDG", DiffusionLDG::create);
DiffusionLDG::DiffusionLDG()
{
}
DiffusionLDG::DiffusionLDG()
{
}
void DiffusionLDG::v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
m_session = pSession;
void DiffusionLDG::v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
m_session = pSession;
m_session->LoadSolverInfo("ShockCaptureType",
m_shockCaptureType, "Off");
m_session->LoadSolverInfo("ShockCaptureType", m_shockCaptureType, "Off");
// Setting up the normals
int i;
int nDim = pFields[0]->GetCoordim(0);
int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
// Set up penalty term for LDG
m_session->LoadParameter("LDGc11", m_C11, 1.0);
m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(nDim);
for(i = 0; i < nDim; ++i)
{
m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
}
pFields[0]->GetTrace()->GetNormals(m_traceNormals);
}
// Setting up the normals
std::size_t nDim = pFields[0]->GetCoordim(0);
std::size_t nTracePts = pFields[0]->GetTrace()->GetTotPoints();
void DiffusionLDG::v_Diffuse(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd,
const Array<OneD, Array<OneD, NekDouble> > &pBwd)
{
int nCoeffs = fields[0]->GetNcoeffs();
Array<OneD, Array<OneD, NekDouble> > tmp(nConvectiveFields);
for(int i=0; i< nConvectiveFields; i++)
{
tmp[i] = Array<OneD, NekDouble>(nCoeffs,0.0);
}
m_traceNormals = Array<OneD, Array<OneD, NekDouble>>{nDim};
for (std::size_t i = 0; i < nDim; ++i)
{
m_traceNormals[i] = Array<OneD, NekDouble>{nTracePts};
}
pFields[0]->GetTrace()->GetNormals(m_traceNormals);
}
DiffusionLDG::v_Diffuse_coeff(nConvectiveFields,fields,inarray,tmp,pFwd,pBwd);
for (int i = 0; i < nConvectiveFields; ++i)
{
fields[i]->BwdTrans (tmp[i], outarray[i]);
}
}
void DiffusionLDG::v_Diffuse(
const std::size_t nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble>> &inarray,
Array<OneD, Array<OneD, NekDouble>> &outarray,
const Array<OneD, Array<OneD, NekDouble>> &pFwd,
const Array<OneD, Array<OneD, NekDouble>> &pBwd)
{
std::size_t nDim = fields[0]->GetCoordim(0);
std::size_t nPts = fields[0]->GetTotPoints();
std::size_t nCoeffs = fields[0]->GetNcoeffs();
std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
void DiffusionLDG::v_Diffuse_coeff(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
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);
int nPts = fields[0]->GetTotPoints();
int nCoeffs = fields[0]->GetNcoeffs();
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
Array<OneD, NekDouble> tmp{nCoeffs};
Array<OneD, Array<OneD, Array<OneD, NekDouble>>> flux{nDim};
Array<OneD, Array<OneD, Array<OneD, NekDouble>>> qfield{nDim};
Array<OneD, NekDouble> tmp(nCoeffs);
for (std::size_t j = 0; j < nDim; ++j)
{
qfield[j] = Array<OneD, Array<OneD, NekDouble>>{nConvectiveFields};
flux[j] = Array<OneD, Array<OneD, NekDouble>>{nConvectiveFields};
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > volumeFlux;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > qfield(nDim);
for (j = 0; j < nDim; ++j)
{
qfield[j] = Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
for (i = 0; i < nConvectiveFields; ++i)
{
qfield[j][i] = Array<OneD, NekDouble>(nPts, 0.0);
}
}
for (std::size_t i = 0; i < nConvectiveFields; ++i)
{
qfield[j][i] = Array<OneD, NekDouble>{nPts, 0.0};
flux[j][i] = Array<OneD, NekDouble>{nTracePts, 0.0};
}
}
Array<OneD, Array<OneD, NekDouble > > traceflux(nConvectiveFields);
for (int i = 0; i < nConvectiveFields; ++i)
{
traceflux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
}
// Compute q_{\eta} and q_{\xi}
// Obtain numerical fluxes
DiffuseCalculateDerivative(nConvectiveFields,fields,inarray,qfield,pFwd,pBwd);
DiffuseVolumeFlux(nConvectiveFields,fields,inarray,qfield,volumeFlux);
DiffuseTraceFlux(nConvectiveFields,fields,inarray,qfield,volumeFlux,traceflux,pFwd,pBwd);
NumFluxforScalar(fields, inarray, flux, pFwd, pBwd);
Array<OneD, Array<OneD, NekDouble> > qdbase(nDim);
for (std::size_t j = 0; j < nDim; ++j)
{
for (std::size_t i = 0; i < nConvectiveFields; ++i)
{
fields[i]->IProductWRTDerivBase(j, inarray[i], tmp);
Vmath::Neg(nCoeffs, tmp, 1);
fields[i]->AddTraceIntegral(flux[j][i], tmp);
fields[i]->SetPhysState(false);
fields[i]->MultiplyByElmtInvMass(tmp, tmp);
fields[i]->BwdTrans(tmp, qfield[j][i]);
}
}
for (i = 0; i < nConvectiveFields; ++i)
{
for (j = 0; j < nDim; ++j)
{
qdbase[j] = qfield[j][i];
}
fields[i]->IProductWRTDerivBase(qdbase,tmp);
// Evaulate <\phi, \hat{F}\cdot n> - outarray[i]
Vmath::Neg (nCoeffs, tmp, 1);
fields[i]->AddTraceIntegral (traceflux[i], tmp);
fields[i]->SetPhysState (false);
fields[i]->MultiplyByElmtInvMass(tmp, outarray[i]);
}
// Initialize viscous tensor
Array<OneD, Array<OneD, Array<OneD, NekDouble>>> viscTensor{nDim};
for (std::size_t j = 0; j < nDim; ++j)
{
viscTensor[j] = Array<OneD, Array<OneD, NekDouble>>{nConvectiveFields};
for (std::size_t i = 0; i < nConvectiveFields; ++i)
{
viscTensor[j][i] = Array<OneD, NekDouble>{nPts, 0.0};
}
}
// Get viscous tensor
m_fluxVector(inarray, qfield, viscTensor);
void DiffusionLDG::v_DiffuseCalculateDerivative(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD,Array<OneD, Array<OneD, NekDouble> > > &inarrayderivative,
const Array<OneD, Array<OneD, NekDouble> > &pFwd,
const Array<OneD, Array<OneD, NekDouble> > &pBwd)
// Compute u from q_{\eta} and q_{\xi}
// Obtain numerical fluxes
NumFluxforVector(fields, inarray, viscTensor, flux[0]);
Array<OneD, Array<OneD, NekDouble>> qdbase{nDim};
for (std::size_t i = 0; i < nConvectiveFields; ++i)
{
for (std::size_t j = 0; j < nDim; ++j)
{
int nDim = fields[0]->GetCoordim(0);
int nCoeffs = fields[0]->GetNcoeffs();
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
qdbase[j] = viscTensor[j][i];
}
fields[i]->IProductWRTDerivBase(qdbase, tmp);
// Evaulate <\phi, \hat{F}\cdot n> - outarray[i]
Vmath::Neg(nCoeffs, tmp, 1);
fields[i]->AddTraceIntegral(flux[0][i], tmp);
fields[i]->SetPhysState(false);
fields[i]->MultiplyByElmtInvMass(tmp, tmp);
fields[i]->BwdTrans(tmp, outarray[i]);
}
}
Array<OneD, NekDouble> tmp(nCoeffs);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > flux (nDim);
for (int j = 0; j < nDim; ++j)
{
flux[j] = Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
for (int i = 0; i < nConvectiveFields; ++i)
{
flux[j][i] = Array<OneD, NekDouble>(nTracePts, 0.0);
}
}
void DiffusionLDG::NumFluxforScalar(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble>> &ufield,
Array<OneD, Array<OneD, Array<OneD, NekDouble>>> &uflux,
const Array<OneD, Array<OneD, NekDouble>> &pFwd,
const Array<OneD, Array<OneD, NekDouble>> &pBwd)
{
std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
std::size_t nvariables = fields.num_elements();
std::size_t nDim = fields[0]->GetCoordim(0);
Array<OneD, NekDouble> Fwd{nTracePts};
Array<OneD, NekDouble> Bwd{nTracePts};
Array<OneD, NekDouble> fluxtemp{nTracePts, 0.0};
// Get the sign of (v \cdot n), v = an arbitrary vector
// Evaluate upwind flux:
// uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
for (std::size_t i = 0; i < nvariables; ++i)
{
// Compute Fwd and Bwd value of ufield of i direction
if (pFwd == NullNekDoubleArrayofArray ||
pBwd == NullNekDoubleArrayofArray)
{
fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
}
else
{
Fwd = pFwd[i];
Bwd = pBwd[i];
}
v_NumFluxforScalar(fields, inarray, flux, pFwd, pBwd);
// Upwind
Vmath::Vcopy(nTracePts, Fwd, 1, fluxtemp, 1);
for (int j = 0; j < nDim; ++j)
{
for (int i = 0; i < nConvectiveFields; ++i)
{
fields[i]->IProductWRTDerivBase(j, inarray[i], tmp);
Vmath::Neg (nCoeffs, tmp, 1);
fields[i]->AddTraceIntegral (flux[j][i], tmp);
fields[i]->SetPhysState (false);
fields[i]->MultiplyByElmtInvMass(tmp, tmp);
fields[i]->BwdTrans (tmp, inarrayderivative[j][i]);
}
}
// Imposing weak boundary condition with flux
if (fields[0]->GetBndCondExpansions().num_elements())
{
ApplyScalarBCs(fields, i, ufield[i], Fwd, Bwd, fluxtemp);
}
void DiffusionLDG::v_DiffuseVolumeFlux(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble>> &inarray,
Array<OneD,Array<OneD, Array<OneD, NekDouble> > > &inarrayderivative,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &VolumeFlux,
Array< OneD, int > &nonZeroIndex)
for (std::size_t j = 0; j < nDim; ++j)
{
if(VolumeFlux.num_elements())
{
int nDim = inarrayderivative.num_elements();
int nPts = inarray[nConvectiveFields-1].num_elements();
for (int j = 0; j < nDim; ++j)
{
for (int i = 0; i < nConvectiveFields; ++i)
{
Vmath::Vcopy(nPts,inarrayderivative[j][i],1,VolumeFlux[j][i],1);
}
}
}
else
{
VolumeFlux = inarrayderivative;
}
Vmath::Vmul(nTracePts, m_traceNormals[j], 1, fluxtemp, 1,
uflux[j][i], 1);
}
}
}
void DiffusionLDG::v_DiffuseTraceFlux(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble>> &inarray,
Array<OneD,Array<OneD, Array<OneD, NekDouble> > > &qfield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &VolumeFlux,
Array<OneD, Array<OneD, NekDouble> > &TraceFlux,
const Array<OneD, Array<OneD, NekDouble>> &pFwd,
const Array<OneD, Array<OneD, NekDouble>> &pBwd,
Array< OneD, int > &nonZeroIndex)
void DiffusionLDG::ApplyScalarBCs(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const std::size_t var, const Array<OneD, const NekDouble> &ufield,
const Array<OneD, const NekDouble> &Fwd,
const Array<OneD, const NekDouble> &Bwd,
Array<OneD, NekDouble> &penaltyflux)
{
boost::ignore_unused(ufield, Bwd);
// Number of boundary regions
std::size_t nBndRegions =
fields[var]->GetBndCondExpansions().num_elements();
std::size_t cnt = 0;
for (std::size_t i = 0; i < nBndRegions; ++i)
{
if (fields[var]->GetBndConditions()[i]->GetBoundaryConditionType() ==
SpatialDomains::ePeriodic)
{
int i,j,k,e;
int nBndEdgePts;
int nDim = qfield.num_elements();
int nPts = inarray[nConvectiveFields-1].num_elements();
int nTracePts = TraceFlux[nConvectiveFields-1].num_elements();
// Compute u from q_{\eta} and q_{\xi}
// Obtain numerical fluxes
v_NumFluxforVector(fields, inarray, qfield, TraceFlux);
if (m_ArtificialDiffusionVector)
{
Array<OneD, NekDouble> muvar(nPts, 0.0);
m_ArtificialDiffusionVector(inarray, muvar);
int numConvFields = nConvectiveFields;
if (m_shockCaptureType == "Smooth")
{
numConvFields = nConvectiveFields - 1;
}
for (j = 0; j < nDim; ++j)
{
for (i = 0; i < numConvFields; ++i)
{
Vmath::Vmul(nPts,VolumeFlux[j][i],1,muvar,1,VolumeFlux[j][i],1);
}
}
Array<OneD, NekDouble> FwdMuVar(nTracePts, 0.0);
Array<OneD, NekDouble> BwdMuVar(nTracePts, 0.0);
fields[0]->GetFwdBwdTracePhys(muvar,FwdMuVar,BwdMuVar);
int nBndRegions = fields[0]->GetBndCondExpansions().
num_elements();
int cnt = 0;
for (i = 0; i < nBndRegions; ++i)
{
if (fields[0]->GetBndConditions()[i]->GetBoundaryConditionType()
== SpatialDomains::ePeriodic)
{
continue;
}
// Number of boundary expansion related to that region
int nBndEdges = fields[0]->
GetBndCondExpansions()[i]->GetExpSize();
// Weakly impose boundary conditions by modifying flux
// values
for (e = 0; e < nBndEdges ; ++e)
{
nBndEdgePts = fields[0]->GetBndCondExpansions()[i]
->GetExp(e)->GetTotPoints();
int id2 = fields[0]->GetTrace()->GetPhys_Offset(
fields[0]->GetTraceMap()
->GetBndCondTraceToGlobalTraceMap(cnt++));
for (k = 0; k < nBndEdgePts; ++k)
{
BwdMuVar[id2+k] = 0.0;
}
}
}
for(i = 0; i < numConvFields; ++i)
{
for(k = 0; k < nTracePts; ++k)
{
TraceFlux[i][k] = 0.5 * (FwdMuVar[k] + BwdMuVar[k]) * TraceFlux[i][k];
}
}
}
continue;
}
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,
const Array<OneD, Array<OneD, NekDouble> > &pFwd,
const Array<OneD, Array<OneD, NekDouble> > &pBwd)
// Number of boundary expansion related to that region
std::size_t nBndEdges =
fields[var]->GetBndCondExpansions()[i]->GetExpSize();
// Weakly impose boundary conditions by modifying flux values
for (std::size_t e = 0; e < nBndEdges; ++e)
{
int i, j;
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
int nvariables = fields.num_elements();
int nDim = fields[0]->GetCoordim(0);;
Array<OneD, NekDouble > Fwd (nTracePts);
Array<OneD, NekDouble > Bwd (nTracePts);
Array<OneD, NekDouble > Vn (nTracePts, 0.0);
Array<OneD, NekDouble > fluxtemp(nTracePts, 0.0);
// Get the normal velocity Vn
for(i = 0; i < nDim; ++i)
std::size_t nBndEdgePts = fields[var]
->GetBndCondExpansions()[i]
->GetExp(e)
->GetTotPoints();
std::size_t id1 =
fields[var]->GetBndCondExpansions()[i]->GetPhys_Offset(e);
std::size_t id2 = fields[0]->GetTrace()->GetPhys_Offset(
fields[0]->GetTraceMap()->GetBndCondTraceToGlobalTraceMap(
cnt++));
// AV boundary conditions
if ( boost::iequals(fields[var]->GetBndConditions()[i]->
GetUserDefined(),"Wall") ||
boost::iequals(fields[var]->GetBndConditions()[i]->
GetUserDefined(),"Symmetry") ||
boost::iequals(fields[var]->GetBndConditions()[i]->
GetUserDefined(),"WallViscous") ||
boost::iequals(fields[var]->GetBndConditions()[i]->
GetUserDefined(),"WallAdiabatic"))
{
Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
Vn, 1, Vn, 1);
Vmath::Vcopy(nBndEdgePts, &Fwd[id2], 1, &penaltyflux[id2], 1);
}
// Get the sign of (v \cdot n), v = an arbitrary vector
// Evaluate upwind flux:
// uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
for (i = 0; i < nvariables ; ++i)
// For Dirichlet boundary condition: uflux = g_D
else if (fields[var]
->GetBndConditions()[i]
->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
{
// Compute Fwd and Bwd value of ufield of i direction
if (pFwd == NullNekDoubleArrayofArray ||
pBwd == NullNekDoubleArrayofArray)
{
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
// 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
Array<OneD, NekDouble > uplus(nTracePts);
fields[i]->ExtractTracePhys(ufield[i], uplus);
if(fields[0]->GetBndCondExpansions().num_elements())
{
v_WeakPenaltyforScalar(fields, i, ufield[i], uplus, 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})
// edge::eBackward, uFwd \(\tan_{\xi}^Bwd \cdot \vec{n})
// else if Vn < 0, flux = uBwd*(tan_{\xi}^- \cdot \vec{n}),
// i.e,
// edge::eForward, uBwd \(\tan_{\xi}^Fwd \cdot \vec{n})
// edge::eBackward, uBwd \(\tan_{\xi}^Bwd \cdot \vec{n})
Vmath::Vmul(nTracePts,
m_traceNormals[j], 1,
fluxtemp, 1,
uflux[j][i], 1);
}
Vmath::Vcopy(
nBndEdgePts,
&(fields[var]->GetBndCondExpansions()[i]->GetPhys())[id1],
1, &penaltyflux[id2], 1);
}
// For Neumann boundary condition: uflux = u+
else if ((fields[var]->GetBndConditions()[i])
->GetBoundaryConditionType() ==
SpatialDomains::eNeumann)
{
Vmath::Vcopy(nBndEdgePts, &Fwd[id2], 1, &penaltyflux[id2], 1);
}
}
}
}
/**
* @brief Build the numerical flux for the 2nd order derivatives
* todo: add variable coeff and h dependence to penalty term
*/
void DiffusionLDG::NumFluxforVector(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble>> &ufield,
Array<OneD, Array<OneD, Array<OneD, NekDouble>>> &qfield,
Array<OneD, Array<OneD, NekDouble>> &qflux)
{
std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
std::size_t nvariables = fields.num_elements();
std::size_t nDim = qfield.num_elements();
Array<OneD, NekDouble> Fwd{nTracePts};
Array<OneD, NekDouble> Bwd{nTracePts};
Array<OneD, NekDouble> qFwd{nTracePts};
Array<OneD, NekDouble> qBwd{nTracePts};
Array<OneD, NekDouble> qfluxtemp{nTracePts, 0.0};
Array<OneD, NekDouble> uterm{nTracePts};
// Evaulate upwind flux:
// qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
for (std::size_t i = 0; i < nvariables; ++i)
{