Skip to content

GitLab

Commit 05ad0d4b authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'master' into feature/ann_remove

parents 20db9d1a e52b1536
Hide whitespace changes
Inline Side-by-side
docs/user-guide/solvers/compressible-flow.tex
View file @ 05ad0d4b
......@@ -226,7 +226,9 @@ conditions are employed (i.e. $T_{w}$). Default value = 300.15$K$;
\item \inltt{wInf} farfield $Z$-component of the velocity (i.e. $w_{\infty}$). Default value = 0.0 $m/s$;
\item \inltt{mu} dynamic viscosity (i.e. $\mu_{\infty}$). Default value = 1.78e-05 $Pa s$;
\item \inltt{Pr} Prandtl number. Default value = 0.72;
\item \inltt{thermalConductivity} thermal conductivity (i.e. $\kappa_{\infty}$). Default value = 0.0257 $W / (K m)$;
\item \inltt{thermalConductivity} thermal conductivity (i.e. $\kappa_{\infty}$). This can be set as an
alternative to \inltt{Pr}, in which case the Prandtl number is calculated from $\kappa_{\infty}$
(it is only possible to set one of them). By default, this is obtained from the Prandtl number;
\end{itemize}
\subsection*{Solver info}
......@@ -248,6 +250,8 @@ Under this section it is possible to set the solver information.
\begin{itemize}
\item \inltt{NavierStokesCFE} (Compressible Navier-Stokes equations);
\item \inltt{EulerCFE} (Compressible Euler equations).
\item \inltt{IsentropicVortex} (Isentropic vortex test case).
\item \inltt{RinglebFlow} (Ringleb flow test case).
\end{itemize}
\item \inltt{Projection} is the type of projection we want to use:
\begin{itemize}
......@@ -263,14 +267,15 @@ Note that the Continuous projection is not supported in the Compressible Flow So
\item \inltt{FRcmin} (Flux-Reconstruction with $c = c_{min}$);
\item \inltt{FRcinf} (Flux-Reconstruction with $c = \infty$).
\end{itemize}
\item \inltt{DiffusionType} is the diffusion operator we want to use:
\item \inltt{DiffusionType} is the diffusion operator we want to use
for the Navier-Stokes equations:
\begin{itemize}
\item \inltt{WeakDG} (classical DG in weak form);
\item \inltt{FRDG} (Flux-Reconstruction recovering nodal DG scheme);
\item \inltt{FRSD} (Flux-Reconstruction recovering a spectral difference (SD) scheme);
\item \inltt{FRHU} (Flux-Reconstruction recovering Huynh (G2) scheme);
\item \inltt{FRcmin} (Flux-Reconstruction with $c = c_{min}$);
\item \inltt{FRcinf} (Flux-Reconstruction with $c = \infty$).
\item \inltt{LDGNS} (LDG);
\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);
\item \inltt{LFRcminNS} (Flux-Reconstruction with $c = c_{min}$);
\item \inltt{LFRcinfNS} (Flux-Reconstruction with $c = \infty$).
\end{itemize}
\item \inltt{TimeIntegrationMethod} is the time-integration scheme we want to use.
Note that only an explicit discretisation is supported:
......@@ -293,14 +298,6 @@ we want to use for the advection operator:
\item \inltt{LaxFriedrichs};
\item \inltt{Roe}.
\end{itemize}
\item \inltt{ProblemType} is the problem type we want to solve.
This tag is supported for solving ad hoc problems such as the
isentropic vortex or the Ringleb flow.
\begin{itemize}
\item \inltt{General};
\item \inltt{IsentropicVortex};
\item \inltt{RinglebFlow};\\[0.2em]
\end{itemize}
\item \inltt{ViscosityType} is the viscosity type we want to use:
\begin{itemize}
\item \inltt{Constant} (Constant viscosity);
......
library/LibUtilities/BasicUtils/NekFactory.hpp
View file @ 05ad0d4b
......@@ -55,7 +55,7 @@
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#ifndef MAX_PARAM
#define MAX_PARAM 5 // default maximum number of parameters to support
#define MAX_PARAM 6 // default maximum number of parameters to support
#endif
namespace Nektar
......
library/SolverUtils/Advection/Advection.cpp
View file @ 05ad0d4b
......@@ -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);
}
......
library/SolverUtils/Advection/Advection.h
View file @ 05ad0d4b
......@@ -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(
......
library/SolverUtils/Advection/Advection3DHomogeneous1D.cpp
View file @ 05ad0d4b
......@@ -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();
......
library/SolverUtils/Advection/Advection3DHomogeneous1D.h
View file @ 05ad0d4b
......@@ -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(
......
library/SolverUtils/Advection/AdvectionFR.cpp
View file @ 05ad0d4b
......@@ -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;
......
library/SolverUtils/Advection/AdvectionFR.h
View file @ 05ad0d4b
......@@ -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,
......
library/SolverUtils/Advection/AdvectionNonConservative.cpp
View file @ 05ad0d4b
......@@ -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();
......
library/SolverUtils/Advection/AdvectionNonConservative.h
View file @ 05ad0d4b
......@@ -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);
};
}
}
......
library/SolverUtils/Advection/AdvectionWeakDG.cpp
View file @ 05ad0d4b
......@@ -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,25 @@ 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];
numflux[i] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
}
}
m_riemann->Solve(m_spaceDim, Fwd, Bwd, numflux);
......
library/SolverUtils/Advection/AdvectionWeakDG.h
View file @ 05ad0d4b
......@@ -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);
};
}
}
......
library/SolverUtils/Diffusion/Diffusion.cpp
View file @ 05ad0d4b
......@@ -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);
}
}
}
library/SolverUtils/Diffusion/Diffusion.h
View file @ 05ad0d4b
......@@ -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> > >
......@@ -113,11 +115,6 @@ namespace Nektar
m_fluxVectorNS = fluxVector;
}
inline void SetRiemannSolver(RiemannSolverSharedPtr riemann)
{
m_riemann = riemann;
}
inline void SetHomoDerivs(Array<OneD, Array<OneD, NekDouble> > &deriv)
{
v_SetHomoDerivs(deriv);
......@@ -131,7 +128,6 @@ namespace Nektar
protected:
DiffusionFluxVecCB m_fluxVector;
DiffusionFluxVecCBNS m_fluxVectorNS;
RiemannSolverSharedPtr m_riemann;
DiffusionArtificialDiffusion m_ArtificialDiffusionVector;
virtual void v_InitObject(
......@@ -145,7 +141,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)
......
library/SolverUtils/Diffusion/Diffusion3DHomogeneous1D.cpp
View file @ 05ad0d4b
......@@ -115,28 +115,9 @@ namespace Nektar
m_planeCounter = 0;
m_planeDiff->SetFluxVectorNS(m_fluxVectorNS);
if (m_riemann)
{
// Set Riemann solver and flux vector callback for this plane.
m_planeDiff->SetRiemannSolver(m_riemann);
// Override Riemann solver scalar and vector callbacks.
map<string, RSScalarFuncType>::iterator it1;
map<string, RSScalarFuncType> scalars = m_riemann->GetScalars();
for (it1 = scalars.begin(); it1 != scalars.end(); ++it1)
{
boost::shared_ptr<HomoRSScalar> tmp =
MemoryManager<HomoRSScalar>
::AllocateSharedPtr(it1->second, m_numPlanes);
m_riemann->SetScalar(it1->first, &HomoRSScalar::Exec, tmp);
}
}
m_fieldsPlane = Array<OneD, MultiRegions::ExpListSharedPtr>
(nConvectiveFields);
if (m_fluxVectorNS)
{
m_inarrayPlane = Array<OneD, Array<OneD, NekDouble> >
......@@ -190,7 +171,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;
......
library/SolverUtils/Diffusion/Diffusion3DHomogeneous1D.h
View file @ 05ad0d4b
......@@ -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);
};
}
}
......
library/SolverUtils/Diffusion/DiffusionLDG.cpp
View file @ 05ad0d4b
......@@ -43,21 +43,21 @@ namespace Nektar
{
std::string DiffusionLDG::type = GetDiffusionFactory().
RegisterCreatorFunction("LDG", DiffusionLDG::create);
DiffusionLDG::DiffusionLDG()
{
}
void DiffusionLDG::v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
m_session = pSession;
m_session->LoadSolverInfo("ShockCaptureType",
m_session->LoadSolverInfo("ShockCaptureType",
m_shockCaptureType, "Off");
// Setting up the normals
// Setting up the normals
int i;
int nDim = pFields[0]->GetCoordim(0);
int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
......@@ -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);
......@@ -83,7 +85,6 @@ namespace Nektar
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
Array<OneD, NekDouble> qcoeffs(nCoeffs);
Array<OneD, NekDouble> temp (nCoeffs);
Array<OneD, Array<OneD, NekDouble> > fluxvector(nDim);
......@@ -116,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)
{
......@@ -220,13 +221,15 @@ 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();
int nvariables = fields.num_elements();
int nDim = uflux.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);
......@@ -242,48 +245,55 @@ 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);
// 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
{
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
if(fields[0]->GetBndCondExpansions().num_elements())
{
v_WeakPenaltyforScalar(fields, i, ufield[i], 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
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})
// 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,
......@@ -308,14 +318,14 @@ namespace Nektar
int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
Array<OneD, NekDouble > uplus(nTracePts);
fields[var]->ExtractTracePhys(ufield, uplus);
for (i = 0; i < nBndRegions; ++i)
{
// Number of boundary expansion related to that region
nBndEdges = fields[var]->
GetBndCondExpansions()[i]->GetExpSize();
// Weakly impose boundary conditions by modifying flux values
for (e = 0; e < nBndEdges ; ++e)
{
......@@ -392,6 +402,17 @@ namespace Nektar
// qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
for (i = 0; i < nvariables; ++i)
{
// Generate Stability term = - C11 ( u- - u+ )
fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
Vmath::Vsub(nTracePts,
Fwd, 1, Bwd, 1,
uterm, 1);
Vmath::Smul(nTracePts,
-1.0 * C11, uterm, 1,
uterm, 1);
qflux[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
for (j = 0; j < nDim; ++j)
{
......@@ -418,18 +439,7 @@ namespace Nektar
m_traceNormals[j], 1,
qfluxtemp, 1,
qfluxtemp, 1);
// Generate Stability term = - C11 ( u- - u+ )
fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
Vmath::Vsub(nTracePts,
Fwd, 1, Bwd, 1,
uterm, 1);
Vmath::Smul(nTracePts,
-1.0 * C11, uterm, 1,
uterm, 1);
// Flux = {Fwd, Bwd} * (nx, ny, nz) + uterm * (nx, ny)
Vmath::Vadd(nTracePts,
uterm, 1,
......@@ -474,11 +484,10 @@ namespace Nektar
int nBndEdges, nBndEdgePts;
int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
Array<OneD, NekDouble > uterm(nTracePts);
Array<OneD, NekDouble > qtemp(nTracePts);
int cnt = 0;
/*
// Setting up the normals
m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(nDim);
......
library/SolverUtils/Diffusion/DiffusionLDG.h
View file @ 05ad0d4b
......@@ -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,