Commit d3948e5b authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'feature/SpechpSVV' of localhost:nektar

parents c12b0358 b5a56d79
......@@ -649,23 +649,6 @@ namespace Nektar
}
}
void Transposition::SetSpecVanVisc(Array<OneD, NekDouble> visc)
{
m_specVanVisc = visc;
}
NekDouble Transposition::GetSpecVanVisc(const int k)
{
NekDouble returnval = 0.0;
if(m_specVanVisc.num_elements())
{
returnval = m_specVanVisc[k];
}
return returnval;
}
// TODO: Impelement 2D and 3D transposition routines
}
......
......@@ -160,9 +160,6 @@ namespace Nektar
/// MPI_Alltoallv offset map of send/recv buffer in global vector.
Array<OneD,int> m_OffsetMap;
/// Spectral vanishing Viscosity coefficient for stabilisation
Array<OneD, NekDouble> m_specVanVisc;
};
typedef boost::shared_ptr<Transposition> TranspositionSharedPtr;
......
......@@ -1413,7 +1413,6 @@ namespace Nektar
return returnval;
}
DNekMatSharedPtr QuadExp::v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
{
LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
......@@ -1522,11 +1521,11 @@ namespace Nektar
case StdRegions::eLaplacian:
{
if( (m_metricinfo->GetGtype() == SpatialDomains::eDeformed) ||
(mkey.GetNVarCoeff() > 0) )
(mkey.GetNVarCoeff() > 0)||(mkey.ConstFactorExists(StdRegions::eFactorSVVCutoffRatio)))
{
NekDouble one = 1.0;
DNekMatSharedPtr mat = GenMatrix(mkey);
returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
}
else
......@@ -1901,11 +1900,11 @@ namespace Nektar
}
void QuadExp::v_LaplacianMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
void QuadExp::v_LaplacianMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
{
if(mkey.GetNVarCoeff() == 0)
if(mkey.GetNVarCoeff() == 0 &&!mkey.ConstFactorExists(StdRegions::eFactorSVVCutoffRatio))
{
// This implementation is only valid when there are no
// coefficients associated to the Laplacian operator
......@@ -1917,18 +1916,18 @@ namespace Nektar
int nmodes0 = m_base[0]->GetNumModes();
int nmodes1 = m_base[1]->GetNumModes();
int wspsize = max(max(max(nqtot,m_ncoeffs),nquad1*nmodes0),nquad0*nmodes1);
const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
const Array<TwoD, const NekDouble>& metric = m_metricinfo->GetLaplacianMetrics();
// Allocate temporary storage
Array<OneD,NekDouble> wsp0(3*wspsize);
Array<OneD,NekDouble> wsp1(wsp0+wspsize);
Array<OneD,NekDouble> wsp2(wsp0+2*wspsize);
if(!(m_base[0]->Collocation() && m_base[1]->Collocation()))
{
// LAPLACIAN MATRIX OPERATION
......@@ -1942,7 +1941,8 @@ namespace Nektar
{
StdExpansion2D::PhysTensorDeriv(inarray,wsp1,wsp2);
}
// wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
// wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
// where g0, g1 and g2 are the metric terms set up in the GeomFactors class
......@@ -1960,12 +1960,12 @@ namespace Nektar
Vmath::Vmul(nqtot,&metric[0][0],1,&wsp1[0],1,&wsp0[0],1);
Vmath::Vmul(nqtot,&metric[2][0],1,&wsp2[0],1,&wsp2[0],1);
}
// outarray = m = (D_xi1 * B)^T * k
// wsp1 = n = (D_xi2 * B)^T * l
IProductWRTBase_SumFacKernel(dbase0,base1,wsp0,outarray,wsp1,false,true);
IProductWRTBase_SumFacKernel(base0,dbase1,wsp2,wsp1, wsp0,true,false);
// outarray = outarray + wsp1
// = L * u_hat
Vmath::Vadd(m_ncoeffs,wsp1.get(),1,outarray.get(),1,outarray.get(),1);
......@@ -1978,12 +1978,12 @@ namespace Nektar
int nmodes0 = m_base[0]->GetNumModes();
int nmodes1 = m_base[1]->GetNumModes();
int wspsize = max(max(max(nqtot,m_ncoeffs),nquad1*nmodes0),nquad0*nmodes1);
const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
// Allocate temporary storage
Array<OneD,NekDouble> wsp0(6*wspsize);
Array<OneD,NekDouble> wsp1(wsp0+wspsize);
......@@ -1991,7 +1991,7 @@ namespace Nektar
Array<OneD,NekDouble> wsp3(wsp0+3*wspsize);
Array<OneD,NekDouble> wsp4(wsp0+4*wspsize);
Array<OneD,NekDouble> wsp5(wsp0+5*wspsize);
if(!(m_base[0]->Collocation() && m_base[1]->Collocation()))
{
// LAPLACIAN MATRIX OPERATION
......@@ -2005,7 +2005,7 @@ namespace Nektar
{
StdExpansion2D::PhysTensorDeriv(inarray,wsp1,wsp2);
}
// wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
// wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
// where g0, g1 and g2 are the metric terms set up in the GeomFactors class
......@@ -2020,7 +2020,7 @@ namespace Nektar
Vmath::Vvtvvtp(nqtot,&gmat[0][0],1,&gmat[1][0],1,&gmat[2][0],1,&gmat[3][0],1,&wsp4[0],1);
// wsp5 = g1*g1 + g3*g3;
Vmath::Vvtvvtp(nqtot,&gmat[1][0],1,&gmat[1][0],1,&gmat[3][0],1,&gmat[3][0],1,&wsp5[0],1);
// If 3D coordinates, tag on extra terms
if(dim == 3)
{
......@@ -2031,7 +2031,7 @@ namespace Nektar
// wsp5 += g5*g5
Vmath::Vvtvp(nqtot,&gmat[5][0],1,&gmat[5][0],1,&wsp5[0],1,&wsp5[0],1);
}
Vmath::Vvtvvtp(nqtot,&wsp3[0],1,&wsp1[0],1,&wsp4[0],1,&wsp2[0],1,&wsp0[0],1);
Vmath::Vvtvvtp(nqtot,&wsp4[0],1,&wsp1[0],1,&wsp5[0],1,&wsp2[0],1,&wsp2[0],1);
}
......@@ -2040,14 +2040,14 @@ namespace Nektar
NekDouble g0 = gmat[0][0]*gmat[0][0] + gmat[2][0]*gmat[2][0];
NekDouble g1 = gmat[0][0]*gmat[1][0] + gmat[2][0]*gmat[3][0];
NekDouble g2 = gmat[1][0]*gmat[1][0] + gmat[3][0]*gmat[3][0];
if(dim == 3)
{
g0 += gmat[4][0]*gmat[4][0];
g1 += gmat[4][0]*gmat[5][0];
g2 += gmat[5][0]*gmat[5][0];
}
if(fabs(g1) < NekConstants::kGeomFactorsTol)
{
Vmath::Smul(nqtot,g0,&wsp1[0],1,&wsp0[0],1);
......@@ -2059,15 +2059,15 @@ namespace Nektar
Vmath::Svtsvtp(nqtot,g1,&wsp1[0],1,g2,&wsp2[0],1,&wsp2[0],1);
}
}
MultiplyByQuadratureMetric(wsp0,wsp0);
MultiplyByQuadratureMetric(wsp2,wsp2);
// outarray = m = (D_xi1 * B)^T * k
// wsp1 = n = (D_xi2 * B)^T * l
IProductWRTBase_SumFacKernel(dbase0,base1,wsp0,outarray,wsp1,false,true);
IProductWRTBase_SumFacKernel(base0,dbase1,wsp2,wsp1, wsp0,true,false);
// outarray = outarray + wsp1
// = L * u_hat
Vmath::Vadd(m_ncoeffs,wsp1.get(),1,outarray.get(),1,outarray.get(),1);
......@@ -2080,6 +2080,7 @@ namespace Nektar
}
void QuadExp::v_HelmholtzMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
......
......@@ -1327,7 +1327,7 @@ namespace Nektar
case StdRegions::eLaplacian:
{
if( (m_metricinfo->GetGtype() == SpatialDomains::eDeformed) ||
(mkey.GetNVarCoeff() > 0) )
(mkey.GetNVarCoeff() > 0)||(mkey.ConstFactorExists(StdRegions::eFactorSVVCutoffRatio)))
{
NekDouble one = 1.0;
DNekMatSharedPtr mat = GenMatrix(mkey);
......@@ -1715,7 +1715,8 @@ namespace Nektar
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
{
if(mkey.GetNVarCoeff() == 0)
if(mkey.GetNVarCoeff() == 0 && !mkey.ConstFactorExists(StdRegions::eFactorSVVCutoffRatio))
{
// This implementation is only valid when there are no coefficients
// associated to the Laplacian operator
......@@ -1746,6 +1747,7 @@ namespace Nektar
BwdTrans_SumFacKernel(base0,base1,inarray,wsp0,wsp1);
StdExpansion2D::PhysTensorDeriv(wsp0,wsp1,wsp2);
// wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
// wsp2 = l = g1 * wsp1 + g2 * wsp2 = g1 * du_dxi1 + g2 * du_dxi2
// where g0, g1 and g2 are the metric terms set up in the GeomFactors class
......@@ -1796,6 +1798,7 @@ namespace Nektar
BwdTrans_SumFacKernel(base0,base1,inarray,wsp0,wsp1);
StdExpansion2D::PhysTensorDeriv(wsp0,wsp1,wsp2);
// wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
// wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
const Array<TwoD, const NekDouble>& gmat = m_metricinfo->GetGmat();
......@@ -1881,6 +1884,7 @@ namespace Nektar
// wsp2 = l = wsp4 * wsp1 + wsp5 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
Vmath::Vvtvvtp(nqtot,&wsp3[0],1,&wsp1[0],1,&wsp4[0],1,&wsp2[0],1,&wsp0[0],1);
Vmath::Vvtvvtp(nqtot,&wsp4[0],1,&wsp1[0],1,&wsp5[0],1,&wsp2[0],1,&wsp2[0],1);
// substep 4: multiply by jacobian and quadrature weights
MultiplyByQuadratureMetric(wsp0,wsp0);
MultiplyByQuadratureMetric(wsp2,wsp2);
......
......@@ -233,7 +233,7 @@ namespace Nektar
beta = 2*M_PI*(m_transposition->GetK(n))/m_lhom;
new_factors = factors;
// add in Homogeneous Fourier direction and SVV if turned on
new_factors[StdRegions::eFactorLambda] += beta*beta*(1+m_transposition->GetSpecVanVisc(n));
new_factors[StdRegions::eFactorLambda] += beta*beta*(1+GetSpecVanVisc(n));
m_planes[n]->HelmSolve(fce + cnt,
e_out = outarray + cnt1,
......
......@@ -244,7 +244,7 @@ namespace Nektar
beta = 2*M_PI*(m_transposition->GetK(n))/m_lhom;
new_factors = factors;
// add in Homogeneous Fourier direction and SVV if turned on
new_factors[StdRegions::eFactorLambda] += beta*beta*(1+m_transposition->GetSpecVanVisc(n));
new_factors[StdRegions::eFactorLambda] += beta*beta*(1+GetSpecVanVisc(n));
m_planes[n]->HelmSolve(fce + cnt,
e_out = outarray + cnt1,
flags, new_factors, varcoeff, dirForcing);
......
......@@ -401,16 +401,16 @@ namespace Nektar
/// Set the \a i th coefficiient in \a m_coeffs to value \a val
inline void SetCoeff(int i, NekDouble val);
/// Set the coefficiient in \a m_coeffs to value \a val (0D Exapnsion)
/// Set the coefficiient in \a m_coeffs to value \a val (0D Exapnsion)
inline void SetCoeff(NekDouble val);
/// Set the physical value in \a m_coeffs to value \a val (0D Exapnsion)
/// Set the physical value in \a m_coeffs to value \a val (0D Exapnsion)
inline void SetPhys(NekDouble val);
inline const SpatialDomains::VertexComponentSharedPtr &GetGeom(void) const;
inline const SpatialDomains::VertexComponentSharedPtr &GetVertex(void) const;
inline const SpatialDomains::VertexComponentSharedPtr &GetGeom(void) const;
inline const SpatialDomains::VertexComponentSharedPtr &GetVertex(void) const;
/// Set the \a i th coefficiient in #m_coeffs to value \a val
inline void SetCoeffs(int i, NekDouble val);
......@@ -468,13 +468,20 @@ namespace Nektar
return v_L2();
}
/// This function calculates the energy associated with each one of the modes
/// of a 3D homogeneous nD expansion
/// This function calculates the energy associated with
/// each one of the modesof a 3D homogeneous nD expansion
Array<OneD, const NekDouble> HomogeneousEnergy (void)
{
return v_HomogeneousEnergy();
}
/// This function sets the Spectral Vanishing Viscosity
/// in homogeneous1D expansion.
void SetHomo1DSpecVanVisc(Array<OneD, NekDouble> visc)
{
v_SetHomo1DSpecVanVisc(visc);
}
/// This function returns a vector containing the wave
/// numbers in z-direction associated
/// with the 3D homogenous expansion. Required if a
......@@ -525,6 +532,7 @@ namespace Nektar
/// This function returns the number of elements in the expansion.
inline int GetExpSize(void);
/// This function returns the number of elements in the
/// expansion which may be different for a homogeoenous extended
/// expansionp.
......@@ -1212,6 +1220,14 @@ namespace Nektar
return LibUtilities::NullBasisSharedPtr;
}
// wrapper function to set viscosity for Homo1D expansion
virtual void v_SetHomo1DSpecVanVisc(Array<OneD, NekDouble> visc)
{
ASSERTL0(false,
"This method is not defined or valid for this class type");
}
virtual boost::shared_ptr<ExpList> &v_GetPlane(int n);
};
......
......@@ -158,7 +158,25 @@ namespace Nektar
DNekBlkMatSharedPtr GetHomogeneous1DBlockMatrix(Homogeneous1DMatType mattype, CoeffState coeffstate = eLocal) const;
NekDouble GetSpecVanVisc(const int k)
{
NekDouble returnval = 0.0;
if(m_specVanVisc.num_elements())
{
returnval = m_specVanVisc[k];
}
return returnval;
}
// virtual functions
virtual void v_SetHomo1DSpecVanVisc(Array<OneD, NekDouble> visc)
{
m_specVanVisc = visc;
}
virtual int v_GetNumElmts(void)
{
return m_planes[0]->GetExpSize();
......@@ -189,8 +207,7 @@ namespace Nektar
virtual void v_IProductWRTBase_IterPerExp(const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray);
virtual std::vector<SpatialDomains::FieldDefinitionsSharedPtr> v_GetFieldDefinitions(void);
virtual void v_GetFieldDefinitions(std::vector<SpatialDomains::FieldDefinitionsSharedPtr> &fielddef);
......@@ -259,6 +276,9 @@ namespace Nektar
//Padding operations variables
bool m_dealiasing;
int m_padsize;
/// Spectral vanishing Viscosity coefficient for stabilisation
Array<OneD, NekDouble> m_specVanVisc;
};
inline void ExpListHomogeneous1D::HomogeneousFwdTrans(const Array<OneD, const NekDouble> &inarray,
......
This diff is collapsed.
......@@ -1058,6 +1058,13 @@ namespace Nektar
v_LaplacianMatrixOp(inarray,outarray,mkey);
}
void SVVLaplacianFilter(Array<OneD,NekDouble> &array,
const StdMatrixKey &mkey)
{
v_SVVLaplacianFilter(array,mkey);
}
void LaplacianMatrixOp(const int k1, const int k2,
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
......@@ -1795,6 +1802,9 @@ namespace Nektar
Array<OneD,NekDouble> &outarray,
const StdMatrixKey &mkey);
STD_REGIONS_EXPORT virtual void v_SVVLaplacianFilter(Array<OneD,NekDouble> &array,
const StdMatrixKey &mkey);
STD_REGIONS_EXPORT virtual void v_LaplacianMatrixOp(const int k1, const int k2,
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
......
......@@ -109,7 +109,7 @@ namespace Nektar
return false;
}
if(lhs.m_ncoeffs < rhs.m_ncoeffs)
if(lhs.m_ncoeffs < rhs.m_ncoeffs) // probably do not need this check since checking the m_base below?
{
return true;
}
......@@ -144,7 +144,7 @@ namespace Nektar
{
ConstFactorMap::const_iterator x, y;
for(x = lhs.m_factors.begin(), y = rhs.m_factors.begin();
x != lhs.m_factors.end(); ++x, ++y)
x != lhs.m_factors.end(); ++x, ++y)
{
if (x->second < y->second)
{
......
......@@ -128,6 +128,17 @@ namespace Nektar
return x->second;
}
inline bool ConstFactorExists(const ConstFactorType& factor) const
{
ConstFactorMap::const_iterator x = m_factors.find(factor);
if(x != m_factors.end())
{
return true;
}
return false;
}
inline const ConstFactorMap& GetConstFactors() const
{
return m_factors;
......
......@@ -1293,15 +1293,15 @@ namespace Nektar
DNekMat &Imass = *GetStdMatrix(imasskey);
(*Mat) = Imass*Iprod;
}
break;
default:
}
break;
default:
{
Mat = StdExpansion::CreateGeneralMatrix(mkey);
}
break;
}
return Mat;
}
......@@ -1398,6 +1398,82 @@ namespace Nektar
}
}
void StdQuadExp::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 qtot = qa*qb;
int nmodes_a = m_base[0]->GetNumModes();
int nmodes_b = m_base[1]->GetNumModes();
// Declare orthogonal basis.
LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A,nmodes_b,pb);
StdQuadExp OrthoExp(Ba,Bb);
Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
int j,k;
int cuttoff = (int) (mkey.GetConstFactor(eFactorSVVCutoffRatio)*min(nmodes_a,nmodes_b));
NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
// project onto modal space.
OrthoExp.FwdTrans(array,orthocoeffs);
#if 0 // Filter just linear space
int nmodes = min(nmodes_a,nmodes_b);
// apply SVV filter.
for(j = 0; j < nmodes_a; ++j)
{
for(k = 0; k < nmodes_b; ++k)
{
if(j + k >= cuttoff)
{
orthocoeffs[j*nmodes_b+k] *= (1.0+SvvDiffCoeff*exp(-(j+k-nmodes)*(j+k-nmodes)/((NekDouble)((j+k-cuttoff+1)*(j+k-cuttoff+1)))));
}
}
}
#else // Filter just bilinear space
int nmodes = max(nmodes_a,nmodes_b);
Array<OneD, NekDouble> fac(nmodes,0.0);
for(j = cuttoff; j < nmodes; ++j)
{
fac[j] = exp(-(j-nmodes)*(j-nmodes)/((NekDouble) (j-cuttoff+1.0)*(j-cuttoff+1.0)));
}
for(j = cuttoff; j < nmodes_a; ++j)
{
for(k = 0; k < cuttoff; ++k)
{
orthocoeffs[j*nmodes_b+k] *= (1.0+SvvDiffCoeff*fac[j]);
}
for(k = cuttoff; k < nmodes_b; ++k)
{
orthocoeffs[j*nmodes_b+k] *= (1.0+SvvDiffCoeff*fac[j]*fac[k]);
}
}
for(j = 0; j < nmodes_a; ++j)
{
for(k = cuttoff; k < nmodes_b; ++k)
{
orthocoeffs[j*nmodes_b+k] *= (1.0+SvvDiffCoeff*fac[k]);
}
}
#endif
// backward transform to physical space
OrthoExp.BwdTrans(orthocoeffs,array);
}
void StdQuadExp::v_HelmholtzMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
......
......@@ -247,6 +247,9 @@ namespace Nektar
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdMatrixKey &mkey);
STD_REGIONS_EXPORT virtual void v_SVVLaplacianFilter(
Array<OneD, NekDouble> &array,
const StdMatrixKey &mkey);
STD_REGIONS_EXPORT virtual void v_HelmholtzMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
......
......@@ -177,7 +177,8 @@ namespace Nektar
"HybridDGLamToQ1",
"HybridDGLamToQ2",
"HybridDGLamToU",
"FwdTrans"
"FwdTrans",
"Preconditioner",
};
enum VarCoeffType
......@@ -215,13 +216,17 @@ namespace Nektar
{
eFactorLambda,
eFactorTau,
eFactorTime
eFactorTime,
eFactorSVVCutoffRatio,
eFactorSVVDiffCoeff
};
const char* const ConstFactorTypeMap[] = {
"FactorLambda",
"FactorTau",
"FactorTime"
"FactorTime",
"FactorSVVCutoffRatio",