Commit 77ea6bba authored by Dave Moxey's avatar Dave Moxey
Browse files

Adding new advection and Riemann solver classes.


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@3910 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent 4584e859
......@@ -247,7 +247,7 @@ using namespace Nektar;
}
/// \brief vvtvm (vector times vector plus vector): z = w*x - y
/// \brief vvtvm (vector times vector minus vector): z = w*x - y
template<class T> void Vvtvm(int n, const Array<OneD,const T> &w, const int incw, const Array<OneD,const T> &x, const int incx, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
{
ASSERTL1(n*incw <= w.num_elements()+w.GetOffset(),"Array out of bounds");
......@@ -258,7 +258,26 @@ using namespace Nektar;
Vvtvm(n,&w[0],incw,&x[0],incx,&y[0],incy,&z[0],incz);
}
/// \brief vvtvvtp (vector times vector plus vector times vector): z = v*w + y*z
template<class T> void Vvtvvtp (
int n,
const Array<OneD,const T> &v, int incv,
const Array<OneD,const T> &w, int incw,
const Array<OneD,const T> &x, int incx,
const Array<OneD,const T> &y, int incy,
Array<OneD, T> &z, int incz)
{
ASSERTL1(n*incv <= v.num_elements()+v.GetOffset(),"Array out of bounds");
ASSERTL1(n*incw <= w.num_elements()+w.GetOffset(),"Array out of bounds");
ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),"Array out of bounds");
Vvtvvtp(n,&v[0],incv,&w[0],incw,&x[0],incx,&y[0],incy,&z[0],incz);
}
/// \brief svtsvtp (scalar times vector plus scalar times vector): z = alpha*x + beta*y
template<class T> void Svtsvtp(int n, const T alpha, const Array<OneD,const T> &x, const int incx, const T beta, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
{
ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
......@@ -269,6 +288,8 @@ using namespace Nektar;
}
/************ Misc routine from Veclib (and extras) ************/
/// \brief Gather vector z[i] = x[y[i]]
......
///////////////////////////////////////////////////////////////////////////////
//
// File: Advection.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: Abstract base class for advection.
//
///////////////////////////////////////////////////////////////////////////////
#include <SolverUtils/Advection.h>
namespace Nektar
{
namespace SolverUtils
{
AdvectionFactory& GetAdvectionFactory()
{
typedef Loki::SingletonHolder<AdvectionFactory,
Loki::CreateUsingNew,
Loki::NoDestroy > Type;
return Type::Instance();
}
void Advection::Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
v_Advect(nConvectiveFields, fields, advVel, inarray, outarray);
}
}
}
///////////////////////////////////////////////////////////////////////////////
//
// File: Advection.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: Abstract base class for advection.
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_SOLVERUTILS_ADVECTION
#define NEKTAR_SOLVERUTILS_ADVECTION
#include <string>
#include <boost/function.hpp>
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <MultiRegions/ExpList.h>
#include <SolverUtils/RiemannSolvers/RiemannSolver.h>
namespace Nektar
{
namespace SolverUtils
{
typedef boost::function<void (
const int,
const Array<OneD, Array<OneD, NekDouble> >&,
Array<OneD, Array<OneD, NekDouble> >&)> AdvectionFluxVecCB;
class Advection
{
public:
void Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
template<typename FuncPointerT, typename ObjectPointerT>
inline void SetFluxVector(FuncPointerT func, ObjectPointerT obj)
{
m_fluxVector = boost::bind(func, obj, _1, _2, _3);
}
inline void SetRiemannSolver(RiemannSolverSharedPtr riemann)
{
m_riemann = riemann;
}
protected:
virtual void v_Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray) = 0;
AdvectionFluxVecCB m_fluxVector;
RiemannSolverSharedPtr m_riemann;
};
/// A shared pointer to an EquationSystem object
typedef boost::shared_ptr<Advection> AdvectionSharedPtr;
/// Datatype of the NekFactory used to instantiate classes derived
/// from the Advection class.
typedef LibUtilities::NekFactory<std::string, Advection>
AdvectionFactory;
AdvectionFactory& GetAdvectionFactory();
}
}
#endif
///////////////////////////////////////////////////////////////////////////////
//
// File: AdvectionNonConservative.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: Non-conservative advection class.
//
///////////////////////////////////////////////////////////////////////////////
#include <SolverUtils/AdvectionNonConservative.h>
namespace Nektar
{
namespace SolverUtils
{
std::string AdvectionNonConservative::type = GetAdvectionFactory().
RegisterCreatorFunction("NonConservative", AdvectionNonConservative::create);
AdvectionNonConservative::AdvectionNonConservative()
{
}
void AdvectionNonConservative::v_Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
int nDim = advVel.num_elements();
int nPointsTot = fields[0]->GetNpoints();
Array<OneD, NekDouble> grad0,grad1,grad2;
grad0 = Array<OneD, NekDouble> (nPointsTot);
if (nDim > 1)
{
grad1 = Array<OneD,NekDouble>(nPointsTot);
}
if (nDim > 2)
{
grad2 = Array<OneD,NekDouble>(nPointsTot);
}
for (int i = 0; i < nConvectiveFields; ++i)
{
// Evaluate V\cdot Grad(u)
switch(nDim)
{
case 1:
fields[0]->PhysDeriv(inarray[i], grad0);
Vmath::Vmul(nPointsTot,grad0,1,advVel[0],1,outarray[i],1);
break;
case 2:
fields[0]->PhysDeriv(inarray[i],grad0,grad1);
Vmath::Vmul (nPointsTot,grad0,1,advVel[0],1,outarray[i],1);
Vmath::Vvtvp(nPointsTot,grad1,1,advVel[1],1,outarray[i],1,outarray[i],1);
break;
case 3:
fields[0]->PhysDeriv(inarray[i],grad0,grad1,grad2);
Vmath::Vmul (nPointsTot,grad0,1,advVel[0],1,outarray[i],1);
Vmath::Vvtvvtp(nPointsTot,grad1,1,advVel[1],1,grad2,1,advVel[2],1,
outarray[i],1);
break;
default:
ASSERTL0(false,"dimension unknown");
}
}
}
}
}
///////////////////////////////////////////////////////////////////////////////
//
// File: AdvectionNonConservative.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: Non-conservative advection class.
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_SOLVERUTILS_ADVECTIONNONCONSERVATIVE
#define NEKTAR_SOLVERUTILS_ADVECTIONNONCONSERVATIVE
#include <SolverUtils/Advection.h>
namespace Nektar
{
namespace SolverUtils
{
class AdvectionNonConservative : public Advection
{
public:
static AdvectionSharedPtr create()
{
return AdvectionSharedPtr(new AdvectionNonConservative());
}
static std::string type;
protected:
AdvectionNonConservative();
virtual void v_Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
};
}
}
#endif
///////////////////////////////////////////////////////////////////////////////
//
// File: AdvectionWeakDG.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: Weak DG advection class.
//
///////////////////////////////////////////////////////////////////////////////
#include <SolverUtils/AdvectionWeakDG.h>
namespace Nektar
{
namespace SolverUtils
{
std::string AdvectionWeakDG::type = GetAdvectionFactory().
RegisterCreatorFunction("WeakDG", AdvectionWeakDG::create);
AdvectionWeakDG::AdvectionWeakDG()
{
}
void AdvectionWeakDG::v_Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
int i, j;
int nVelDim = fields[0]->GetCoordim(0);
int nPointsTot = fields[0]->GetTotPoints();
int nCoeffs = fields[0]->GetNcoeffs();
int nTracePointsTot = fields[0]->GetTrace()->GetTotPoints();
Array<OneD, Array<OneD, NekDouble> > fluxvector(nVelDim);
Array<OneD, Array<OneD, NekDouble> > tmp (nConvectiveFields);
ASSERTL1(m_riemann,
"Riemann solver must be provided for AdvectionWeakDG.");
for(i = 0; i < nVelDim; ++i)
{
fluxvector[i] = Array<OneD, NekDouble>(nPointsTot);
}
// Get the advection part (without numerical flux)
for(i = 0; i < nConvectiveFields; ++i)
{
tmp[i] = Array<OneD,NekDouble>(nCoeffs, 0.0);
// Get the ith component of the flux vector in (physical space)
m_fluxVector(i, inarray, fluxvector);
for (j = 0; j < nVelDim; ++j)
{
fields[i]->IProductWRTDerivBase(j, fluxvector[j], outarray[i]);
Vmath::Vadd(nCoeffs, outarray[i], 1, tmp[i], 1, tmp[i], 1);
}
}
// Store forwards/backwards space along trace space.
Array<OneD, Array<OneD, NekDouble> > Fwd (nConvectiveFields);
Array<OneD, Array<OneD, NekDouble> > Bwd (nConvectiveFields);
Array<OneD, Array<OneD, NekDouble> > numflux(nConvectiveFields);
for(i = 0; i < nConvectiveFields; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePointsTot);
Bwd[i] = Array<OneD, NekDouble>(nTracePointsTot);
numflux[i] = Array<OneD, NekDouble>(nTracePointsTot);
fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
}
m_riemann->Solve(Fwd, Bwd, numflux);
// Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
for(i = 0; i < nConvectiveFields; ++i)
{
Vmath::Neg (nCoeffs, tmp[i], 1);
fields[i]->AddTraceIntegral (numflux[i], tmp[i]);
fields[i]->MultiplyByElmtInvMass(tmp[i], tmp[i]);
fields[i]->BwdTrans (tmp[i], outarray[i]);
}
}
}
}
///////////////////////////////////////////////////////////////////////////////
//
// File: AdvectionWeakDG.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: Upwind Riemann solver.
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_SOLVERUTILS_ADVECTIONWEAKDG
#define NEKTAR_SOLVERUTILS_ADVECTIONWEAKDG
#include <SolverUtils/Advection.h>
namespace Nektar
{
namespace SolverUtils
{
class AdvectionWeakDG : public Advection
{
public:
static AdvectionSharedPtr create()
{
return AdvectionSharedPtr(new AdvectionWeakDG());
}
static std::string type;
protected:
AdvectionWeakDG();
Array<OneD, Array<OneD, NekDouble> > m_traceNormals;
virtual void v_Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
};
}
}
#endif
SET(SOLVER_UTILS_SOURCES
Advection.cpp
AdvectionNonConservative.cpp
AdvectionWeakDG.cpp
Driver.cpp
DriverArnoldi.cpp
DriverModifiedArnoldi.cpp
......@@ -8,10 +11,15 @@ SET(SOLVER_UTILS_SOURCES
Filters/FilterCheckpoint.cpp
Filters/FilterHistoryPoints.cpp
Filters/FilterThresholdMax.cpp
RiemannSolvers/RiemannSolver.cpp
RiemannSolvers/UpwindSolver.cpp
UnsteadySystem.cpp
)
SET(SOLVER_UTILS_HEADERS
Advection.h
AdvectionNonConservative.h
AdvectionWeakDG.h
Driver.h
DriverArnoldi.h
DriverModifiedArnoldi.h
......@@ -21,6 +29,8 @@ SET(SOLVER_UTILS_HEADERS
Filters/FilterCheckpoint.h
Filters/FilterHistoryPoints.h
Filters/FilterThresholdMax.h
RiemannSolvers/RiemannSolver.h
RiemannSolvers/UpwindSolver.h
SolverUtils.hpp
UnsteadySystem.h
)
......
///////////////////////////////////////////////////////////////////////////////
//
// File: RiemannSolver.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License