Commit 9371a8ef authored by Pavel Burovskiy's avatar Pavel Burovskiy
Browse files

Merge branch 'master' into fix/cg-svtvp

parents e7f2af74 cab0917c
......@@ -1159,7 +1159,7 @@ namespace Nektar
* direction X2.
* @param numericalFlux Riemann flux in the physical space.
* @param divCFlux Divergence of the corrective flux for 2D
* Problems.
* problems.
*
* \todo: Switch on shapes eventually here.
*/
......
......@@ -6,6 +6,8 @@ SET(SOLVER_UTILS_SOURCES
Diffusion/Diffusion.cpp
Diffusion/DiffusionLDG.cpp
Diffusion/DiffusionLDGNS.cpp
Diffusion/DiffusionLFR.cpp
Diffusion/DiffusionLFRNS.cpp
Driver.cpp
DriverArnoldi.cpp
DriverModifiedArnoldi.cpp
......@@ -18,6 +20,7 @@ SET(SOLVER_UTILS_SOURCES
Filters/FilterThresholdMax.cpp
RiemannSolvers/RiemannSolver.cpp
RiemannSolvers/UpwindSolver.cpp
RiemannSolvers/UpwindLDGSolver.cpp
UnsteadySystem.cpp
)
......@@ -29,6 +32,8 @@ SET(SOLVER_UTILS_HEADERS
Diffusion/Diffusion.h
Diffusion/DiffusionLDG.h
Diffusion/DiffusionLDGNS.h
Diffusion/DiffusionLFR.h
Diffusion/DiffusionLFRNS.h
Driver.h
DriverArnoldi.h
DriverModifiedArnoldi.h
......@@ -41,6 +46,7 @@ SET(SOLVER_UTILS_HEADERS
Filters/FilterThresholdMax.h
RiemannSolvers/RiemannSolver.h
RiemannSolvers/UpwindSolver.h
RiemannSolvers/UpwindLDGSolver.h
SolverUtils.hpp
SolverUtilsDeclspec.h
UnsteadySystem.h
......
......@@ -48,9 +48,10 @@ namespace Nektar
}
void Diffusion::InitObject(
const LibUtilities::SessionReaderSharedPtr pSession)
const LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
v_InitObject(pSession);
v_InitObject(pSession, pFields);
}
void Diffusion::Diffuse(
......
......@@ -58,7 +58,6 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> >&)> DiffusionFluxVecCB;
typedef boost::function<void (
const int,
const Array<OneD, Array<OneD, NekDouble> >&,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&)>
......@@ -68,7 +67,8 @@ namespace Nektar
{
public:
SOLVER_UTILS_EXPORT void InitObject(
LibUtilities::SessionReaderSharedPtr pSession);
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
SOLVER_UTILS_EXPORT void Diffuse(
const int nConvectiveFields,
......@@ -85,7 +85,7 @@ namespace Nektar
template<typename FuncPointerT, typename ObjectPointerT>
void SetFluxVectorNS(FuncPointerT func, ObjectPointerT obj)
{
m_fluxVectorNS = boost::bind(func, obj, _1, _2, _3, _4);
m_fluxVectorNS = boost::bind(func, obj, _1, _2, _3);
}
inline void SetRiemannSolver(RiemannSolverSharedPtr riemann)
......@@ -95,7 +95,8 @@ namespace Nektar
protected:
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession)
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
};
......
......@@ -34,6 +34,8 @@
///////////////////////////////////////////////////////////////////////////////
#include <SolverUtils/Diffusion/DiffusionLDG.h>
#include <iostream>
#include <iomanip>
namespace Nektar
{
......@@ -47,7 +49,8 @@ namespace Nektar
}
void DiffusionLDG::v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession)
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
m_session = pSession;
}
......@@ -58,6 +61,7 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
//cout<<setprecision(16);
int i, j, k;
int nDim = fields[0]->GetCoordim(0);
int nPts = fields[0]->GetTotPoints();
......@@ -72,17 +76,23 @@ namespace Nektar
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > flux (nDim);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > qfield(nDim);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > qfieldStd(nDim);
for (j = 0; j < nDim; ++j)
{
qfield[j] =
Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
qfieldStd[j] =
Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
flux[j] =
Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
for (i = 0; i < nConvectiveFields; ++i)
{
qfield[j][i] = Array<OneD, NekDouble>(nPts, 0.0);
qfieldStd[j][i] = Array<OneD, NekDouble>(nPts, 0.0);
flux[j][i] = Array<OneD, NekDouble>(nTracePts, 0.0);
}
}
......@@ -159,9 +169,15 @@ namespace Nektar
}
fields[0]->GetTrace()->GetNormals(m_traceNormals);
// Get the sign of (v \cdot n), v = an arbitrary vector
// Get the normal velocity Vn
for(i = 0; i < nDim; ++i)
{
Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
Vn, 1, Vn, 1);
}
// Evaulate upwind flux:
// Get the sign of (v \cdot n), v = an arbitrary vector
// Evaluate upwind flux:
// uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
for (j = 0; j < nDim; ++j)
{
......@@ -178,7 +194,7 @@ namespace Nektar
// edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
// edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
fields[i]->GetTrace()->Upwind(m_traceNormals[j],
fields[i]->GetTrace()->Upwind(/*m_traceNormals[j]*/Vn,
Fwd, Bwd, fluxtemp);
// Imposing weak boundary condition with flux
......@@ -192,8 +208,7 @@ namespace Nektar
if(fields[0]->GetBndCondExpansions().num_elements())
{
v_WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp,
time);
v_WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
}
// if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}),
......@@ -220,47 +235,31 @@ namespace Nektar
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const int var,
const Array<OneD, const NekDouble> &ufield,
Array<OneD, NekDouble> &penaltyflux,
NekDouble time)
Array<OneD, NekDouble> &penaltyflux)
{
int i, j, e, npoints, id1, id2;
int i, j, e, id1, id2;
// Number of boundary regions
int Nfps, numBDEdge;
int cnt = 0;
int nbnd = fields[var]->GetBndCondExpansions().num_elements();
int nDim = fields[0]->GetCoordim(0);
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
int nBndEdgePts, nBndEdges;
int cnt = 0;
int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
int nDim = fields[0]->GetCoordim(0);
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
Array<OneD, NekDouble > uplus(nTracePts);
fields[var]->ExtractTracePhys(ufield, uplus);
for (i = 0; i < nbnd; ++i)
for (i = 0; i < nBndRegions; ++i)
{
// Number of boundary expansion related to that region
numBDEdge = fields[var]->
nBndEdges = fields[var]->
GetBndCondExpansions()[i]->GetExpSize();
// Evaluate boundary values g_D or g_N from input files
LibUtilities::EquationSharedPtr ifunc =
m_session->GetFunction("InitialConditions", 0);
npoints = fields[var]->
GetBndCondExpansions()[i]->GetNpoints();
Array<OneD,NekDouble> BDphysics(npoints);
Array<OneD,NekDouble> x0(npoints, 0.0);
Array<OneD,NekDouble> x1(npoints, 0.0);
Array<OneD,NekDouble> x2(npoints, 0.0);
fields[var]->GetBndCondExpansions()[i]->GetCoords(x0, x1, x2);
ifunc->Evaluate(x0, x1, x2, time, BDphysics);
// Weakly impose boundary conditions by modifying flux values
for (e = 0; e < numBDEdge ; ++e)
for (e = 0; e < nBndEdges ; ++e)
{
// Number of points on the expansion
Nfps = fields[var]->
nBndEdgePts = fields[var]->
GetBndCondExpansions()[i]->GetExp(e)->GetNumPoints(0);
id1 = fields[var]->
......@@ -274,18 +273,17 @@ namespace Nektar
if (fields[var]->GetBndConditions()[i]->
GetBoundaryConditionType() == SpatialDomains::eDirichlet)
{
Vmath::Vcopy(Nfps,
Vmath::Vcopy(nBndEdgePts,
&(fields[var]->
GetBndCondExpansions()[i]->
UpdatePhys())[id1], 1,
GetPhys())[id1], 1,
&penaltyflux[id2], 1);
}
// For Neumann boundary condition: uflux = u+
else if ((fields[var]->GetBndConditions()[i])->
GetBoundaryConditionType() == SpatialDomains::eNeumann)
{
Vmath::Vcopy(Nfps,
Vmath::Vcopy(nBndEdgePts,
&uplus[id2], 1,
&penaltyflux[id2], 1);
}
......@@ -305,9 +303,8 @@ namespace Nektar
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
int nvariables = fields.num_elements();
int nDim = qfield.num_elements();
NekDouble time = 0.0;
NekDouble C11 = 1.0;
NekDouble C11 = 0.0;
Array<OneD, NekDouble > Fwd(nTracePts);
Array<OneD, NekDouble > Bwd(nTracePts);
Array<OneD, NekDouble > Vn (nTracePts, 0.0);
......@@ -326,6 +323,13 @@ namespace Nektar
}
fields[0]->GetTrace()->GetNormals(m_traceNormals);
// Get the normal velocity Vn
for(i = 0; i < nDim; ++i)
{
Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
Vn, 1, Vn, 1);
}
// Evaulate upwind flux:
// qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
for (i = 0; i < nvariables; ++i)
......@@ -334,7 +338,7 @@ namespace Nektar
for (j = 0; j < nDim; ++j)
{
// Compute Fwd and Bwd value of ufield of jth direction
fields[i]->GetFwdBwdTracePhys(qfield[j][i],qFwd,qBwd);
fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);
// if Vn >= 0, flux = uFwd, i.e.,
// edge::eForward, if V*n>=0 <=> V*n_F>=0, pick
......@@ -348,9 +352,9 @@ namespace Nektar
// edge::eBackward, if V*n<0 <=> V*n_B>=0, pick
// qflux = qFwd = q+
fields[i]->GetTrace()->Upwind(m_traceNormals[j],
qBwd, qFwd,
qfluxtemp);
fields[i]->GetTrace()->Upwind(/*m_traceNormals[j]*/Vn,
qBwd, qFwd,
qfluxtemp);
Vmath::Vmul(nTracePts,
m_traceNormals[j], 1,
......@@ -377,12 +381,9 @@ namespace Nektar
// Imposing weak boundary condition with flux
if (fields[0]->GetBndCondExpansions().num_elements())
{
v_WeakPenaltyforVector(fields,
i, j,
qfield[j][i],
qfluxtemp,
C11,
time);
v_WeakPenaltyforVector(fields, i, j,
qfield[j][i],
qfluxtemp, C11);
}
// q_hat \cdot n = (q_xi \cdot n_xi) or (q_eta \cdot n_eta)
......@@ -409,14 +410,13 @@ namespace Nektar
const int dir,
const Array<OneD, const NekDouble> &qfield,
Array<OneD, NekDouble> &penaltyflux,
NekDouble C11,
NekDouble time)
NekDouble C11)
{
int i, j, e, npoints, id1, id2;
int numBDEdge, Nfps;
int nbnd = fields[var]->GetBndCondExpansions().num_elements();
int nDim = fields[0]->GetCoordim(0);
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
int i, j, e, id1, id2;
int nBndEdges, nBndEdgePts;
int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
int nDim = fields[0]->GetCoordim(0);
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
Array<OneD, NekDouble > uterm(nTracePts);
Array<OneD, NekDouble > qtemp(nTracePts);
......@@ -432,30 +432,15 @@ namespace Nektar
fields[var]->ExtractTracePhys(qfield, qtemp);
for (i = 0; i < nbnd; ++i)
for (i = 0; i < nBndRegions; ++i)
{
numBDEdge = fields[var]->
nBndEdges = fields[var]->
GetBndCondExpansions()[i]->GetExpSize();
// Evaluate boundary values g_D or g_N from input files
LibUtilities::EquationSharedPtr ifunc =
m_session->GetFunction("InitialConditions", 0);
npoints = fields[var]->
GetBndCondExpansions()[i]->GetNpoints();
Array<OneD,NekDouble> BDphysics(npoints);
Array<OneD,NekDouble> x0(npoints, 0.0);
Array<OneD,NekDouble> x1(npoints, 0.0);
Array<OneD,NekDouble> x2(npoints, 0.0);
fields[var]->GetBndCondExpansions()[i]->GetCoords(x0, x1, x2);
ifunc->Evaluate(x0, x1, x2, time, BDphysics);
// Weakly impose boundary conditions by modifying flux values
for (e = 0; e < numBDEdge ; ++e)
for (e = 0; e < nBndEdges ; ++e)
{
Nfps = fields[var]->
nBndEdgePts = fields[var]->
GetBndCondExpansions()[i]->GetExp(e)->GetNumPoints(0);
id1 = fields[var]->
......@@ -470,7 +455,7 @@ namespace Nektar
if(fields[var]->GetBndConditions()[i]->
GetBoundaryConditionType() == SpatialDomains::eDirichlet)
{
Vmath::Vmul(Nfps,
Vmath::Vmul(nBndEdgePts,
&m_traceNormals[dir][id2], 1,
&qtemp[id2], 1,
&penaltyflux[id2], 1);
......@@ -479,11 +464,11 @@ namespace Nektar
else if((fields[var]->GetBndConditions()[i])->
GetBoundaryConditionType() == SpatialDomains::eNeumann)
{
Vmath::Vmul(Nfps,
Vmath::Vmul(nBndEdgePts,
&m_traceNormals[dir][id2], 1,
&(fields[var]->
GetBndCondExpansions()[i]->
UpdatePhys())[id1], 1,
GetPhys())[id1], 1,
&penaltyflux[id2], 1);
}
}
......
......@@ -60,7 +60,8 @@ namespace Nektar
LibUtilities::SessionReaderSharedPtr m_session;
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession);
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_Diffuse(
const int nConvective,
......@@ -77,8 +78,7 @@ namespace Nektar
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const int var,
const Array<OneD, const NekDouble> &ufield,
Array<OneD, NekDouble> &penaltyflux,
NekDouble time);
Array<OneD, NekDouble> &penaltyflux);
virtual void v_NumFluxforVector(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
......@@ -92,8 +92,7 @@ namespace Nektar
const int dir,
const Array<OneD, const NekDouble> &qfield,
Array<OneD, NekDouble> &penaltyflux,
NekDouble C11,
NekDouble time);
NekDouble C11);
};
}
}
......
......@@ -55,11 +55,21 @@ namespace Nektar
protected:
DiffusionLDGNS();
Array<OneD, Array<OneD, NekDouble> > m_traceNormals;
LibUtilities::SessionReaderSharedPtr m_session;
Array<OneD, Array<OneD, NekDouble> > m_traceNormals;
LibUtilities::SessionReaderSharedPtr m_session;
NekDouble m_gamma;
NekDouble m_gasConstant;
NekDouble m_Twall;
std::string m_ViscosityType;
NekDouble m_mu;
NekDouble m_thermalConductivity;
NekDouble m_rhoInf;
NekDouble m_pInf;
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession);
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_Diffuse(
const int nConvective,
......@@ -67,34 +77,32 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
virtual void v_NumFluxforScalar(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &ufield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&uflux);
virtual void v_NumericalFluxO1(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
&numericalFluxO1);
virtual void v_WeakPenaltyforScalar(
virtual void v_WeakPenaltyO1(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const int var,
const Array<OneD, const NekDouble> &ufield,
Array<OneD, NekDouble> &penaltyflux,
NekDouble time);
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &penaltyfluxO1);
virtual void v_NumFluxforVector(
virtual void v_NumericalFluxO2(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &ufield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&qfield,
Array<OneD, Array<OneD, NekDouble> > &qflux);
virtual void v_WeakPenaltyforVector(
virtual void v_WeakPenaltyO2(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const int var,
const int dir,
const Array<OneD, const NekDouble> &qfield,
Array<OneD, NekDouble> &penaltyflux,
NekDouble C11,
NekDouble time);
Array<OneD, NekDouble> &penaltyflux);
};
}
}
#endif
This diff is collapsed.
///////////////////////////////////////////////////////////////////////////////
//
// File: DiffusionLFR.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: LFR diffusion class.
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_SOLVERUTILS_DIFFUSIONLFR
#define NEKTAR_SOLVERUTILS_DIFFUSIONLFR
#include <SolverUtils/Diffusion/Diffusion.h>
namespace Nektar
{
namespace SolverUtils
{
class DiffusionLFR : public Diffusion
{
public:
static DiffusionSharedPtr create(std::string diffType)
{
return DiffusionSharedPtr(new DiffusionLFR(diffType));
}
static std::string type[];
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e0;
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e1;
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e2;
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e3;
Array<OneD, Array<OneD, NekDouble> > m_dGL_xi1;
Array<OneD, Array<OneD, NekDouble> > m_dGR_xi1;
Array<OneD, Array<OneD, NekDouble> > m_dGL_xi2;
Array<OneD, Array<OneD, NekDouble> > m_dGR_xi2;
Array<OneD, Array<OneD, NekDouble> > m_dGL_xi3;
Array<OneD, Array<OneD, NekDouble> > m_dGR_xi3;
DNekMatSharedPtr m_Ixm;
DNekMatSharedPtr m_Ixp;
protected:
DiffusionLFR(std::string diffType);
Array<OneD, Array<OneD, NekDouble> > m_traceNormals;
Array<OneD, Array<OneD, Array<OneD,NekDouble> > > m_tanbasis;
LibUtilities::SessionReaderSharedPtr m_session;
std::string m_diffType;
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_SetupMetrics(
LibUtilities::SessionReaderSharedPtr pSession,
</