Commit 18659d42 authored by Spencer Sherwin's avatar Spencer Sherwin
Browse files

Updates to add in SVV capability to Pyramids using the modified expansion....

Updates to add in SVV capability to Pyramids using the modified expansion. This also required the introduction of IProductWRTDerivBasis into PyrExp.cpp. A regression test was also included
parent 5cb1a2f7
......@@ -335,6 +335,21 @@ int main(int argc, char *argv[]){
E = F;
E->GetCoords(x,y,z);
const LibUtilities::BasisKey OBkey1(eModified_A,order1,Pkey1);
const LibUtilities::BasisKey OBkey2(eModified_A,order2,Pkey2);
const LibUtilities::BasisKey OBkey3(eModifiedPyr_C,order3,Pkey3);
StdRegions::StdPyrExp *F1 = new StdRegions::StdPyrExp(OBkey1,OBkey2,OBkey3);
Array<OneD, NekDouble> coeffs1(order1*order2*order3,0.0);
coeffs1[0] = 0.0;
coeffs1[4] = 1.0;
F1->BwdTrans(coeffs1,sol);
E->FwdTrans(sol,coeffs1);
E->BwdTrans(coeffs1,sol);
//----------------------------------------------
// Define solution to be projected
for(i = 0; i < nq1*nq2*nq3; ++i)
......@@ -486,7 +501,7 @@ NekDouble Pyr_sol(NekDouble x, NekDouble y, NekDouble z,
{
for(l = 0; l < order2-k; ++l)
{
for(m = 0; m < order3-k-l; ++m)
for(m = 0; m < order3-k-l-1; ++m)
{
sol += pow(x,k)*pow(y,l)*pow(z,m);
}
......
......@@ -289,7 +289,7 @@ namespace Nektar
/** \brief Orthogonal basis C
\f$\tilde \psi_{pqr}^c = \left ( {1 - \eta_3} \over 2 \right)^{p+q} P_r^{2p+2q+2, 0}(\eta_3)\f$ \\
\f$\tilde \psi_{pqr}^c = \left ( {1 - \eta_3} \over 2 \right)^{p+q} P_r^{2p+2q+2, 0}(\eta_3)\f$ \ \
*/
......@@ -329,15 +329,13 @@ namespace Nektar
break;
/** \brief Orthogonal basis C for Pyramid expansion (which is richer than tets)
\f$\tilde \psi_{pqr}^c = \left ( {1 - \eta_3} \over 2 \right)^{p+q} P_r^{2p+2q+2, 0}(\eta_3)\f$ \\
*/
// This is tilde psi_pqr in Spen's book, page 105
// The 4-dimensional array is laid out in memory such that
// 1) Eta_z values are the changing the fastest, then r, q, and finally p.
// 2) r index increases by the stride of numPoints.
\f$\tilde \psi_{pqr}^c = \left ( {1 - \eta_3} \over 2\right)^{pq} P_r^{2pq+2, 0}(\eta_3)\f$
\f$ \mbox{where }pq = max(p+q-1,0) \f$
*/
// 1) Eta_z values are the changing the fastest, then r, q, and finally p.
// 2) r index increases by the stride of numPoints.
case eOrthoPyr_C:
{
int P = numModes - 1, Q = numModes - 1, R = numModes - 1;
......@@ -349,14 +347,23 @@ namespace Nektar
{
for( int r = 0; r <= R - max(p,q); ++r, mode += numPoints )
{
Polylib::jacobfd(numPoints, z.data(), mode, NULL, r, 2*p + 2*q + 2.0, 0.0);
// this offset allows for orthogonal
// expansion to span linear FE space
// of modified basis but means that
// the cartesian polynomial space
// spanned by the expansion is one
// order lower.
//int pq = max(p + q -1,0);
int pq = max(p + q,0);
Polylib::jacobfd(numPoints, z.data(), mode, NULL, r, 2*pq + 2.0, 0.0);
for( int k = 0; k < numPoints; ++k )
{
// Note factor of 0.5 is part of normalisation
mode[k] *= pow(0.5*(1.0 - z[k]), p+q);
mode[k] *= pow(0.5*(1.0 - z[k]), pq);
// finish normalisation
mode[k] *= sqrt(r+p+q+1.5);
mode[k] *= sqrt(r+pq+1.5);
}
}
}
......
......@@ -273,6 +273,51 @@ namespace Nektar
* \sum_{k=0}^{nq_0} \psi_{p}^a (\xi_{3k}) g_{pq} (\xi_{3k}) = {\bf
* B_1 G} \f$
*/
void PyrExp::v_IProductWRTBase(
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble>& outarray)
{
v_IProductWRTBase_SumFac(inarray, outarray);
}
void PyrExp::v_IProductWRTBase_SumFac(
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble>& outarray,
bool multiplybyweights)
{
const int nquad0 = m_base[0]->GetNumPoints();
const int nquad1 = m_base[1]->GetNumPoints();
const int nquad2 = m_base[2]->GetNumPoints();
const int order0 = m_base[0]->GetNumModes();
const int order1 = m_base[1]->GetNumModes();
Array<OneD, NekDouble> wsp(order0*nquad2*(nquad1+order1));
if(multiplybyweights)
{
Array<OneD, NekDouble> tmp(nquad0*nquad1*nquad2);
MultiplyByQuadratureMetric(inarray, tmp);
IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
m_base[1]->GetBdata(),
m_base[2]->GetBdata(),
tmp,outarray,wsp,
true,true,true);
}
else
{
IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
m_base[1]->GetBdata(),
m_base[2]->GetBdata(),
inarray,outarray,wsp,
true,true,true);
}
}
#if 0
void PyrExp::v_IProductWRTBase(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
......@@ -295,8 +340,165 @@ namespace Nektar
StdPyrExp::v_IProductWRTBase(tmp,outarray);
}
#endif
/**
* @brief Calculates the inner product \f$ I_{pqr} = (u,
* \partial_{x_i} \phi_{pqr}) \f$.
*
* The derivative of the basis functions is performed using the chain
* rule in order to incorporate the geometric factors. Assuming that
* the basis functions are a tensor product
* \f$\phi_{pqr}(\eta_1,\eta_2,\eta_3) =
* \phi_1(\eta_1)\phi_2(\eta_2)\phi_3(\eta_3)\f$, this yields the
* result
*
* \f[
* I_{pqr} = \sum_{j=1}^3 \left(u, \frac{\partial u}{\partial \eta_j}
* \frac{\partial \eta_j}{\partial x_i}\right)
* \f]
*
* In the pyramid element, we must also incorporate a second set
* of geometric factors which incorporate the collapsed co-ordinate
* system, so that
*
* \f[ \frac{\partial\eta_j}{\partial x_i} = \sum_{k=1}^3
* \frac{\partial\eta_j}{\partial\xi_k}\frac{\partial\xi_k}{\partial
* x_i} \f]
*
* These derivatives can be found on p152 of Sherwin & Karniadakis.
*
* @param dir Direction in which to take the derivative.
* @param inarray The function \f$ u \f$.
* @param outarray Value of the inner product.
*/
void PyrExp::v_IProductWRTDerivBase(
const int dir,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
{
v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
}
void PyrExp::v_IProductWRTDerivBase_SumFac(
const int dir,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
{
const int nquad0 = m_base[0]->GetNumPoints();
const int nquad1 = m_base[1]->GetNumPoints();
const int nquad2 = m_base[2]->GetNumPoints();
const int order0 = m_base[0]->GetNumModes ();
const int order1 = m_base[1]->GetNumModes ();
const int nqtot = nquad0*nquad1*nquad2;
int i;
const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
Array<OneD, NekDouble> gfac0(nquad0 );
Array<OneD, NekDouble> gfac1(nquad1 );
Array<OneD, NekDouble> gfac2(nquad2 );
Array<OneD, NekDouble> tmp1 (nqtot );
Array<OneD, NekDouble> tmp2 (nqtot );
Array<OneD, NekDouble> tmp3 (nqtot );
Array<OneD, NekDouble> tmp4 (nqtot );
Array<OneD, NekDouble> tmp5 (nqtot );
Array<OneD, NekDouble> tmp6 (m_ncoeffs);
Array<OneD, NekDouble> wsp (std::max(nqtot,order0*nquad2*(nquad1+order1)));
const Array<TwoD, const NekDouble>& df =
m_metricinfo->GetDerivFactors(GetPointsKeys());
MultiplyByQuadratureMetric(inarray, tmp1);
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
Vmath::Vmul(nqtot,&df[3*dir][0], 1,tmp1.get(),1,tmp2.get(),1);
Vmath::Vmul(nqtot,&df[3*dir+1][0],1,tmp1.get(),1,tmp3.get(),1);
Vmath::Vmul(nqtot,&df[3*dir+2][0],1,tmp1.get(),1,tmp4.get(),1);
}
else
{
Vmath::Smul(nqtot, df[3*dir][0], tmp1.get(),1,tmp2.get(), 1);
Vmath::Smul(nqtot, df[3*dir+1][0],tmp1.get(),1,tmp3.get(), 1);
Vmath::Smul(nqtot, df[3*dir+2][0],tmp1.get(),1,tmp4.get(), 1);
}
// set up geometric factor: (1+z0)/2
for (i = 0; i < nquad0; ++i)
{
gfac0[i] = 0.5*(1+z0[i]);
}
// set up geometric factor: (1+z1)/2
for(i = 0; i < nquad1; ++i)
{
gfac1[i] = 0.5*(1+z1[i]);
}
// Set up geometric factor: 2/(1-z2)
for (i = 0; i < nquad2; ++i)
{
gfac2[i] = 2.0/(1-z2[i]);
}
const int nq01 = nquad0*nquad1;
for (i = 0; i < nquad2; ++i)
{
Vmath::Smul(nq01,gfac2[i],&tmp2[0]+i*nq01,1,&tmp2[0]+i*nq01,1); // 2/(1-z2) for d/dxi_0
Vmath::Smul(nq01,gfac2[i],&tmp3[0]+i*nq01,1,&tmp3[0]+i*nq01,1); // 2/(1-z2) for d/dxi_1
Vmath::Smul(nq01,gfac2[i],&tmp4[0]+i*nq01,1,&tmp5[0]+i*nq01,1); // 2/(1-z2) for d/dxi_2
}
// (1+z0)/(1-z2) for d/d eta_0
for(i = 0; i < nquad1*nquad2; ++i)
{
Vmath::Vmul(nquad0,&gfac0[0],1,
&tmp5[0]+i*nquad0,1,
&wsp[0]+i*nquad0,1);
}
Vmath::Vadd(nqtot, &tmp2[0], 1, &wsp[0], 1, &tmp2[0], 1);
// (1+z1)/(1-z2) for d/d eta_1
for(i = 0; i < nquad1*nquad2; ++i)
{
Vmath::Smul(nquad0,gfac1[i%nquad1],
&tmp5[0]+i*nquad0,1,
&tmp5[0]+i*nquad0,1);
}
Vmath::Vadd(nqtot, &tmp3[0], 1, &tmp5[0], 1, &tmp3[0], 1);
IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
m_base[1]->GetBdata (),
m_base[2]->GetBdata (),
tmp2,outarray,wsp,
false,true,true);
IProductWRTBase_SumFacKernel(m_base[0]->GetBdata (),
m_base[1]->GetDbdata(),
m_base[2]->GetBdata (),
tmp3,tmp6,wsp,
true,false,true);
Vmath::Vadd(m_ncoeffs, tmp6, 1, outarray, 1, outarray, 1);
IProductWRTBase_SumFacKernel(m_base[0]->GetBdata (),
m_base[1]->GetBdata (),
m_base[2]->GetDbdata(),
tmp4,tmp6,wsp,
true,true,false);
Vmath::Vadd(m_ncoeffs, tmp6, 1, outarray, 1, outarray, 1);
}
//---------------------------------------
// Evaluation functions
//---------------------------------------
......@@ -743,6 +945,36 @@ namespace Nektar
}
}
void PyrExp::v_SVVLaplacianFilter(
Array<OneD, NekDouble> &array,
const StdRegions::StdMatrixKey &mkey)
{
int nq = GetTotPoints();
// Calculate sqrt of the Jacobian
Array<OneD, const NekDouble> jac =
m_metricinfo->GetJac(GetPointsKeys());
Array<OneD, NekDouble> sqrt_jac(nq);
if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
Vmath::Vsqrt(nq,jac,1,sqrt_jac,1);
}
else
{
Vmath::Fill(nq,sqrt(jac[0]),sqrt_jac,1);
}
// Multiply array by sqrt(Jac)
Vmath::Vmul(nq,sqrt_jac,1,array,1,array,1);
// Apply std region filter
StdPyrExp::v_SVVLaplacianFilter( array, mkey);
// Divide by sqrt(Jac)
Vmath::Vdiv(nq,array,1,sqrt_jac,1,array,1);
}
//---------------------------------------
// Matrix creation functions
//---------------------------------------
......@@ -1027,7 +1259,6 @@ namespace Nektar
returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
}
}
return returnval;
......
......@@ -96,6 +96,18 @@ namespace Nektar
LOCAL_REGIONS_EXPORT virtual void v_IProductWRTBase(
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble>& outarray);
LOCAL_REGIONS_EXPORT virtual void v_IProductWRTBase_SumFac(
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble>& outarray,
bool multiplybyweights = true);
LOCAL_REGIONS_EXPORT void v_IProductWRTDerivBase(
const int dir,
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble>& outarray);
LOCAL_REGIONS_EXPORT void v_IProductWRTDerivBase_SumFac(
const int dir,
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble>& outarray);
//---------------------------------------
......@@ -131,6 +143,12 @@ namespace Nektar
Array<OneD, int> &outarray);
LOCAL_REGIONS_EXPORT void v_ComputeFaceNormal(const int face);
LOCAL_REGIONS_EXPORT virtual void v_SVVLaplacianFilter(
Array<OneD, NekDouble> &array,
const StdRegions::StdMatrixKey &mkey);
//---------------------------------------
// Matrix creation functions
//---------------------------------------
......
......@@ -3305,7 +3305,7 @@ namespace Nektar
returnval.push_back(bkey);
returnval.push_back(bkey);
const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha2Beta0);
const LibUtilities::PointsKey pkey1(nummodes+quadoffset, LibUtilities::eGaussRadauMAlpha2Beta0);
LibUtilities::BasisKey bkey1(LibUtilities::eModifiedPyr_C, nummodes, pkey1);
returnval.push_back(bkey1);
}
......
......@@ -420,8 +420,7 @@ namespace Nektar
int order0 = m_base[0]->GetNumModes();
int order1 = m_base[1]->GetNumModes();
Array<OneD, NekDouble> wsp(nquad1*nquad2*order0 +
nquad2*order0*order1);
Array<OneD, NekDouble> wsp(order0*nquad1*(nquad2 + order1));
if(multiplybyweights)
{
......@@ -1789,7 +1788,7 @@ namespace Nektar
const int Q = m_base[1]->GetNumModes()-1;
const int R = m_base[2]->GetNumModes()-1;
int i,j,l;
int i,l;
int cnt = 0;
// Traverse to q-r plane number I
......@@ -1869,5 +1868,139 @@ namespace Nektar
break;
}
}
void StdPyrExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
const StdMatrixKey &mkey)
{
// Generate an orthonogal expansion
int qa = m_base[0]->GetNumPoints();
int qb = m_base[1]->GetNumPoints();
int qc = m_base[2]->GetNumPoints();
int nmodes_a = m_base[0]->GetNumModes();
int nmodes_b = m_base[1]->GetNumModes();
int nmodes_c = m_base[2]->GetNumModes();
// Declare orthogonal basis.
LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A,nmodes_b,pb);
LibUtilities::BasisKey Bc(LibUtilities::eOrthoPyr_C,nmodes_c,pc);
StdPyrExp OrthoExp(Ba,Bb,Bc);
Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
int i,j,k,cnt = 0;
//SVV filter paramaters (how much added diffusion relative to physical one
// and fraction of modes from which you start applying this added diffusion)
//
NekDouble SvvDiffCoeff = mkey.GetConstFactor(StdRegions::eFactorSVVDiffCoeff);
NekDouble SVVCutOff = mkey.GetConstFactor(StdRegions::eFactorSVVCutoffRatio);
//Defining the cut of mode
int cutoff_a = (int) (SVVCutOff*nmodes_a);
int cutoff_b = (int) (SVVCutOff*nmodes_b);
int cutoff_c = (int) (SVVCutOff*nmodes_c);
//To avoid the fac[j] from blowing up
NekDouble epsilon = 1;
// project onto modal space.
OrthoExp.FwdTrans(array,orthocoeffs);
int nmodes = min(min(nmodes_a,nmodes_b),nmodes_c);
NekDouble cutoff = min(min(cutoff_a,cutoff_b),cutoff_c);
for(i = 0; i < nmodes_a; ++i)//P
{
for(j = 0; j < nmodes_b; ++j) //Q
{
int maxij = max(i,j);
for(k = 0; k < nmodes_c-maxij; ++k) //R
{
if(j + k >= cutoff || i + k >= cutoff)
{
orthocoeffs[cnt] *= (SvvDiffCoeff*exp(-(i+k-nmodes)*(i+k-nmodes)/((NekDouble)((i+k-cutoff+epsilon)*(i+k-cutoff+epsilon))))*exp(-(j-nmodes)*(j-nmodes)/((NekDouble)((j-cutoff+epsilon)*(j-cutoff+epsilon)))));
}
else
{
orthocoeffs[cnt] *= 0.0;
}
cnt++;
}
}
}
// backward transform to physical space
OrthoExp.BwdTrans(orthocoeffs,array);
}
void StdPyrExp::v_ReduceOrderCoeffs(
int numMin,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
{
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
int nqtot = nquad0*nquad1*nquad2;
int nmodes0 = m_base[0]->GetNumModes();
int nmodes1 = m_base[1]->GetNumModes();
int nmodes2 = m_base[2]->GetNumModes();
int numMax = nmodes0;
Array<OneD, NekDouble> coeff (m_ncoeffs);
Array<OneD, NekDouble> coeff_tmp1(m_ncoeffs, 0.0);
Array<OneD, NekDouble> phys_tmp (nqtot, 0.0);
Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
const LibUtilities::PointsKey Pkey2 = m_base[2]->GetPointsKey();
LibUtilities::BasisKey bortho0(
LibUtilities::eOrtho_A, nmodes0, Pkey0);
LibUtilities::BasisKey bortho1(
LibUtilities::eOrtho_A, nmodes1, Pkey1);
LibUtilities::BasisKey bortho2(
LibUtilities::eOrthoPyr_C, nmodes2, Pkey2);
int cnt = 0;
int u = 0;
int i = 0;
StdRegions::StdPyrExpSharedPtr OrthoPyrExp;
OrthoPyrExp = MemoryManager<StdRegions::StdPyrExp>
::AllocateSharedPtr(bortho0, bortho1, bortho2);
BwdTrans(inarray,phys_tmp);
OrthoPyrExp->FwdTrans(phys_tmp, coeff);
// filtering
for (u = 0; u < numMin; ++u)
{
for (i = 0; i < numMin; ++i)
{
int maxui = max(u,i);
Vmath::Vcopy(numMin - maxui, tmp = coeff + cnt, 1,
tmp2 = coeff_tmp1 + cnt, 1);
cnt += nmodes2 - maxui;
}
for (i = numMin; i < nmodes1; ++i)
{
int maxui = max(u,i);
cnt += numMax - maxui;
}
}
OrthoPyrExp->BwdTrans(coeff_tmp1, phys_tmp);
StdPyrExp::FwdTrans(phys_tmp, outarray);
}
}
}
......@@ -251,6 +251,16 @@ namespace Nektar
STD_REGIONS_EXPORT virtual DNekMatSharedPtr v_CreateStdMatrix(
const StdMatrixKey &mkey);
STD_REGIONS_EXPORT virtual void v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
const StdMatrixKey &mkey);
//---------------------------------------
// Method for applying sensors
//---------------------------------------
STD_REGIONS_EXPORT virtual void v_ReduceOrderCoeffs(
int numMin,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray);
private:
//---------------------------------------
// Private helper functions
......
......@@ -58,6 +58,7 @@ IF( NEKTAR_SOLVER_INCNAVIERSTOKES )
ADD_NEKTAR_TEST(Hex_channel_varP)
ADD_NEKTAR_TEST(Pyr_channel_m3)
ADD_NEKTAR_TEST(Pyr_channel_varP)
ADD_NEKTAR_TEST(Pyr_channel_SVV)
ADD_NEKTAR_TEST(Hex_channel_m6_nodalRestart)
ADD_NEKTAR_TEST_LENGTHY(Hex_channel_m8)
ADD_NEKTAR_TEST_LENGTHY(Hex_channel_m8_srhs)
......
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>3D channel flow, Pyramidic elements, using SVV</description>
<executable>IncNavierStokesSolver</executable>
<parameters>Pyr_channel_SVV.xml</parameters>
<files>
<file description="Session File">Pyr_channel_SVV.xml</file>