Commit 74d54580 authored by Spencer Sherwin's avatar Spencer Sherwin
Browse files

Added in Collection methods for BwdTrans and PhysDeriv and setup unit test

parent 4ffa7606
......@@ -918,7 +918,7 @@ class BwdTrans_SumFac_Prism : public Operator
{
Blas::Dgemm('N', 'N', m_nquad2, m_numElmt, m_nmodes2-i,
1.0, m_base2.get()+mode*m_nquad2, m_nquad2,
&input[0]+mode1, totmodes, 0.0,
input.get()+mode1, totmodes, 0.0,
&wsp[j*m_nquad2*m_numElmt*m_nmodes0+ cnt],
m_nquad2);
mode1 += m_nmodes2-i;
......@@ -1018,5 +1018,168 @@ OperatorKey BwdTrans_SumFac_Prism::m_type = GetOperatorFactory().
OperatorKey(ePrism, eBwdTrans, eSumFac,false),
BwdTrans_SumFac_Prism::create, "BwdTrans_SumFac_Prism");
/**
* @brief Backward transform operator using sum-factorisation (Pyr)
*/
class BwdTrans_SumFac_Pyr : public Operator
{
public:
OPERATOR_CREATE(BwdTrans_SumFac_Pyr)
virtual ~BwdTrans_SumFac_Pyr()
{
}
virtual void operator()(
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &output1,
Array<OneD, NekDouble> &output2,
Array<OneD, NekDouble> &wsp)
{
ASSERTL1(wsp.num_elements() == m_wspSize,
"Incorrect workspace size");
// Assign second half of workspace for 2nd DGEMM operation.
int totmodes = m_stdExp->GetNcoeffs();
Array<OneD, NekDouble> wsp2
= wsp + m_nmodes0*m_nmodes1*m_nquad2*m_numElmt;
Vmath::Zero(m_nmodes0*m_nmodes1*m_nquad2*m_numElmt, wsp, 1);
int i = 0;
int j = 0;
int mode = 0;
int mode1 = 0;
int cnt = 0;
for (i = mode = mode1 = 0; i < m_nmodes0; ++i)
{
for (j = 0; j < m_nmodes1; ++j, ++cnt)
{
int ijmax = max(i,j);
Blas::Dgemm('N', 'N', m_nquad2, m_numElmt, m_nmodes2-ijmax,
1.0, m_base2.get()+mode*m_nquad2, m_nquad2,
input.get()+mode1, totmodes, 0.0,
wsp.get() + cnt*m_nquad2*m_numElmt, m_nquad2);
mode += m_nmodes2-ijmax;
mode1 += m_nmodes2-ijmax;
}
//increment mode in case order1!=order2
for(j = m_nmodes1; j < m_nmodes2-i; ++j)
{
int ijmax = max(i,j);
mode += m_nmodes2-ijmax;
}
}
// vertex mode - currently (1+c)/2 x (1-b)/2 x (1-a)/2
// component is evaluated
if(m_sortTopVertex)
{
for(i = 0; i < m_numElmt; ++i)
{
// top singular vertex
// (1+c)/2 x (1+b)/2 x (1-a)/2 component
Blas::Daxpy(m_nquad2, input[1+i*totmodes],
m_base2.get() + m_nquad2, 1,
&wsp[m_nquad2*m_numElmt] + i*m_nquad2, 1);
// top singular vertex
// (1+c)/2 x (1-b)/2 x (1+a)/2 component
Blas::Daxpy(m_nquad2, input[1+i*totmodes],
m_base2.get() + m_nquad2, 1,
&wsp[m_nmodes1*m_nquad2*m_numElmt]
+ i*m_nquad2, 1);
// top singular vertex
// (1+c)/2 x (1+b)/2 x (1+a)/2 component
Blas::Daxpy(m_nquad2, input[1+i*totmodes],
m_base2.get() + m_nquad2, 1,
&wsp[(m_nmodes1+1)*m_nquad2*m_numElmt]
+ i*m_nquad2, 1);
}
}
// Perform summation over '1' direction
mode = 0;
for(i = 0; i < m_nmodes0; ++i)
{
Blas::Dgemm('N', 'T', m_nquad1, m_nquad2*m_numElmt, m_nmodes1,
1.0, m_base1.get(), m_nquad1,
wsp.get() + mode*m_nquad2*m_numElmt,
m_nquad2*m_numElmt,
0.0, wsp2.get() + i*m_nquad1*m_nquad2*m_numElmt,
m_nquad1);
mode += m_nmodes1;
}
// Perform summation over '0' direction
Blas::Dgemm('N', 'T', m_nquad0, m_nquad1*m_nquad2*m_numElmt,
m_nmodes0, 1.0, m_base0.get(), m_nquad0,
wsp2.get(), m_nquad1*m_nquad2*m_numElmt,
0.0, output.get(), m_nquad0);
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nquad1;
const int m_nquad2;
const int m_nmodes0;
const int m_nmodes1;
const int m_nmodes2;
Array<OneD, const NekDouble> m_base0;
Array<OneD, const NekDouble> m_base1;
Array<OneD, const NekDouble> m_base2;
bool m_sortTopVertex;
private:
BwdTrans_SumFac_Pyr(
vector<StdRegions::StdExpansionSharedPtr> pCollExp,
CoalescedGeomDataSharedPtr pGeomData)
: Operator (pCollExp, pGeomData),
m_nquad0 (m_stdExp->GetNumPoints(0)),
m_nquad1 (m_stdExp->GetNumPoints(1)),
m_nquad2 (m_stdExp->GetNumPoints(2)),
m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
m_nmodes2 (m_stdExp->GetBasisNumModes(2)),
m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
m_base2 (m_stdExp->GetBasis(2)->GetBdata())
{
m_wspSize = m_numElmt*m_nmodes0*m_nquad2*(m_nmodes1 + m_nquad1);
if(m_stdExp->GetBasis(0)->GetBasisType()
== LibUtilities::eModified_A)
{
m_sortTopVertex = true;
}
else
{
m_sortTopVertex = false;
}
}
};
/// Factory initialisation for the BwdTrans_SumFac_Pyr operator
OperatorKey BwdTrans_SumFac_Pyr::m_type = GetOperatorFactory().
RegisterCreatorFunction(
OperatorKey(ePyramid, eBwdTrans, eSumFac,false),
BwdTrans_SumFac_Pyr::create, "BwdTrans_SumFac_Pyr");
}
}
......@@ -1537,5 +1537,228 @@ OperatorKey PhysDeriv_SumFac_Prism::m_typeArr[] = {
};
/**
* @brief Phys deriv operator using sum-factorisation (Pyramid)
*/
class PhysDeriv_SumFac_Pyr : public Operator
{
public:
OPERATOR_CREATE(PhysDeriv_SumFac_Pyr)
virtual ~PhysDeriv_SumFac_Pyr()
{
}
virtual void operator()(
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output0,
Array<OneD, NekDouble> &output1,
Array<OneD, NekDouble> &output2,
Array<OneD, NekDouble> &wsp)
{
int nPhys = m_stdExp->GetTotPoints();
int ntot = m_numElmt*nPhys;
Array<OneD, NekDouble> tmp0,tmp1,tmp2;
Array<OneD, Array<OneD, NekDouble> > Diff(3);
Array<OneD, Array<OneD, NekDouble> > out(3);
out[0] = output0; out[1] = output1; out[2] = output2;
for(int i = 0; i < m_dim; ++i)
{
Diff[i] = wsp + i*ntot;
}
// dEta0
Blas::Dgemm('N','N', m_nquad0,m_nquad1*m_nquad2*m_numElmt,
m_nquad0,1.0, m_Deriv0,m_nquad0,&input[0],
m_nquad0,0.0,&Diff[0][0],m_nquad0);
int cnt = 0;
for(int i = 0; i < m_numElmt; ++i)
{
// dEta 1
for (int j = 0; j < m_nquad2; ++j)
{
Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1,
1.0, &input[i*nPhys+j*m_nquad0*m_nquad1],
m_nquad0, m_Deriv1, m_nquad1, 0.0,
&Diff[1][i*nPhys+j*m_nquad0*m_nquad1],
m_nquad0);
}
// dEta 2
Blas::Dgemm('N','T',m_nquad0*m_nquad1,m_nquad2,m_nquad2,
1.0, &input[i*nPhys],m_nquad0*m_nquad1,
m_Deriv2,m_nquad2, 0.0,&Diff[2][i*nPhys],
m_nquad0*m_nquad1);
// dxi0 = 2/(1-eta_2) d Eta_0
Vmath::Vmul(nPhys,&m_fac0[0],1,Diff[0].get()+cnt,1,
Diff[0].get()+cnt,1);
// dxi1 = 2/(1-eta_2) d Eta_1
Vmath::Vmul(nPhys,&m_fac0[0],1,Diff[1].get()+cnt,1,
Diff[1].get()+cnt,1);
// dxi2 = (1+eta0)/(1-eta_2) d Eta_0 + d/dEta2;
Vmath::Vvtvp(nPhys,&m_fac1[0],1,Diff[0].get()+cnt,1,
Diff[2].get()+cnt,1,Diff[2].get()+cnt,1);
// dxi2 += (1+eta1)/(1-eta_2) d Eta_1
Vmath::Vvtvp(nPhys,&m_fac2[0],1,Diff[1].get()+cnt,1,
Diff[2].get()+cnt,1,Diff[2].get()+cnt,1);
cnt += nPhys;
}
// calculate full derivative
for(int i = 0; i < m_coordim; ++i)
{
Vmath::Vmul(ntot,m_derivFac[i*m_dim],1,Diff[0],1,out[i],1);
for(int j = 1; j < m_dim; ++j)
{
Vmath::Vvtvp (ntot, m_derivFac[i*m_dim+j], 1,
Diff[j], 1, out[i], 1, out[i], 1);
}
}
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
int nPhys = m_stdExp->GetTotPoints();
int ntot = m_numElmt*nPhys;
Array<OneD, NekDouble> tmp0,tmp1,tmp2;
Array<OneD, Array<OneD, NekDouble> > Diff(3);
for(int i = 0; i < m_dim; ++i)
{
Diff[i] = wsp + i*ntot;
}
// dEta0
Blas::Dgemm('N','N', m_nquad0,m_nquad1*m_nquad2*m_numElmt,
m_nquad0,1.0, m_Deriv0,m_nquad0,&input[0],
m_nquad0,0.0,&Diff[0][0],m_nquad0);
int cnt = 0;
for(int i = 0; i < m_numElmt; ++i)
{
// dEta 1
for (int j = 0; j < m_nquad2; ++j)
{
Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1,
1.0, &input[i*nPhys+j*m_nquad0*m_nquad1],
m_nquad0, m_Deriv1, m_nquad1, 0.0,
&Diff[1][i*nPhys+j*m_nquad0*m_nquad1],
m_nquad0);
}
// dEta 2
Blas::Dgemm('N','T',m_nquad0*m_nquad1,m_nquad2,m_nquad2,
1.0, &input[i*nPhys],m_nquad0*m_nquad1,
m_Deriv2,m_nquad2, 0.0,&Diff[2][i*nPhys],
m_nquad0*m_nquad1);
// dxi0 = 2/(1-eta_2) d Eta_0
Vmath::Vmul(nPhys,&m_fac0[0],1,Diff[0].get()+cnt,1,
Diff[0].get()+cnt,1);
// dxi1 = 2/(1-eta_2) d Eta_1
Vmath::Vmul(nPhys,&m_fac0[0],1,Diff[1].get()+cnt,1,
Diff[1].get()+cnt,1);
// dxi2 = (1+eta0)/(1-eta_2) d Eta_0 + d/dEta2;
Vmath::Vvtvp(nPhys,&m_fac1[0],1,Diff[0].get()+cnt,1,
Diff[2].get()+cnt,1,Diff[2].get()+cnt,1);
// dxi2 = (1+eta1)/(1-eta_2) d Eta_1 + d/dEta2;
Vmath::Vvtvp(nPhys,&m_fac2[0],1,Diff[1].get()+cnt,1,
Diff[2].get()+cnt,1,Diff[2].get()+cnt,1);
cnt += nPhys;
}
// calculate full derivative
Vmath::Vmul(ntot,m_derivFac[dir*m_dim],1,Diff[0],1,output,1);
for(int j = 1; j < m_dim; ++j)
{
Vmath::Vvtvp (ntot, m_derivFac[dir*m_dim+j], 1,
Diff[j], 1, output, 1, output, 1);
}
}
protected:
Array<TwoD, const NekDouble> m_derivFac;
int m_dim;
int m_coordim;
const int m_nquad0;
const int m_nquad1;
const int m_nquad2;
NekDouble *m_Deriv0;
NekDouble *m_Deriv1;
NekDouble *m_Deriv2;
Array<OneD, NekDouble> m_fac0;
Array<OneD, NekDouble> m_fac1;
Array<OneD, NekDouble> m_fac2;
private:
PhysDeriv_SumFac_Pyr(
vector<StdRegions::StdExpansionSharedPtr> pCollExp,
CoalescedGeomDataSharedPtr pGeomData)
: Operator(pCollExp, pGeomData),
m_nquad0 (m_stdExp->GetNumPoints(0)),
m_nquad1 (m_stdExp->GetNumPoints(1)),
m_nquad2 (m_stdExp->GetNumPoints(2))
{
LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
m_dim = PtsKey.size();
m_coordim = m_stdExp->GetCoordim();
m_derivFac = pGeomData->GetDerivFactors(pCollExp);
const Array<OneD, const NekDouble>& z0
= m_stdExp->GetBasis(0)->GetZ();
const Array<OneD, const NekDouble>& z1
= m_stdExp->GetBasis(1)->GetZ();
const Array<OneD, const NekDouble>& z2
= m_stdExp->GetBasis(2)->GetZ();
m_fac0 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
m_fac1 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
m_fac2 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
for (int i = 0; i < m_nquad0; ++i)
{
for(int j = 0; j < m_nquad1; ++j)
{
for(int k = 0; k < m_nquad2; ++k)
{
m_fac0[i+j*m_nquad0 + k*m_nquad0*m_nquad1] =
2.0/(1-z2[k]);
m_fac1[i+j*m_nquad0 + k*m_nquad0*m_nquad1] =
0.5*(1+z0[i]);
m_fac2[i+j*m_nquad0 + k*m_nquad0*m_nquad1] =
0.5*(1+z1[j]);
}
}
}
m_Deriv0 = &((m_stdExp->GetBasis(0)->GetD())->GetPtr())[0];
m_Deriv1 = &((m_stdExp->GetBasis(1)->GetD())->GetPtr())[0];
m_Deriv2 = &((m_stdExp->GetBasis(2)->GetD())->GetPtr())[0];
m_wspSize = 3*m_nquad0*m_nquad1*m_nquad2*m_numElmt;
}
};
/// Factory initialisation for the PhysDeriv_SumFac_Pyr operators
OperatorKey PhysDeriv_SumFac_Pyr::m_typeArr[] = {
GetOperatorFactory().RegisterCreatorFunction(
OperatorKey(ePyramid, ePhysDeriv, eSumFac, false),
PhysDeriv_SumFac_Pyr::create, "PhysDeriv_SumFac_Pyr")
};
}
}
......@@ -124,20 +124,50 @@ namespace Nektar
eta_z = m_base[2]->GetZ();
int i, j, k, n;
if (out_dxi1.num_elements() > 0)
{
for (k = 0, n = 0; k < Qz; ++k)
{
NekDouble fac = 2.0/(1.0 - eta_z[k]);
for (j = 0; j < Qy; ++j)
{
for (i = 0; i < Qx; ++i, ++n)
{
out_dxi1[n] = fac * dEta_bar1[n];
}
}
}
}
if (out_dxi2.num_elements() > 0)
{
for (k = 0, n = 0; k < Qz; ++k)
{
NekDouble fac = 2.0/(1.0 - eta_z[k]);
for (j = 0; j < Qy; ++j)
{
for (i = 0; i < Qx; ++i, ++n)
{
out_dxi2[n] = fac * dXi2[n];
}
}
}
}
for (k = 0, n = 0; k < Qz; ++k)
if (out_dxi3.num_elements() > 0)
{
for (j = 0; j < Qy; ++j)
for (k = 0, n = 0; k < Qz; ++k)
{
for (i = 0; i < Qx; ++i, ++n)
NekDouble fac = 1.0/(1.0 - eta_z[k]);
for (j = 0; j < Qy; ++j)
{
if (out_dxi1.num_elements() > 0)
out_dxi1[n] = 2.0/(1.0 - eta_z[k]) * dEta_bar1[n];
if (out_dxi2.num_elements() > 0)
out_dxi2[n] = 2.0/(1.0 - eta_z[k]) * dXi2[n];
if (out_dxi3.num_elements() > 0)
out_dxi3[n] = (1.0+eta_x[i])/(1.0-eta_z[k])*dEta_bar1[n] +
(1.0+eta_y[j])/(1.0-eta_z[k])*dXi2[n] + dEta3[n];
NekDouble fac1 = (1.0+eta_y[j]);
for (i = 0; i < Qx; ++i, ++n)
{
out_dxi3[n] = (1.0+eta_x[i])*fac*dEta_bar1[n] +
fac1*fac*dXi2[n] + dEta3[n];
}
}
}
}
......
......@@ -3,6 +3,7 @@ SET(Sources
TestHexCollection.cpp
TestQuadCollection.cpp
TestPrismCollection.cpp
TestPyrCollection.cpp
TestSegCollection.cpp
TestTetCollection.cpp
TestTriCollection.cpp
......
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment