Commit a29b079c authored by Gianmarco Mengaldo's avatar Gianmarco Mengaldo
Browse files

Restructuring of AdvectioFR class. Added some regression tests for FR 2D....

Restructuring of AdvectioFR class. Added some regression tests for FR 2D. Added a method to compute the contravariant fluxes.
parent f7d34ccc
......@@ -660,18 +660,21 @@ namespace Nektar
}//end scope
break;
/** \brief Left derivative of the correction function for FR DG method (see J Sci Comput (2011) 47: 50–72)
/** \brief Left derivative of the correction function for FR DG
method (see J Sci Comput (2011) 47: 50–72)
\f$\tilde dGL_{DG} = (-1)^{p}/2( d(L_(p)) - d(L_{p+1}) )\f$
*/
case eDG_DG_Left:
{
// Number of modes (here the left polynomial DG_DG_Left of degree p is stored)
// Number of modes (here the left polynomial DG_DG_Left
// of degree p is stored)
mode = m_bdata.data();
int p = numModes - 1;
// Auxiliary vectors to build up the auxiliary derivatives of the Legendre polynomials
// Auxiliary vectors to build up the auxiliary derivatives
// of the Legendre polynomials
Array<OneD,NekDouble> dLp (numPoints, 0.0);
Array<OneD,NekDouble> dLpp (numPoints, 0.0);
Array<OneD,NekDouble> dLpm (numPoints, 0.0);
......@@ -702,18 +705,21 @@ namespace Nektar
break;
/** \brief Right derivative of the correction function for FR DG method (see J Sci Comput (2011) 47: 50–72)
/** \brief Right derivative of the correction function for FR DG
method (see J Sci Comput (2011) 47: 50–72)
\f$\tilde dGR_{DG} = 1/2( d(L_{p}) + d(L_{p+1}) ) \f$
*/
case eDG_DG_Right:
{
// Number of modes (here the right polynomial DG_DG_Right of degree p is stored)
// Number of modes (here the right polynomial DG_DG_Right
// of degree p is stored)
mode = m_bdata.data();
int p = numModes - 1;
// Auxiliary vectors to build up the auxiliary derivatives of the Legendre polynomials
// Auxiliary vectors to build up the auxiliary derivatives
// of the Legendre polynomials
Array<OneD,NekDouble> dLp (numPoints, 0.0);
Array<OneD,NekDouble> dLpp (numPoints, 0.0);
Array<OneD,NekDouble> dLpm (numPoints, 0.0);
......@@ -745,18 +751,22 @@ namespace Nektar
break;
/** \brief Left derivative of the correction function for FR SD method (see J Sci Comput (2011) 47: 50–72)
/** \brief Left derivative of the correction function for FR SD
method (see J Sci Comput (2011) 47: 50–72)
\f$\tilde dGL_{SD} = (-1)^{p}/2( d(L_{p}) - (pd(L_{p-1}) + (p+1)d(L_{p+1}))/(2p + 1) ) \f$
\f$\tilde dGL_{SD} = (-1)^{p}/2( d(L_{p}) - (pd(L_{p-1}) +
(p+1)d(L_{p+1}))/(2p + 1) ) \f$
*/
case eDG_SD_Left:
{
// Number of modes (here the left polynomial DG_SD_Left of degree p is stored)
// Number of modes (here the left polynomial DG_SD_Left of
// degree p is stored)
mode = m_bdata.data();
int p = numModes - 1;
// Auxiliary vectors to build up the auxiliary derivatives of the Legendre polynomials
// Auxiliary vectors to build up the auxiliary derivatives
// of the Legendre polynomials
Array<OneD,NekDouble> dLp (numPoints, 0.0);
Array<OneD,NekDouble> dLpp (numPoints, 0.0);
Array<OneD,NekDouble> dLpm (numPoints, 0.0);
......@@ -790,18 +800,22 @@ namespace Nektar
}//end scope
break;
/** \brief Right derivative of the correction function for FR SD method (see J Sci Comput (2011) 47: 50–72)
/** \brief Right derivative of the correction function for FR SD
method (see J Sci Comput (2011) 47: 50–72)
\f$\tilde dGR_{SD} = 1/2( d(L_{p}) + (pd(L_{p-1}) + (p+1)d(L_{p+1}))/(2p + 1) ) \f$
\f$\tilde dGR_{SD} = 1/2( d(L_{p}) + (pd(L_{p-1}) +
(p+1)d(L_{p+1}))/(2p + 1) ) \f$
*/
case eDG_SD_Right:
{
// Number of modes (here the right polynomials DG_SD_Right of degree p is stored)
// Number of modes (here the right polynomials DG_SD_Right
// of degree p is stored)
mode = m_bdata.data();
int p = numModes - 1;
// Auxiliary vectors to build up the auxiliary derivatives of the Legendre polynomials
// Auxiliary vectors to build up the auxiliary derivatives
// of the Legendre polynomials
Array<OneD,NekDouble> dLp (numPoints, 0.0);
Array<OneD,NekDouble> dLpp (numPoints, 0.0);
Array<OneD,NekDouble> dLpm (numPoints, 0.0);
......@@ -836,18 +850,22 @@ namespace Nektar
break;
/** \brief Left derivative of the correction function for FR HU method (see J Sci Comput (2011) 47: 50–72)
/** \brief Left derivative of the correction function for FR HU
method (see J Sci Comput (2011) 47: 50–72)
\f$\tilde dGL_{HU} = (-1)^{p}/2( d(L_{p}) - ((p+1)d(L_{p-1}) + pd(L_{p+1}))/(2p + 1) ) \f$
\f$\tilde dGL_{HU} = (-1)^{p}/2( d(L_{p}) - ((p+1)d(L_{p-1})
+ pd(L_{p+1}))/(2p + 1) ) \f$
*/
case eDG_HU_Left:
{
// Number of modes (here the left polynomial DG_HU_Left of degree p is stored)
// Number of modes (here the left polynomial DG_HU_Left of
// degree p is stored)
mode = m_bdata.data();
int p = numModes - 1;
// Auxiliary vectors to build up the auxiliary derivatives of the Legendre polynomials
// Auxiliary vectors to build up the auxiliary derivatives
// of the Legendre polynomials
Array<OneD,NekDouble> dLp (numPoints, 0.0);
Array<OneD,NekDouble> dLpp (numPoints, 0.0);
Array<OneD,NekDouble> dLpm (numPoints, 0.0);
......@@ -881,18 +899,22 @@ namespace Nektar
}//end scope
break;
/** \brief Right derivative of the correction function for FR HU method (see J Sci Comput (2011) 47: 50–72)
/** \brief Right derivative of the correction function for FR HU
method (see J Sci Comput (2011) 47: 50–72)
\f$\tilde dGR_{SD} = 1/2( d(L_{p}) + ((p+1)d(L_{p-1}) + pd(L_{p+1}))/(2p + 1) ) \f$
\f$\tilde dGR_{SD} = 1/2( d(L_{p}) + ((p+1)d(L_{p-1}) +
pd(L_{p+1}))/(2p + 1) ) \f$
*/
case eDG_HU_Right:
{
// Number of modes (here the right polynomial DG_HU_Right of degree p is stored)
// Number of modes (here the right polynomial DG_HU_Right of
// degree p is stored)
mode = m_bdata.data();
int p = numModes - 1;
// Auxiliary vectors to build up the auxiliary derivatives of the Legendre polynomials
// Auxiliary vectors to build up the auxiliary derivatives
// of the Legendre polynomials
Array<OneD,NekDouble> dLp (numPoints, 0.0);
Array<OneD,NekDouble> dLpp (numPoints, 0.0);
Array<OneD,NekDouble> dLpm (numPoints, 0.0);
......@@ -925,11 +947,207 @@ namespace Nektar
}
}//end scope
break;
/** \brief Left derivative of the correction function for a
method with a selected c (see J Sci Comput (2011) 47: 50–72)
*/
case eDG_c_Left:
{
// Number of modes
mode = m_bdata.data();
int p = numModes - 1;
// Auxiliary vectors to build up the auxiliary derivatives
// of the Legendre polynomials
Array<OneD,NekDouble> dLp (numPoints, 0.0);
Array<OneD,NekDouble> dLpp (numPoints, 0.0);
Array<OneD,NekDouble> dLpm (numPoints, 0.0);
// Function sign
NekDouble sign = pow(-1.0, p);
// Factors to build the scheme
NekDouble ap;
ap = tgamma(2 * p + 1) / (pow(2.0, p) * tgamma(p + 1) * tgamma(p + 1));
double c = 2 * p / ((2 * p + 1) * (p + 1) * tgamma(p + 1) * tgamma(p + 1));
//double c = 2 * (p + 1) / ((2 * p + 1) * p * tgamma(p + 1) * tgamma(p + 1));
double c_p = -2/((2 * p + 1) * (ap * tgamma(p + 1))*(ap * tgamma(p + 1)));
if (c < c_p)
{
ASSERTL0(false, "c out of bounds");
break;
}
NekDouble etap;
etap = 0.5 * c * (2 * p + 1) * pow (1.0 * ap * tgamma(p + 1), 2 );
NekDouble overeta = 1.0 / (1.0 + etap);
// Derivative of the Legendre polynomials
// dLp = derivative of the Legendre polynomial of order p
// dLpp = derivative of the Legendre polynomial of order p+1
// dLpm = derivative of the Legendre polynomial of order p-1
Polylib::jacobd(numPoints, z.data(), &(dLp[0]), p, 0.0, 0.0);
Polylib::jacobd(numPoints, z.data(), &(dLpp[0]), p+1, 0.0, 0.0);
Polylib::jacobd(numPoints, z.data(), &(dLpm[0]), p-1, 0.0, 0.0);
// Building the DG_c_Left
for(i = 0; i < numPoints; ++i)
{
mode[i] = etap * dLpm[i];
mode[i] += dLpp[i];
mode[i] *= overeta;
mode[i] = dLp[i] - mode[i];
mode[i] = 0.5 * sign * mode[i];
}
}//end scope
break;
/** \brief Left derivative of the correction function for a
method with a selected c (see J Sci Comput (2011) 47: 50–72)
*/
case eDG_c_Right:
{
// Number of modes (here the right polynomial DG_c_Right of
// degree p is stored)
mode = m_bdata.data();
int p = numModes - 1;
// Auxiliary vectors to build up the auxiliary derivatives
// of the Legendre polynomials
Array<OneD,NekDouble> dLp (numPoints, 0.0);
Array<OneD,NekDouble> dLpp (numPoints, 0.0);
Array<OneD,NekDouble> dLpm (numPoints, 0.0);
// Function sign
NekDouble sign = pow(-1.0, p);
// Factors to build the scheme
NekDouble ap;
ap = tgamma(2 * p + 1) / (pow(2.0, p) * tgamma(p + 1) * tgamma(p + 1));
double c = 2 * p / ((2 * p + 1) * (p + 1) * tgamma(p + 1) * tgamma(p + 1));
//double c = 2 * (p + 1) / ((2 * p + 1) * p * tgamma(p + 1) * tgamma(p + 1));
double c_p = -2/((2 * p + 1) * (ap * tgamma(p + 1))*(ap * tgamma(p + 1)));
if (c < c_p)
{
ASSERTL0(false, "c out of bounds");
break;
}
NekDouble etap;
etap = 0.5 * c * (2 * p + 1) * pow (1.0 * ap * tgamma(p + 1), 2 );
NekDouble overeta = 1.0 / (1.0 + etap);
// Derivative of the Legendre polynomials
// dLp = derivative of the Legendre polynomial of order p
// dLpp = derivative of the Legendre polynomial of order p+1
// dLpm = derivative of the Legendre polynomial of order p-1
Polylib::jacobd(numPoints, z.data(), &(dLp[0]), p, 0.0, 0.0);
Polylib::jacobd(numPoints, z.data(), &(dLpp[0]), p+1, 0.0, 0.0);
Polylib::jacobd(numPoints, z.data(), &(dLpm[0]), p-1, 0.0, 0.0);
// Building the DG_c_Right
for(i = 0; i < numPoints; ++i)
{
mode[i] = etap * dLpm[i];
mode[i] += dLpp[i];
mode[i] *= overeta;
mode[i] += dLp[i];
mode[i] = 0.5 * mode[i];
}
}//end scope
break;
/** \brief Left derivative of the correction function for a
method with c equal to infinity (see J Sci Comput (2011)
47: 50–72)
*/
case eDG_inf_Left:
{
// Number of modes (here the left polynomial DG_DG_Left of
// degree p is stored)
mode = m_bdata.data();
int p = numModes - 2;
// Auxiliary vectors to build up the auxiliary derivatives
// of the Legendre polynomials
Array<OneD,NekDouble> dLp (numPoints, 0.0);
Array<OneD,NekDouble> dLpp (numPoints, 0.0);
Array<OneD,NekDouble> dLpm (numPoints, 0.0);
// Function sign to build up DG_DG_Left
NekDouble sign = pow(-1.0, p - 1);
// Derivative of the Legendre polynomials
// dLp = derivative of the Legendre polynomial of order p
// dLpp = derivative of the Legendre polynomial of order p+1
// dLpm = derivative of the Legendre polynomial of order p-1
Polylib::jacobd(numPoints, z.data(), &(dLp[0]), p, 0.0, 0.0);
Polylib::jacobd(numPoints, z.data(), &(dLpp[0]), p+1, 0.0, 0.0);
Polylib::jacobd(numPoints, z.data(), &(dLpm[0]), p-1, 0.0, 0.0);
// Building the DG_DG_Left
for(int i = 0; i < numPoints; ++i)
{
mode[i] = dLpm[i];
mode[i] = mode[i] - dLp[i];
mode[i] = 0.5 * sign * m_bdata[i];
}
}//end scope
break;
/** \brief Right derivative of the correction function for a
method with c equal to infinity
(see J Sci Comput (2011) 47: 50–72)
*/
case eDG_inf_Right:
{
// Number of modes (here the right polynomial DG_DG_Right
// of degree p is stored)
mode = m_bdata.data();
int p = numModes - 2;
// Auxiliary vectors to build up the auxiliary derivatives
// of the Legendre polynomials
Array<OneD,NekDouble> dLp (numPoints, 0.0);
Array<OneD,NekDouble> dLpp (numPoints, 0.0);
Array<OneD,NekDouble> dLpm (numPoints, 0.0);
// Function sign to build up DG_DG_Right
NekDouble sign = pow(-1.0, p);
// Derivative of the Legendre polynomials
// dLp = derivative of the Legendre polynomial of order p
// dLpp = derivative of the Legendre polynomial of order p+1
// dLpm = derivative of the Legendre polynomial of order p-1
Polylib::jacobd(numPoints, z.data(), &(dLp[0]), p, 0.0, 0.0);
Polylib::jacobd(numPoints, z.data(), &(dLpp[0]), p+1, 0.0, 0.0);
Polylib::jacobd(numPoints, z.data(), &(dLpm[0]), p-1, 0.0, 0.0);
// Building the DG_DG_Right
for(int i = 0; i < numPoints; ++i)
{
mode[i] += dLpm[i];
mode[i] += dLp[i];
mode[i] = 0.5 * mode[i];
}
}//end scope
break;
default:
ASSERTL0(false, "Basis Type not known or "
"not implemented at this time.");
"not implemented at this time.");
}
}
......
......@@ -117,6 +117,10 @@ namespace Nektar
case eDG_SD_Right:
case eDG_HU_Left:
case eDG_HU_Right:
case eDG_c_Left:
case eDG_c_Right:
case eDG_inf_Left:
case eDG_inf_Right:
case eFourierSingleMode:
case eFourierHalfModeRe:
case eFourierHalfModeIm:
......
......@@ -64,6 +64,10 @@ namespace Nektar
eDG_SD_Right, //!< Derivative of the Right correction function for SD FR \f$ dGR_{p}(z_i) \f$
eDG_HU_Left, //!< Derivative of the left correction function for HU FR \f$ dGL_{p}(z_i) \f$
eDG_HU_Right, //!< Derivative of the Right correction function for HU FR \f$ dGR_{p}(z_i) \f$
eDG_c_Left, //!< Derivative of the left correction function for c FR \f$ dGL_{p}(z_i) \f$
eDG_c_Right, //!< Derivative of the Right correction function for c FR \f$ dGR_{p}(z_i) \f$
eDG_inf_Left, //!< Derivative of the left correction function for c_inf FR \f$ dGL_{p}(z_i) \f$
eDG_inf_Right, //!< Derivative of the Right correction function for c_inf FR \f$ dGR_{p}(z_i) \f$
SIZE_BasisType //!< Length of enum list
};
}
......
......@@ -66,8 +66,12 @@ namespace Nektar
"DG_DG_Right",
"DG_SD_Left",
"DG_SD_Right",
"DG_HU_Left"
"DG_HU_Right"
"DG_HU_Left",
"DG_HU_Right",
"DG_c_Left",
"DG_c_Right",
"DG_inf_Left",
"DG_inf_Right"
};
......
......@@ -765,6 +765,137 @@ namespace Nektar
&outarray[0],1);
}
}
void QuadExp::v_GetEdgeQFactors(
const int edge,
Array<OneD, NekDouble> &outarray)
{
int i;
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
const Array<OneD, const NekDouble>& jac = m_metricinfo->GetJac();
const Array<TwoD, const NekDouble>& gmat = m_metricinfo->GetGmat();
Array<OneD, NekDouble> j (max(nquad0, nquad1), 0.0);
Array<OneD, NekDouble> g0(max(nquad0, nquad1), 0.0);
Array<OneD, NekDouble> g1(max(nquad0, nquad1), 0.0);
Array<OneD, NekDouble> g2(max(nquad0, nquad1), 0.0);
Array<OneD, NekDouble> g3(max(nquad0, nquad1), 0.0);
if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
switch (edge)
{
case 0:
Vmath::Vcopy(nquad0, &(gmat[1][0]), 1, &(g1[0]), 1);
Vmath::Vcopy(nquad0, &(gmat[3][0]), 1, &(g3[0]), 1);
Vmath::Vcopy(nquad0, &(jac[0]), 1, &(j[0]), 1);
for (i = 0; i < nquad0; ++i)
{
outarray[i] = j[i]*sqrt(g1[i]*g1[i] + g3[i]*g3[i]);
}
break;
case 1:
Vmath::Vcopy(nquad1,
&(gmat[0][0])+(nquad0-1), nquad0,
&(g0[0]), 1);
Vmath::Vcopy(nquad1,
&(gmat[2][0])+(nquad0-1), nquad0,
&(g2[0]), 1);
Vmath::Vcopy(nquad1,
&(jac[0])+(nquad0-1), nquad0,
&(j[0]), 1);
for (i = 0; i < nquad0; ++i)
{
outarray[i] = j[i]*sqrt(g0[i]*g0[i] + g2[i]*g2[i]);
}
break;
case 2:
Vmath::Vcopy(nquad0,
&(gmat[1][0])+(nquad0*nquad1-1), -1,
&(g1[0]), 1);
Vmath::Vcopy(nquad0,
&(gmat[3][0])+(nquad0*nquad1-1), -1,
&(g3[0]), 1);
Vmath::Vcopy(nquad0,
&(jac[0])+(nquad0*nquad1-1), -1,
&(j[0]), 1);
for (i = 0; i < nquad0; ++i)
{
outarray[i] = j[i]*sqrt(g1[i]*g1[i] + g3[i]*g3[i]);
}
break;
case 3:
Vmath::Vcopy(nquad1,
&(gmat[0][0])+nquad0*(nquad1-1), -nquad0,
&(g0[0]), 1);
Vmath::Vcopy(nquad1,
&(gmat[2][0])+nquad0*(nquad1-1), -nquad0,
&(g2[0]), 1);
Vmath::Vcopy(nquad1,
&(jac[0])+nquad0*(nquad1-1), -nquad0,
&(j[0]), 1);
for (i = 0; i < nquad0; ++i)
{
outarray[i] = j[i]*sqrt(g0[i]*g0[i] + g2[i]*g2[i]);
}
break;
default:
ASSERTL0(false,"edge value (< 3) is out of range");
break;
}
}
else
{
switch (edge)
{
case 0:
for (i = 0; i < nquad0; ++i)
{
outarray[i] = jac[0]*sqrt(gmat[1][0]*gmat[1][0] +
gmat[3][0]*gmat[3][0]);
}
break;
case 1:
for (i = 0; i < nquad1; ++i)
{
outarray[i] = jac[0]*sqrt(gmat[0][0]*gmat[0][0] +
gmat[2][0]*gmat[2][0]);
}
break;
case 2:
for (i = 0; i < nquad0; ++i)
{
outarray[i] = jac[0]*sqrt(gmat[1][0]*gmat[1][0] +
gmat[3][0]*gmat[3][0]);
}
break;
case 3:
for (i = 0; i < nquad1; ++i)
{
outarray[i] = jac[0]*sqrt(gmat[0][0]*gmat[0][0] +
gmat[2][0]*gmat[2][0]);
}
break;
default:
ASSERTL0(false,"edge value (< 3) is out of range");
break;
}
}
}
void QuadExp::v_ComputeEdgeNormal(const int edge)
......
......@@ -58,9 +58,9 @@ namespace Nektar
* points and order definition
*/
LOCAL_REGIONS_EXPORT QuadExp(
const LibUtilities::BasisKey &Ba,
const LibUtilities::BasisKey &Bb,
const SpatialDomains::QuadGeomSharedPtr &geom);
const LibUtilities::BasisKey &Ba,
const LibUtilities::BasisKey &Bb,
const SpatialDomains::QuadGeomSharedPtr &geom);
LOCAL_REGIONS_EXPORT QuadExp(const QuadExp &T);
......@@ -71,101 +71,108 @@ namespace Nektar
// Integration Methods
//-------------------------------
LOCAL_REGIONS_EXPORT virtual NekDouble v_Integral(
const Array<OneD, const NekDouble> &inarray);
const Array<OneD, const NekDouble> &inarray);
//----------------------------
// Differentiation Methods
//----------------------------
LOCAL_REGIONS_EXPORT virtual void v_PhysDeriv(const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &out_d0,
Array<OneD, NekDouble> &out_d1,
Array<OneD, NekDouble> &out_d2 = NullNekDouble1DArray);
LOCAL_REGIONS_EXPORT virtual void v_PhysDeriv(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &out_d0,
Array<OneD, NekDouble> &out_d1,
Array<OneD, NekDouble> &out_d2 = NullNekDouble1DArray);
LOCAL_REGIONS_EXPORT virtual void v_PhysDeriv(const int dir,
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble> &outarray);
LOCAL_REGIONS_EXPORT virtual void v_PhysDeriv(
const int dir,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray);
LOCAL_REGIONS_EXPORT virtual void v_PhysDirectionalDeriv(const Array<OneD, const NekDouble> &inarray,
const Array<OneD, const NekDouble>& direction,
Array<OneD, NekDouble> &out);
LOCAL_REGIONS_EXPORT virtual void v_PhysDirectionalDeriv(