Commit 7598e614 authored by Dave Moxey's avatar Dave Moxey

Merge branch 'fix/LDG-penalty' into 'master'

Fix/ldg penalty

See merge request !1080
parents 36b04d7e 6cb5ebd5
Pipeline #1181 failed with stages
in 7 minutes and 43 seconds
......@@ -14,3 +14,6 @@ ThirdParty
*.pdf
# Kdevelop project files
*.kdev4
# Sublime project files
*.sublime-project
*.sublime-workspace
......@@ -132,6 +132,7 @@ v5.0.0
- Fix surface extraction, added regression test (!994)
- Fix 2D meshing running out of memory due to missing else (!1012)
- Add support for .msh v4.1 file input (!1054)
- Added penalty term to LDG and LDGNS, slight generalization of LDG (!1080)
**FieldConvert**:
- Add input module for Semtex field files (!777)
......
......@@ -11,14 +11,14 @@ The ADRSolver is designed to solve partial differential equations of the form:
\begin{equation}
\alpha \dfrac{\partial u}{\partial t} + \lambda u + \nu \nabla u + \epsilon \nabla \cdot (D \nabla u) = f
\end{equation}
in either discontinuous or continuous projections of the solution field.
For a full list of the equations which are supported, and the capabilities of each equation,
in either discontinuous or continuous projections of the solution field.
For a full list of the equations which are supported, and the capabilities of each equation,
see the table below.
\begin{table}[h!]
\begin{center}
\tiny
\renewcommand\arraystretch{2.2}
\renewcommand\arraystretch{2.2}
\begin{tabular}{llll}
\toprule
\textbf{Equation to solve} & \textbf{EquationType} & \textbf{Dimensions} &
......@@ -26,22 +26,22 @@ see the table below.
\midrule
$u = f$ &
\inltt{Projection} & All & Continuous/Discontinuous \\
$\nabla^2 u = 0$ &
$\nabla^2 u = 0$ &
\inltt{Laplace} & All & Continuous/Discontinuous \\
$\nabla^2 u = f$ &
$\nabla^2 u = f$ &
\inltt{Poisson} & All & Continuous/Discontinuous \\
$\nabla^2 u + \lambda u = f$ &
$\nabla^2 u + \lambda u = f$ &
\inltt{Helmholtz} & All & Continuous/Discontinuous \\
$\epsilon \nabla^2 u + \mathbf{V}\nabla u = f$ &
$\epsilon \nabla^2 u + \mathbf{V}\nabla u = f$ &
\inltt{SteadyAdvectionDiffusion} & 2D only & Continuous/Discontinuous \\
$\epsilon \nabla^2 u + \lambda u = f$ &
$\epsilon \nabla^2 u + \lambda u = f$ &
\inltt{SteadyDiffusionReaction} & 2D only & Continuous/Discontinuous \\
$\epsilon \nabla^2 u + \mathbf{V}\nabla u + \lambda u = f$ &
\inltt{SteadyAdvectionDiffusionReaction} & 2D only &
$\epsilon \nabla^2 u + \mathbf{V}\nabla u + \lambda u = f$ &
\inltt{SteadyAdvectionDiffusionReaction} & 2D only &
Continuous/Discontinuous \\
$ \dfrac{\partial u}{\partial t} + \mathbf{V}\nabla u = f$ &
\inltt{UnsteadyAdvection} & All & Continuous/Discontinuous \\
$\dfrac{\partial u}{\partial t} = \epsilon \nabla^2 u$ &
$\dfrac{\partial u}{\partial t} = \epsilon \nabla^2 u$ &
\inltt{UnsteadyDiffusion} & All & Continuous/Discontinuous \\
$\dfrac{\partial u}{\partial t} = \epsilon \nabla^2 u + R(u)$ &
\inltt{UnsteadyReactionDiffusion} & All & Continuous \\
......@@ -64,7 +64,7 @@ ADRSolver session.xml
\section{Session file configuration}
The type of equation which is to be solved is specified through the EquationType
The type of equation which is to be solved is specified through the EquationType
SOLVERINFO option in the session file. This can be set as in table \ref{t:ADR1}.
At present, the Steady non-symmetric solvers cannot be used in parallel. \\
......@@ -76,7 +76,7 @@ The solver info are listed below:
\item \textbf{TimeIntegrationMethod}: The following types of time integration methods have been tested with each solver:
\begin{center}
\footnotesize
\renewcommand\arraystretch{1.2}
\renewcommand\arraystretch{1.2}
\begin{tabular}{lcccc}
\toprule
\textbf{EqType} & \textbf{Explicit} & \textbf{Diagonally Implicit} &
......@@ -116,7 +116,7 @@ The solver info are listed below:
\end{itemize}
\item \textbf{DiffusionType}:
\begin{itemize}
\item \inltt{LDG}.
\item \inltt{LDG} (The penalty term is proportional to an optional parameter \inltt{LDGc11} which is by default set to one; proportionality to polynomial order can be manually imposed by setting the parameter \inltt{LDGc11} equal to $p^2$).
\end{itemize}
\item \textbf{UpwindType}:
\begin{itemize}
......@@ -129,13 +129,13 @@ The solver info are listed below:
The following parameters can be specified in the \inltt{PARAMETERS} section of
the session file:
\begin{itemize}
\item \inltt{epsilon}: sets the diffusion coefficient $\epsilon$.\\
\item \inltt{epsilon}: sets the diffusion coefficient $\epsilon$.\\
\textit{Can be used} in: SteadyDiffusionReaction, SteadyAdvectionDiffusionReaction, UnsteadyDiffusion, UnsteadyAdvectionDiffusion. \\
\textit{Default value}: 0.
\item \inltt{d00}, \inltt{d11}, \inltt{d22}: sets the diagonal entries of the
diffusion tensor $D$. \\
\textit{Can be used in}: UnsteadyDiffusion \\
\textit{Default value}: All set to 1 (i.e. identity matrix).
\textit{Default value}: All set to 1 (i.e. identity matrix).
\item \inltt{lambda}: sets the reaction coefficient $\lambda$. \\
\textit{Can be used in}: SteadyDiffusionReaction, Helmholtz, SteadyAdvectionDiffusionReaction\\
\textit{Default value}: 0.
......@@ -215,7 +215,7 @@ and then define the actual advection function as
Two boundary regions are defined, one at each end of the domain, and periodicity
is enforced
\begin{lstlisting}[style=XMLStyle]
\begin{lstlisting}[style=XMLStyle]
<BOUNDARYREGIONS>
<B ID="0"> C[1] </B>
<B ID="1"> C[2] </B>
......@@ -268,7 +268,7 @@ We consider the elliptic partial differential equation:
where $\nabla^2$ is the Laplacian and $\lambda$ is a real positive constant.
\subsubsection{Input file}
\subsubsection{Input file}
The input for this example is given in the example file
\inlsh{Helmholtz2D\_modal.xml}
......@@ -300,7 +300,7 @@ basis with $7$ modes (maximum polynomial order is $6$).
\end{lstlisting}
Only one parameter is needed for this problem. In this example $\lambda = 1$ and
the Continuous Galerkin Method is used as projection scheme to solve the
the Continuous Galerkin Method is used as projection scheme to solve the
Helmholtz equation, so we need to specify the following parameters and solver
information.
\begin{lstlisting}[style=XMLStyle]
......@@ -331,7 +331,7 @@ region is then assigned an appropriate boundary condition.
<D VAR="u" VALUE="sin(PI*x)*sin(PI*y)" />
</REGION>
<REGION REF="1">
<R VAR="u" VALUE="sin(PI*x)*sin(PI*y)-PI*sin(PI*x)*cos(PI*y)"
<R VAR="u" VALUE="sin(PI*x)*sin(PI*y)-PI*sin(PI*x)*cos(PI*y)"
PRIMCOEFF="1" />
</REGION>
<REGION REF="2">
......@@ -343,7 +343,7 @@ region is then assigned an appropriate boundary condition.
</BOUNDARYCONDITIONS>
\end{lstlisting}
We know that for $f = -(\lambda + 2 \pi^2)sin(\pi x)cos(\pi y)$, the exact
We know that for $f = -(\lambda + 2 \pi^2)sin(\pi x)cos(\pi y)$, the exact
solution of the two-dimensional Helmholtz equation is $u = sin(\pi x)cos(\pi
y)$. These functions are defined specified to initialise the problem and verify
the correct solution is obtained by evaluating the $L_2$ and $L_{inf}$ errors.
......@@ -363,7 +363,7 @@ the correct solution is obtained by evaluating the $L_2$ and $L_{inf}$ errors.
ADRSolver Test_Helmholtz2D_modal.xml
\end{lstlisting}
This execution should print out a summary of input file, the $L_2$ and
This execution should print out a summary of input file, the $L_2$ and
$L_{inf}$ errors and the time spent on the calculation.
\subsubsection{Post-processing}
......@@ -417,9 +417,9 @@ Schmidt number (the relative thickness of the momentum to mass transfer boundary
layer) is sufficiently large. The analytical solution for the non-dimensional
mass transfer at the wall is given by:
\begin{align*}
S h(z) = \dfrac{2^{4/3}(Pe R/z)^{1/3}}{g^{1/3}\Gamma(4/3)} ,
S h(z) = \dfrac{2^{4/3}(Pe R/z)^{1/3}}{g^{1/3}\Gamma(4/3)} ,
\end{align*}
where $z$ is the streamwise coordinate, $R$ the pipe radius, $\Gamma(4/3)$ an incomplete
where $z$ is the streamwise coordinate, $R$ the pipe radius, $\Gamma(4/3)$ an incomplete
Gamma function and $Pe$ the Peclet number given by:
\begin{align*}
Pe = \dfrac{2 U R}{\epsilon}
......@@ -538,7 +538,7 @@ flow), hence we can define this through an analytical function as follows:
\end{lstlisting}
We assume that the initial domain concentration is uniform everywhere and the
same as the inlet. This is defined by,
same as the inlet. This is defined by,
\begin{lstlisting}[style=XMLStyle]
<FUNCTION NAME="InitialConditions">
<E VAR="u" VALUE="1" />
......
......@@ -272,7 +272,7 @@ Note that only \inltt{WeakDG} is fully supported, the other operators work only
\item \inltt{DiffusionType} is the diffusion operator we want to use
for the Navier-Stokes equations:
\begin{itemize}
\item \inltt{LDGNS} (LDG);
\item \inltt{LDGNS} (LDG with primitive variables. The penalty term is inversely proportional to the element size, proportional to the local viscosity for the momentum equations and to the thermal conductivity for the energy equation, and proportional to an optional parameter \inltt{LDGNSc11} which is by default set to one; proportionality to polynomial order can be manually imposed by setting the parameter \inltt{LDGNSc11} equal to $p^2$);
\item \inltt{LFRDGNS} (Flux-Reconstruction recovering nodal DG scheme);
\item \inltt{LFRSDNS} (Flux-Reconstruction recovering a spectral difference (SD) scheme);
\item \inltt{LFRHUNS} (Flux-Reconstruction recovering Huynh (G2) scheme);
......
......@@ -43,16 +43,16 @@ namespace Nektar
static DiffusionFactory instance;
return instance;
}
void Diffusion::InitObject(
const LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
v_InitObject(pSession, pFields);
}
void Diffusion::Diffuse(
const int nConvectiveFields,
const std::size_t nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
......
......@@ -52,23 +52,28 @@ namespace Nektar
namespace SolverUtils
{
typedef std::function<void (
const int,
const int,
const Array<OneD, Array<OneD, NekDouble> >&,
Array<OneD, Array<OneD, NekDouble> >&,
Array<OneD, Array<OneD, NekDouble> >&)> DiffusionFluxVecCB;
const Array<OneD, Array<OneD, NekDouble> > &,
const Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&)>
DiffusionFluxVecCB;
typedef std::function<void (
const Array<OneD, Array<OneD, NekDouble> >&,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&)>
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&)>
DiffusionFluxVecCBNS;
typedef std::function<void (
const Array<OneD, Array<OneD, NekDouble> >&,
const Array<OneD, Array<OneD, NekDouble> >&,
Array<OneD, Array<OneD, NekDouble> >&)>
DiffusionFluxPenaltyNS;
typedef std::function<void (
const Array<OneD, Array<OneD, NekDouble> >&,
Array<OneD, NekDouble >&)>
DiffusionArtificialDiffusion;
class Diffusion
{
public:
......@@ -79,34 +84,33 @@ namespace Nektar
SOLVER_UTILS_EXPORT void InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
SOLVER_UTILS_EXPORT void Diffuse(
const int nConvectiveFields,
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 = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
SOLVER_UTILS_EXPORT void FluxVec(
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
&fluxvector);
template<typename FuncPointerT, typename ObjectPointerT>
template<typename FuncPointerT, typename ObjectPointerT>
void SetFluxVector(FuncPointerT func, ObjectPointerT obj)
{
m_fluxVector = std::bind(
func, obj, std::placeholders::_1, std::placeholders::_2,
std::placeholders::_3, std::placeholders::_4,
std::placeholders::_5);
std::placeholders::_3);
}
void SetFluxVectorVec(DiffusionFluxVecCB fluxVector)
void SetFluxVector(DiffusionFluxVecCB fluxVector)
{
m_fluxVector = fluxVector;
}
template<typename FuncPointerT, typename ObjectPointerT>
template<typename FuncPointerT, typename ObjectPointerT>
void SetFluxVectorNS(FuncPointerT func, ObjectPointerT obj)
{
m_fluxVectorNS = std::bind(
......@@ -114,16 +118,22 @@ namespace Nektar
std::placeholders::_3);
}
void SetFluxVectorNS(DiffusionFluxVecCBNS fluxVector)
{
m_fluxVectorNS = fluxVector;
}
template<typename FuncPointerT, typename ObjectPointerT>
void SetArtificialDiffusionVector(FuncPointerT func, ObjectPointerT obj)
void SetFluxPenaltyNS(FuncPointerT func, ObjectPointerT obj)
{
m_ArtificialDiffusionVector = std::bind(
func, obj, std::placeholders::_1, std::placeholders::_2);
m_fluxPenaltyNS = std::bind(
func, obj, std::placeholders::_1, std::placeholders::_2,
std::placeholders::_3);
}
void SetFluxVectorNS(DiffusionFluxVecCBNS fluxVector)
void SetFluxPenaltyNS(DiffusionFluxPenaltyNS flux)
{
m_fluxVectorNS = fluxVector;
m_fluxPenaltyNS = flux;
}
inline void SetHomoDerivs(Array<OneD, Array<OneD, NekDouble> > &deriv)
......@@ -135,11 +145,11 @@ namespace Nektar
{
return v_GetFluxTensor();
}
protected:
DiffusionFluxVecCB m_fluxVector;
DiffusionFluxVecCBNS m_fluxVectorNS;
DiffusionArtificialDiffusion m_ArtificialDiffusionVector;
DiffusionFluxPenaltyNS m_fluxPenaltyNS;
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
......@@ -147,9 +157,9 @@ namespace Nektar
{
boost::ignore_unused(pSession, pFields);
};
virtual void v_Diffuse(
const int nConvectiveFields,
const std::size_t nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
......@@ -161,17 +171,17 @@ namespace Nektar
{
boost::ignore_unused(deriv);
}
virtual Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &v_GetFluxTensor()
{
static Array<OneD, Array<OneD, Array<OneD, NekDouble> > > tmp;
return tmp;
}
};
};
/// A shared pointer to an EquationSystem object
typedef std::shared_ptr<Diffusion> DiffusionSharedPtr;
/// Datatype of the NekFactory used to instantiate classes derived
/// from the Diffusion class.
typedef LibUtilities::NekFactory<std::string, Diffusion, std::string> DiffusionFactory;
......
......@@ -72,7 +72,7 @@ namespace Nektar
};
/**
* @brief Diffusion3DHomogeneous1D uses the 2D WeakDG approach
* @brief Diffusion3DHomogeneous1D uses the 2D WeakDG approach
* to compute the diffusion term looping on the planes in the z
* direction and adding the flux in z direction at the end.
*/
......@@ -80,12 +80,12 @@ namespace Nektar
{
// Strip trailing string "3DHomogeneous1D" to determine 2D diffusion
// type, and create a diffusion object for the plane.
string name = diffType.substr(0, diffType.length()-15);
m_planeDiff = GetDiffusionFactory().CreateInstance(name, name);
m_diffType = diffType.substr(0, diffType.length()-15);
m_planeDiff = GetDiffusionFactory().CreateInstance(m_diffType, m_diffType);
}
/**
* @brief Initiliase Diffusion3DHomogeneous1D objects and store
* @brief Initiliase Diffusion3DHomogeneous1D objects and store
* them before starting the time-stepping.
*
* @param pSession Pointer to session reader.
......@@ -114,7 +114,26 @@ namespace Nektar
m_homoLen = pFields[0]->GetHomoLen();
m_trans = pFields[0]->GetTransposition();
m_planeCounter = 0;
m_planeDiff->SetFluxVectorNS(m_fluxVectorNS);
if (m_diffType == "LDG")
{
// Set viscous flux for LDG
m_planeDiff->SetFluxVector(m_fluxVector);
}
else if (m_diffType == "LDGNS")
{
// Set viscous flux for LDGNS
m_planeDiff->SetFluxVectorNS(m_fluxVectorNS);
// Set penalty flux
m_planeDiff->SetFluxPenaltyNS(m_fluxPenaltyNS);
}
else if (m_diffType == "LFRDGNS" ||
m_diffType == "LFRHUNS" ||
m_diffType == "LFRSDNS" )
{
// Set viscous flux for FR cases
m_planeDiff->SetFluxVectorNS(m_fluxVectorNS);
}
m_fieldsPlane = Array<OneD, MultiRegions::ExpListSharedPtr>
(nConvectiveFields);
......@@ -169,7 +188,7 @@ namespace Nektar
* using an LDG interface flux and the the flux in the third direction.
*/
void Diffusion3DHomogeneous1D::v_Diffuse(
const int nConvectiveFields,
const std::size_t nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
......@@ -181,22 +200,21 @@ namespace Nektar
Array<OneD, NekDouble> tmp(m_numPoints), tmp2;
Array<OneD, Array<OneD, NekDouble> > viscHComp;
const int nPointsTot = fields[0]->GetNpoints();
int i, j;
NekDouble beta;
if (m_fluxVectorNS)
{
viscHComp = Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
for (i = 0; i < nConvectiveFields - 1; ++i)
for (int i = 0; i < nConvectiveFields - 1; ++i)
{
fields[0]->PhysDeriv(2, inarray[i], m_homoDerivStore[i]);
viscHComp[i] = Array<OneD, NekDouble>(m_numPoints);
}
}
for (i = 0; i < m_numPlanes; ++i)
for (int i = 0; i < m_numPlanes; ++i)
{
// Set up memory references for fields, inarray and outarray for
// this plane.
......@@ -212,7 +230,7 @@ namespace Nektar
m_outarrayPlane[j] = Array<OneD, NekDouble>(
m_numPointsPlane, tmp2 = outarray[j] + m_planePos[i]);
}
if (m_fluxVectorNS)
{
......@@ -220,15 +238,43 @@ namespace Nektar
}
m_planeDiff->Diffuse(nConvectiveFields,
m_fieldsPlane,
m_inarrayPlane,
m_outarrayPlane);
if (m_diffType == "LDGNS")
{
// Store plane Fwd/Bwd traces
std::size_t nTracePts = m_fieldsPlane[0]->GetTrace()
->GetTotPoints();
std::size_t nScalar = m_inarrayPlane.num_elements();
Array<OneD, Array<OneD, NekDouble> > Fwd(nScalar);
Array<OneD, Array<OneD, NekDouble> > Bwd(nScalar);
{
for(std::size_t k = 0; k < nScalar; ++k)
{
Fwd[k] = Array<OneD, NekDouble>(nTracePts, 0.0);
Bwd[k] = Array<OneD, NekDouble>(nTracePts, 0.0);
m_fieldsPlane[k]->GetFwdBwdTracePhys(
m_inarrayPlane[k], Fwd[k], Bwd[k]);
}
}
m_planeDiff->Diffuse(nConvectiveFields,
m_fieldsPlane,
m_inarrayPlane,
m_outarrayPlane,
Fwd,
Bwd);
}
else
{
m_planeDiff->Diffuse(nConvectiveFields,
m_fieldsPlane,
m_inarrayPlane,
m_outarrayPlane);
}
if (m_fluxVectorNS)
{
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &viscTensor = m_planeDiff->GetFluxTensor();
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
&viscTensor = m_planeDiff->GetFluxTensor();
// Extract H (viscTensor[2])
for (int j = 0; j < nConvectiveFields - 1; ++j)
......@@ -239,24 +285,25 @@ namespace Nektar
}
}
}
if (m_fluxVectorNS)
{
for (j = 0; j < nConvectiveFields - 1; ++j)
for (int j = 0; j < nConvectiveFields - 1; ++j)
{
fields[j+1]->PhysDeriv(2, viscHComp[j], tmp);
Vmath::Vadd(nPointsTot, outarray[j+1], 1, tmp, 1, outarray[j+1], 1);
Vmath::Vadd(nPointsTot, outarray[j+1], 1, tmp, 1,
outarray[j+1], 1);
}
}
else
{
for (j = 0; j < nConvectiveFields; ++j)
for (int j = 0; j < nConvectiveFields; ++j)
{
fields[j]->HomogeneousFwdTrans(inarray[j], tmp);
for (i = 0; i < m_numPlanes; ++i)
for (int i = 0; i < m_numPlanes; ++i)
{
beta = 2*M_PI*m_trans->GetK(i)/m_homoLen;
beta *= beta;
......
......@@ -51,7 +51,7 @@ namespace Nektar
new Diffusion3DHomogeneous1D(diffType));
}
static std::string type[];
protected:
Diffusion3DHomogeneous1D(std::string diffType);
......@@ -59,10 +59,10 @@ namespace Nektar
std::string m_diffType;
SolverUtils::DiffusionSharedPtr m_planeDiff;
NekDouble m_homoLen;
int m_numPoints;
int m_numPointsPlane;
int m_numPlanes;
int m_planeCounter;
std::size_t m_numPoints;
std::size_t m_numPointsPlane;
std::size_t m_numPlanes;
std::size_t m_planeCounter;
Array<OneD, unsigned int> m_planes;
Array<OneD, unsigned int> m_planePos;
Array<OneD, Array<OneD, NekDouble> > m_homoDerivStore;
......@@ -76,14 +76,14 @@ namespace Nektar
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_Diffuse(
const int nConvective,
const std::size_t nConvective,
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 = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
};
};
}
}
#endif
......@@ -51,59 +51,63 @@ namespace Nektar
boost::ignore_unused(diffType);
return DiffusionSharedPtr(new DiffusionLDG());
}
static std::string type;
protected:
DiffusionLDG();
std::string m_shockCaptureType;
std::string m_shockCaptureType;
/// Coefficient of penalty term
NekDouble m_C11;
Array<OneD, Array<OneD, NekDouble> > m_traceNormals;
LibUtilities::SessionReaderSharedPtr m_session;
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_Diffuse(
const int nConvective,
const std::size_t nConvective,
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 = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
virtual void v_NumFluxforScalar(
void 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 = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
virtual void v_WeakPenaltyforScalar(
void ApplyScalarBCs(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const int var,
const std::size_t var,
const Array<OneD, const NekDouble> &ufield,
const Array<OneD, const NekDouble> &uplus,
const Array<OneD, const NekDouble> &Fwd,
const Array<OneD, const NekDouble> &Bwd,
Array<OneD, NekDouble> &penaltyflux);
virtual void v_NumFluxforVector(
void 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);
virtual void v_WeakPenaltyforVector(
void ApplyVectorBCs(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const int var,
const int dir,
const std::size_t var,
const std::size_t dir,
const Array<OneD, const NekDouble> &qfield,
const Array<OneD, const NekDouble> &qtemp,
Array<OneD, NekDouble> &penaltyflux,
NekDouble C11);
};
const Array<OneD, const NekDouble> &qFwd,
const Array<OneD, const NekDouble> &qBwd,
Array<OneD, NekDouble> &penaltyflux);
};
}
}
#endif