Commit 282cc75f authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'feature/dealiasing' of localhost:nektar

parents 94ff5684 c88af532
......@@ -186,27 +186,21 @@ int main(int argc, char *argv[])
// /////////////////////////////////////////////////////////
// /////////////////////////////////////////////////////////
// Interpolation //
// /////////////////////////////////////////////////////////
// /////////////////////////////////////////////////////////
// Generate a list of interpolating nodes
int nNodes = 2*nPts; // Number of interpolating nodes
boost::shared_ptr<Points<NekDouble> > nodes = PointsManager()[PointsKey(nNodes, pointsType)];
// const ptr<Points<NekDouble> > nodes = PointsManager()[PointsKey(nNodes, pointsType)];
//SharedArray<const NekDouble> zNode = nodes->GetZ();
Array<OneD, const NekDouble> zNode = nodes->GetZ();
// Get the interpolation matrix I
// Note that I is 2N rows by N columns
const Points<NekDouble>::MatrixSharedPtrType Iptr = points->GetI(nNodes,zNode);
// const Points<NekDouble>::MatrixSharedPtrType Iptr = points->GetI(nNodes, zNode);
const Points<NekDouble>::MatrixSharedPtrType Iptr = points->GetI(nNodes,zNode);
const NekMatrix<NekDouble> & I = *Iptr;
// Interpolate the data values in the y vector using the interpolation matrix I
// SharedArray<NekDouble> u(I.GetRows());
Array<OneD, NekDouble> u(I.GetRows());
for(int i = 0; i < int(I.GetRows()); ++i)
{
......@@ -216,7 +210,7 @@ int main(int argc, char *argv[])
u[i] += I(i,j) * y[j];
}
}
// Display the original samples
cout << setprecision(3);
cout << "\nOriginal data: \nx = ";
......@@ -282,7 +276,92 @@ int main(int argc, char *argv[])
}
///////////////////////////////////////////////////////////
// Galkerin Projection //
///////////////////////////////////////////////////////////
// Generate a list of projection nodes
Array<OneD, NekDouble> yNode(zNode.num_elements());
for(int i = 0; i < nNodes; ++i)
{
yNode[i] = func( zNode[i], nNodes, pointsType );
}
PointsKey key1(nNodes, pointsType);
// Note that I is 2N rows by N columns
const Points<NekDouble>::MatrixSharedPtrType GPptr = points->GetGalerkinProjection(key1);
const NekMatrix<NekDouble> & GP = *GPptr;
// Project the data values in the yNode vector using the projection matrix GP
for(int i = 0; i < int(GP.GetRows()); ++i)
{
u[i] = 0;
for(int j = 0; j < int(GP.GetColumns()); ++j)
{
u[i] += GP(i,j) * yNode[j];
}
}
// Display the original samples
cout << setprecision(3);
cout << "\n\n\n **** Galerkin Project ****";
cout << "\n\nResults of Galkerin Project with " << kPointsTypeStr[pointsType] << ":";
cout << "\nOriginal data: \nx = ";
for(int i = 0; i < nPts; ++i)
{
cout << setw(6) << z[i] << " ";
}
cout << "\ny = ";
for(int i = 0; i < nPts; ++i)
{
cout << setw(6) << y[i] << " ";
}
cout << "\nproject = ";
for(int i = 0; i < nPts; ++i)
{
cout << setw(6) << u[i] << " ";
}
// Display the pointwise error
cout << setprecision(1);
cout << "\nerror = ";
Linf = 0;
RMS = 0;
for(int i = 0; i < int(GP.GetRows()); ++i)
{
long double exact = func(z[i], nPts, pointsType);
long double error = exact - u[i];
Linf = max(Linf, fabs(error));
RMS += error*error;
long double epsilon = 1e-2;
if( fabs(exact) > epsilon )
{
error /= exact;
}
cout << setw(6) << error << " ";
}
RMS = sqrt(RMS) / int(GP.GetRows());
cout << setprecision(6);
cout << "\nLinf = " << setw(6) << Linf;
cout << "\nRMS = " << setw(6) << RMS << endl;
// Show the projection matrix
cout << "\nI = " << endl;
for(int i = 0; i < int(GP.GetRows()); ++i)
{
cout << " ";
for(int j = 0; j < int(GP.GetColumns()); ++j)
{
printf("% 5.3f ", GP(i,j));
}
cout << "" << endl;
}
// /////////////////////////////////////////////////////////
......
......@@ -69,6 +69,7 @@ SET(FoundationHeaders
./Foundations/NodalTriEvenlySpaced.h
./Foundations/NodalTetEvenlySpaced.h
./Foundations/NodalPrismEvenlySpaced.h
./Foundations/PhysGalerkinProject.h
)
SET(InterpreterHeaders
......@@ -167,6 +168,7 @@ SET(FoundationSources
./Foundations/NodalTriEvenlySpaced.cpp
./Foundations/NodalTetEvenlySpaced.cpp
./Foundations/NodalPrismEvenlySpaced.cpp
./Foundations/PhysGalerkinProject.cpp
)
SET(InterpreterSources
......
......@@ -654,6 +654,23 @@ 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
}
......
......@@ -94,6 +94,10 @@ namespace Nektar
bool UseNumMode = false,
TranspositionDir dir = eNoTrans);
LIB_UTILITIES_EXPORT void SetSpecVanVisc(Array<OneD, NekDouble> visc);
LIB_UTILITIES_EXPORT NekDouble GetSpecVanVisc(const int k);
protected:
CommSharedPtr m_hcomm;
......@@ -156,6 +160,9 @@ 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;
......
......@@ -307,6 +307,7 @@ namespace Nektar
return GetI(numpoints, xpoints);
}
const boost::shared_ptr<NekMatrix<NekDouble> > GaussPoints::GetI(const PointsKey &pkey)
{
ASSERTL0(pkey.GetPointsDim()==1, "Gauss Points can only interp to other 1d point distributions");
......@@ -391,6 +392,49 @@ namespace Nektar
}
return y;
}
const boost::shared_ptr<NekMatrix<NekDouble> > GaussPoints::GetGalerkinProjection(const PointsKey &pkey)
{
return m_GalerkinProjectionManager[pkey];
}
boost::shared_ptr< NekMatrix<NekDouble> > GaussPoints::CreateGPMatrix(const PointsKey &pkey)
{
boost::shared_ptr< NekMatrix<NekDouble> > returnval = CalculateGalerkinProjectionMatrix(pkey);
// Delegate to function below
return returnval;
}
boost::shared_ptr<NekMatrix<NekDouble> > GaussPoints::CalculateGalerkinProjectionMatrix(const PointsKey &pkey)
{
int numpointsfrom = pkey.GetNumPoints();
int numpointsto = GetNumPoints();
Array<OneD, const NekDouble> weightsfrom;
weightsfrom = PointsManager()[pkey]->GetW();
boost::shared_ptr< NekMatrix<NekDouble> > Interp = GetI(pkey);
Array<OneD, NekDouble> GalProj(numpointsfrom*numpointsto);
// set up inner product matrix and multiply by inverse of
// diagaonal mass matrix
for(int i = 0; i < numpointsto; ++i)
{
Vmath::Vmul(numpointsfrom,Interp->GetPtr().get() +i*numpointsfrom,1,
&weightsfrom[0],1,&GalProj[0] +i,numpointsto);
Vmath::Smul(numpointsfrom,1.0/m_weights[i],&GalProj[0]+i,numpointsto,
&GalProj[0]+i,numpointsto);
}
NekDouble* t = GalProj.data();
boost::shared_ptr< NekMatrix<NekDouble> > returnval(MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(numpointsto,numpointsfrom,t));
return returnval;
}
} // end of namespace LibUtilities
} // end of namespace Nektar
......
......@@ -62,6 +62,12 @@ namespace Nektar
LIB_UTILITIES_EXPORT const boost::shared_ptr<NekMatrix<NekDouble> > GetI(const Array<OneD, const NekDouble>& x);
LIB_UTILITIES_EXPORT const boost::shared_ptr<NekMatrix<NekDouble> > GetI(unsigned int numpoints, const Array<OneD, const NekDouble>& x);
LIB_UTILITIES_EXPORT boost::shared_ptr< NekMatrix<NekDouble> > CreateGPMatrix(const PointsKey &pkey);
LIB_UTILITIES_EXPORT const boost::shared_ptr<NekMatrix<NekDouble> > GetGalerkinProjection(const PointsKey &pkey);
GaussPoints(const PointsKey &pkey):PointsBaseType(pkey)
{
m_InterpManager.RegisterCreator(PointsKey(0, eGaussGaussLegendre),
......@@ -101,7 +107,46 @@ namespace Nektar
m_InterpManager.RegisterCreator(PointsKey(0, ePolyEvenlySpaced),
boost::bind(&GaussPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, eBoundaryLayerPoints),
boost::bind(&GaussPoints::CreateMatrix, this, _1)); }
boost::bind(&GaussPoints::CreateMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussGaussLegendre),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauMLegendre),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauPLegendre),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussLobattoLegendre),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussGaussChebyshev),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauMChebyshev),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauPChebyshev),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussLobattoChebyshev),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauMAlpha0Beta1),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauMAlpha0Beta2),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauMAlpha1Beta0),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauMAlpha2Beta0),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussKronrodLegendre),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauKronrodMLegendre),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussRadauKronrodMAlpha1Beta0),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eGaussLobattoKronrodLegendre),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eFourierEvenlySpaced),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, ePolyEvenlySpaced),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1));
m_GalerkinProjectionManager.RegisterCreator(PointsKey(0, eBoundaryLayerPoints),
boost::bind(&GaussPoints::CreateGPMatrix, this, _1)); }
private:
/// These should not be called. All creation is done
......@@ -115,6 +160,10 @@ namespace Nektar
void CalculateDerivMatrix();
void CalculateInterpMatrix(unsigned int npts, const Array<OneD, const NekDouble>& xpoints, Array<OneD, NekDouble>& interp);
boost::shared_ptr<NekMatrix<NekDouble> > CalculateGalerkinProjectionMatrix(const PointsKey &pkey);
/// functions used by the Kronrod points
NekDouble LagrangeInterpolant(NekDouble x, int npts, const Array<OneD, const NekDouble>& xpts, const Array<OneD, const NekDouble>& funcvals);
NekDouble LagrangePoly(NekDouble x, int pt, int npts, const Array<OneD, const NekDouble>& xpts);
......
///////////////////////////////////////////////////////////////////////////////
//
// File PhysGalerkinProject.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Definition of Physical Space Galerkin Projection methods
//
///////////////////////////////////////////////////////////////////////////////
#include <LibUtilities/Foundations/PhysGalerkinProject.h>
#include <LibUtilities/BasicUtils/Vmath.hpp>
#include <LibUtilities/BasicUtils/VmathArray.hpp>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <LibUtilities/LinearAlgebra/NekTypeDefs.hpp>
#include <LibUtilities/Foundations/Points.h>
#include <LibUtilities/Foundations/Basis.h>
namespace Nektar
{
namespace LibUtilities
{
// Physical Space Interpolation methods
// 1D Interpolation
void PhysGalerkinProject1D(const BasisKey &fbasis0,
const Array<OneD, const NekDouble>& from,
const BasisKey &tbasis0,
Array<OneD, NekDouble> &to)
{
PhysGalerkinProject1D(fbasis0.GetPointsKey(),from,tbasis0.GetPointsKey(),to);
}
void PhysGalerkinProject1D(const PointsKey &fpoints0,
const Array<OneD, const NekDouble>& from,
const PointsKey &tpoints0,
Array<OneD, NekDouble> &to)
{
if(fpoints0 == tpoints0) //check to see if the same
{
Vmath::Vcopy(fpoints0.GetNumPoints(),from,1,to,1);
}
else // interpolate
{
DNekMatSharedPtr GP0;
GP0 = PointsManager()[tpoints0]->GetGalerkinProjection(fpoints0);
NekVector<NekDouble> in(fpoints0.GetNumPoints(),from,eWrapper);
NekVector<NekDouble> out(tpoints0.GetNumPoints(),to,eWrapper);
GP0->Transpose();
out = (*GP0)*in;
}
}
void PhysGalerkinProject1D(const BasisKey &fbasis0,
const NekDouble *from,
const BasisKey &tbasis0,
NekDouble *to)
{
PhysGalerkinProject1D(fbasis0.GetPointsKey(),from,tbasis0.GetPointsKey(),to);
}
void PhysGalerkinProject1D(const PointsKey &fpoints0,
const NekDouble *from,
const PointsKey &tpoints0,
NekDouble *to)
{
if(fpoints0 == tpoints0) //check to see if the same
{
Vmath::Vcopy(fpoints0.GetNumPoints(),from,1,to,1);
}
else // interpolate
{
DNekMatSharedPtr GP0;
GP0 = PointsManager()[tpoints0]
->GetGalerkinProjection(fpoints0);
Blas::Dgemv('T', tpoints0.GetNumPoints(), fpoints0.GetNumPoints(),
1.0, GP0->GetPtr().get(), tpoints0.GetNumPoints(),
from, 1, 0.0, to, 1);
}
}
// 2D Interpolation
void PhysGalerkinProject2D(const BasisKey &fbasis0,
const BasisKey &fbasis1,
const Array<OneD, const NekDouble>& from,
const BasisKey &tbasis0,
const BasisKey &tbasis1,
Array<OneD, NekDouble> &to)
{
PhysGalerkinProject2D(fbasis0.GetPointsKey(),fbasis1.GetPointsKey(),from.data(),
tbasis0.GetPointsKey(),tbasis1.GetPointsKey(),to.data());
}
void PhysGalerkinProject2D(const PointsKey &fpoints0,
const PointsKey &fpoints1,
const Array<OneD, const NekDouble>& from,
const PointsKey &tpoints0,
const PointsKey &tpoints1,
Array<OneD, NekDouble> &to)
{
PhysGalerkinProject2D(fpoints0,fpoints1,from.data(),tpoints0,tpoints1,to.data());
}
void PhysGalerkinProject2D(const PointsKey &fpoints0,
const PointsKey &fpoints1,
const NekDouble *from,
const PointsKey &tpoints0,
const PointsKey &tpoints1,
NekDouble *to)
{
DNekMatSharedPtr GP0,GP1;
Array<OneD, NekDouble> wsp(tpoints1.GetNumPoints()*fpoints0.GetNumPoints()); // fnp0*tnp1
int fnp0 = fpoints0.GetNumPoints();
int fnp1 = fpoints1.GetNumPoints();
int tnp0 = tpoints0.GetNumPoints();
int tnp1 = tpoints1.GetNumPoints();
if(fpoints1 == tpoints1)
{
Vmath::Vcopy(fnp0*tnp1,from,1,wsp.get(),1);
}
else
{
GP1 = PointsManager()[tpoints1]->GetGalerkinProjection(fpoints1);
Blas::Dgemm('N', 'T', fnp0, tnp1, fnp1, 1.0, from, fnp0,
GP1->GetPtr().get(), tnp1, 0.0, wsp.get(), fnp0);
}
if(fpoints0 == tpoints0)
{
Vmath::Vcopy(tnp0*tnp1,wsp.get(),1,to,1);
}
else
{
GP0 = PointsManager()[tpoints0]->GetGalerkinProjection(fpoints0);
Blas::Dgemm('N', 'N', tnp0, tnp1, fnp0, 1.0,
GP0->GetPtr().get(),
tnp0, wsp.get(), fnp0, 0.0, to, tnp0);
}
}
} // end of namespace
} // end of namespace
///////////////////////////////////////////////////////////////////////////////
//
// File PhysGalerkinProject.h
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Definition of Physcal space (1D tensor based) Galerkin Projection
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_LIB_UTILTIES_FOUNDATIONS_PHYSGALERKIN_H
#define NEKTAR_LIB_UTILTIES_FOUNDATIONS_PHYSGALERKIN_H
#include <LibUtilities/Foundations/FoundationsFwd.hpp>
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
#include <LibUtilities/LibUtilitiesDeclspec.h>
namespace Nektar { template <typename Dim, typename DataType> class Array; }
namespace Nektar
{
namespace LibUtilities
{
// Physical Space Interpolation methods
LIB_UTILITIES_EXPORT void PhysGalerkinProject1D(const BasisKey &fbasis0,
const Array<OneD, const NekDouble>& from,
const BasisKey &tbasis0,
Array<OneD, NekDouble> &to);
LIB_UTILITIES_EXPORT void PhysGalerkinProject1D(const PointsKey &fpoints0,
const Array<OneD, const NekDouble>& from,
const PointsKey &tpoints0,
Array<OneD, NekDouble> &to);
LIB_UTILITIES_EXPORT void PhysGalerkinProject1D(const BasisKey &fbasis0,
const NekDouble *from,
const BasisKey &tbasis0,
NekDouble *to);
LIB_UTILITIES_EXPORT void PhysGalerkinProject1D(const PointsKey &fpoints0,
const NekDouble *from,
const PointsKey &tpoints0,
NekDouble *to);
// 2D PhysGalkerinProjectolation
LIB_UTILITIES_EXPORT void PhysGalerkinProject2D(const BasisKey &fbasis0,
const BasisKey &fbasis1,
const Array<OneD, const NekDouble>& from,
const BasisKey &tbasis0,
const BasisKey &tbasis1,
Array<OneD, NekDouble> &to);
LIB_UTILITIES_EXPORT void PhysGalerkinProject2D(const PointsKey &fpoints0,