Commit e2001640 authored by Douglas Serson's avatar Douglas Serson

Merge branch 'feature/IncNS_DGKernel' into 'master'

Feature/inc ns dg kernel

Closes #60

See merge request !851
parents 4e36b131 3eaeea6a
......@@ -69,6 +69,7 @@ v5.0.0
**IncNavierStokesSolver**
- Replace steady-state check based on difference of norms by check based on
norm of the difference, to be consistent with the compressible solver (!832)
- Updated SVV to allow for the DGKernel extension (!851)
**CompressibleFlowSolver**
- Add 3D regression tests (!567)
......
......@@ -542,6 +542,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing)
{
......@@ -573,7 +574,7 @@ namespace Nektar
// Solve the system
GlobalLinSysKey key(StdRegions::eHelmholtz,
m_locToGloMap,factors,varcoeff);
m_locToGloMap,factors,varcoeff,varfactors);
if(flags.isSet(eUseGlobal))
{
......
......@@ -208,6 +208,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing);
......
......@@ -874,6 +874,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing)
......@@ -923,7 +924,7 @@ namespace Nektar
// Add weak boundary conditions to forcing
Vmath::Vadd(contNcoeffs, wsp, 1, gamma, 1, wsp, 1);
GlobalLinSysKey key(StdRegions::eHelmholtz,m_locToGloMap,factors,varcoeff);
GlobalLinSysKey key(StdRegions::eHelmholtz,m_locToGloMap,factors,varcoeff,varfactors);
if(flags.isSet(eUseGlobal))
{
......
......@@ -247,6 +247,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing);
......
......@@ -580,6 +580,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing)
{
......@@ -626,7 +627,7 @@ namespace Nektar
Vmath::Vadd(contNcoeffs, wsp, 1, gamma, 1, wsp, 1);
// Solve the system
GlobalLinSysKey key(StdRegions::eHelmholtz, m_locToGloMap, factors,varcoeff);
GlobalLinSysKey key(StdRegions::eHelmholtz, m_locToGloMap, factors,varcoeff,varfactors);
if(flags.isSet(eUseGlobal))
{
......
......@@ -184,6 +184,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing);
virtual void v_GeneralMatrixOp(
......
......@@ -254,6 +254,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing)
{
......@@ -304,7 +305,7 @@ namespace Nektar
m_planes[n]->HelmSolve(wfce,
e_out = outarray + cnt1,
flags, new_factors, varcoeff,
dirForcing,
varfactors, dirForcing,
PhysSpaceForcing);
}
......
......@@ -97,6 +97,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing);
......
......@@ -161,6 +161,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing)
{
......@@ -202,7 +203,8 @@ namespace Nektar
wfce = (PhysSpaceForcing)? fce+cnt:fce+cnt1;
m_lines[l]->HelmSolve(wfce,
e_out = outarray + cnt1,
flags, new_factors, varcoeff, dirForcing,
flags, new_factors, varcoeff, varfactors,
dirForcing,
PhysSpaceForcing);
cnt += m_lines[l]->GetTotPoints();
......
......@@ -92,6 +92,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing);
......
......@@ -1184,6 +1184,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing)
{
......
......@@ -263,6 +263,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing);
......
......@@ -1859,6 +1859,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing)
......
......@@ -232,6 +232,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing);
virtual void v_GeneralMatrixOp(
......
......@@ -2003,6 +2003,7 @@ using namespace boost::assign;
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing)
{
......
......@@ -191,6 +191,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing);
......
......@@ -266,6 +266,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing)
{
......@@ -303,7 +304,8 @@ namespace Nektar
m_planes[n]->HelmSolve(
wfce,
e_out = outarray + cnt1,
flags, new_factors, varcoeff, dirForcing,
flags, new_factors, varcoeff, varfactors,
dirForcing,
PhysSpaceForcing);
}
......
......@@ -278,6 +278,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing);
......
......@@ -215,6 +215,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing)
{
......@@ -255,7 +256,7 @@ namespace Nektar
m_lines[n]->HelmSolve(wfce,
e_out = outarray + cnt1,
flags, new_factors,
varcoeff, dirForcing,
varcoeff, varfactors,dirForcing,
PhysSpaceForcing);
cnt += m_lines[n]->GetTotPoints();
......
......@@ -140,6 +140,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing);
......
......@@ -2537,6 +2537,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing)
{
......@@ -2926,8 +2927,8 @@ namespace Nektar
Array<OneD, NekDouble> &bnd)
{
int n, cnt;
Array<OneD, NekDouble> tmp1, tmp2;
LocalRegions::ExpansionSharedPtr elmt;
Array<OneD, NekDouble> tmp1;
StdRegions::StdExpansionSharedPtr elmt;
Array<OneD, int> ElmtID,EdgeID;
GetBoundaryToElmtMap(ElmtID,EdgeID);
......@@ -2952,8 +2953,8 @@ namespace Nektar
elmt = GetExp(ElmtID[cnt+n]);
elmt->GetTracePhysVals(EdgeID[cnt+n],
GetBndCondExpansions()[i]->GetExp(n),
tmp1 = phys + offsetPhys,
tmp2 = bnd + offsetBnd);
phys + offsetPhys,
tmp1 = bnd + offsetBnd);
}
}
......
......@@ -46,6 +46,7 @@
#include <MultiRegions/MultiRegions.hpp>
#include <MultiRegions/GlobalMatrix.h>
#include <MultiRegions/GlobalMatrixKey.h>
#include <MultiRegions/GlobalLinSysKey.h>
#include <MultiRegions/GlobalOptimizationParameters.h>
#include <MultiRegions/AssemblyMap/AssemblyMap.h>
#include <tinyxml.h>
......@@ -283,6 +284,8 @@ namespace Nektar
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff =
StdRegions::NullVarCoeffMap,
const MultiRegions::VarFactorsMap &varfactors =
MultiRegions::NullVarFactorsMap,
const Array<OneD, const NekDouble> &dirForcing =
NullNekDouble1DArray,
const bool PhysSpaceForcing = true);
......@@ -1193,6 +1196,7 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing);
......@@ -1740,12 +1744,13 @@ namespace Nektar
const FlagList &flags,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varcoeff,
const MultiRegions::VarFactorsMap &varfactors,
const Array<OneD, const NekDouble> &dirForcing,
const bool PhysSpaceForcing)
{
v_HelmSolve(inarray, outarray, flags, factors, varcoeff,
dirForcing, PhysSpaceForcing);
varfactors, dirForcing, PhysSpaceForcing);
}
......
......@@ -213,7 +213,8 @@ namespace Nektar
*
* @param asmMap Assembly map used to construct the global system.
*/
PreconditionerSharedPtr GlobalLinSys::CreatePrecon(AssemblyMapSharedPtr asmMap)
PreconditionerSharedPtr GlobalLinSys::CreatePrecon(AssemblyMapSharedPtr
asmMap)
{
PreconditionerType pType = asmMap->GetPreconType();
std::string PreconType = MultiRegions::PreconditionerTypeMap[pType];
......@@ -233,43 +234,96 @@ namespace Nektar
return m_expList.lock()->GetExpSize();
}
/**
* @brief Retrieves the block matrix from n-th expansion using the
* matrix key provided by the #m_linSysKey.
*
* @param n Number of the expansion.
* @return Block matrix for the specified expansion.
*/
DNekScalMatSharedPtr GlobalLinSys::v_GetBlock(unsigned int n)
{
/**
Assemble the matrix key for each block n
**/
LocalRegions::MatrixKey GlobalLinSys::GetBlockMatrixKey(unsigned int n)
{
std::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
int cnt = 0;
DNekScalMatSharedPtr loc_mat;
LocalRegions::ExpansionSharedPtr vExp =
std::dynamic_pointer_cast<LocalRegions::Expansion>(
expList->GetExp(n));
LocalRegions::ExpansionSharedPtr vExp = expList->GetExp( n );
// need to be initialised with zero size for non variable
// coefficient case
StdRegions::VarCoeffMap vVarCoeffMap;
StdRegions::ConstFactorMap vConstFactorMap =
m_linSysKey.GetConstFactors();
// setup variable factors
if(m_linSysKey.GetNVarFactors() > 0)
{
if(m_linSysKey.GetVarFactors().
count(StdRegions::eFactorSVVDiffCoeff) != 0)
{
vConstFactorMap[StdRegions::eFactorSVVDiffCoeff] =
m_linSysKey.GetVarFactors(
StdRegions::eFactorSVVDiffCoeff)[n];
ASSERTL1(m_linSysKey.GetConstFactors().
count(StdRegions::eFactorSVVCutoffRatio),
"VarCoeffSVVCuroffRatio is set but "
" not FactorSVVCutoffRatio");
vConstFactorMap[StdRegions::eFactorSVVCutoffRatio] =
m_linSysKey.GetVarFactors(
StdRegions::eFactorSVVCutoffRatio)[n];
}
if(m_linSysKey.GetVarFactors().
count(StdRegions::eFactorSVVPowerKerDiffCoeff) != 0)
{
vConstFactorMap[StdRegions::eFactorSVVPowerKerDiffCoeff] =
m_linSysKey.GetVarFactors(
StdRegions::eFactorSVVPowerKerDiffCoeff)[n];
}
if(m_linSysKey.GetVarFactors().
count(StdRegions::eFactorSVVDGKerDiffCoeff) != 0)
{
vConstFactorMap[StdRegions::eFactorSVVDGKerDiffCoeff] =
m_linSysKey.GetVarFactors(
StdRegions::eFactorSVVDGKerDiffCoeff)[n];
}
}
// retrieve variable coefficients
if(m_linSysKey.GetNVarCoeffs() > 0)
{
cnt = expList->GetPhys_Offset(n);
for (auto &x : m_linSysKey.GetVarCoeffs())
{
vVarCoeffMap[x.first] = x.second + cnt;
}
}
LocalRegions::MatrixKey matkey(m_linSysKey.GetMatrixType(),
vExp->DetShapeType(),
*vExp, m_linSysKey.GetConstFactors(),
*vExp,
vConstFactorMap,
vVarCoeffMap);
loc_mat = vExp->GetLocMatrix(matkey);
return matkey;
}
/**
* @brief Retrieves the block matrix from n-th expansion using the
* matrix key provided by the #m_linSysKey.
*
* @param n Number of the expansion.
* @return Block matrix for the specified expansion.
*/
DNekScalMatSharedPtr GlobalLinSys::v_GetBlock(unsigned int n)
{
LocalRegions::ExpansionSharedPtr vExp = m_expList.lock()->GetExp( n );
DNekScalMatSharedPtr loc_mat;
loc_mat = vExp->GetLocMatrix(GetBlockMatrixKey(n));
// apply robin boundary conditions to the matrix.
if(m_robinBCInfo.count(n) != 0) // add robin mass matrix
......@@ -311,37 +365,14 @@ namespace Nektar
DNekScalBlkMatSharedPtr GlobalLinSys::v_GetStaticCondBlock(
unsigned int n)
{
std::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
int cnt = 0;
LocalRegions::ExpansionSharedPtr vExp = m_expList.lock()->GetExp( n );
DNekScalBlkMatSharedPtr loc_mat;
DNekScalMatSharedPtr tmp_mat;
LocalRegions::ExpansionSharedPtr vExp = expList->GetExp( n );
// need to be initialised with zero size for non variable
// coefficient case
StdRegions::VarCoeffMap vVarCoeffMap;
// retrieve variable coefficients
if(m_linSysKey.GetNVarCoeffs() > 0)
{
cnt = expList->GetPhys_Offset(n);
for (auto &x : m_linSysKey.GetVarCoeffs())
{
vVarCoeffMap[x.first] = x.second + cnt;
}
}
LocalRegions::MatrixKey matkey(m_linSysKey.GetMatrixType(),
vExp->DetShapeType(),
*vExp,
m_linSysKey.GetConstFactors(),
vVarCoeffMap);
loc_mat = vExp->GetLocStaticCondMatrix(matkey);
loc_mat = vExp->GetLocStaticCondMatrix(GetBlockMatrixKey(n));
if(m_robinBCInfo.count(n) != 0) // add robin mass matrix
{
DNekScalMatSharedPtr tmp_mat;
RobinBCInfoSharedPtr rBC;
tmp_mat = loc_mat->GetBlock(0,0);
......@@ -388,31 +419,8 @@ namespace Nektar
*/
void GlobalLinSys::v_DropStaticCondBlock(unsigned int n)
{
std::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
LocalRegions::ExpansionSharedPtr vExp = expList->GetExp( n );
// need to be initialised with zero size for non variable
// coefficient case
StdRegions::VarCoeffMap vVarCoeffMap;
// retrieve variable coefficients
if(m_linSysKey.GetNVarCoeffs() > 0)
{
int cnt = expList->GetPhys_Offset(n);
for (auto &x : m_linSysKey.GetVarCoeffs())
{
vVarCoeffMap[x.first] = x.second + cnt;
}
}
LocalRegions::MatrixKey matkey(m_linSysKey.GetMatrixType(),
vExp->DetShapeType(),
*vExp,
m_linSysKey.GetConstFactors(),
vVarCoeffMap);
vExp->DropLocStaticCondMatrix(matkey);
LocalRegions::ExpansionSharedPtr vExp = m_expList.lock()->GetExp( n );
vExp->DropLocStaticCondMatrix(GetBlockMatrixKey(n));
}
void GlobalLinSys::v_InitObject()
......
......@@ -139,6 +139,8 @@ namespace Nektar
PreconditionerSharedPtr CreatePrecon(AssemblyMapSharedPtr asmMap);
private:
LocalRegions::MatrixKey GetBlockMatrixKey(unsigned int n);
/// Solve a linear system based on mapping.
virtual void v_Solve(
const Array<OneD, const NekDouble> &in,
......
......@@ -59,11 +59,23 @@ namespace Nektar
GlobalLinSysKey::GlobalLinSysKey(const StdRegions::MatrixType matrixType,
const AssemblyMapSharedPtr &locToGloMap,
const StdRegions::ConstFactorMap &factors,
const StdRegions::VarCoeffMap &varCoeffs) :
const StdRegions::VarCoeffMap &varCoeffs,
const VarFactorsMap &varFactors) :
GlobalMatrixKey(matrixType, locToGloMap, factors, varCoeffs),
m_solnType(locToGloMap->GetGlobalSysSolnType())
m_solnType(locToGloMap->GetGlobalSysSolnType()),
m_varFactors(varFactors),
m_varFactors_hashes(varFactors.size())
{
// Create hash
int i = 0;
for (VarFactorsMap::const_iterator x = varFactors.begin();
x != varFactors.end(); ++x)
{
m_varFactors_hashes[i] = boost::hash_range(x->second.begin(),
x->second.begin() + x->second.num_elements());
boost::hash_combine(m_varFactors_hashes[i], (int)x->first);
i++;
}
}
......@@ -72,7 +84,9 @@ namespace Nektar
*/
GlobalLinSysKey::GlobalLinSysKey(const GlobalLinSysKey &key):
GlobalMatrixKey(key),
m_solnType(key.m_solnType)
m_solnType(key.m_solnType),
m_varFactors(key.m_varFactors),
m_varFactors_hashes(key.m_varFactors_hashes)
{
}
......@@ -105,6 +119,28 @@ namespace Nektar
return false;
}
if(lhs.m_varFactors.size() < rhs.m_varFactors.size())
{
return true;
}
if(lhs.m_varFactors.size() > rhs.m_varFactors.size())
{
return false;
}
for (unsigned int i = 0; i < lhs.m_varFactors_hashes.size(); ++i)
{
if(lhs.m_varFactors_hashes[i] < rhs.m_varFactors_hashes[i])
{
return true;
}
if(lhs.m_varFactors_hashes[i] > rhs.m_varFactors_hashes[i])
{
return false;
}
}
return (*dynamic_cast<const GlobalMatrixKey*>(&lhs)
< *dynamic_cast<const GlobalMatrixKey*>(&rhs));
}
......@@ -131,6 +167,9 @@ namespace Nektar
}
os << "Number of variable coefficients: "
<< rhs.GetNVarCoeffs() << endl;
os << "Number of variable factors : "
<< rhs.GetNVarFactors() << endl;
return os;
}
......
......@@ -43,6 +43,10 @@ namespace Nektar
{
namespace MultiRegions
{
typedef std::map<StdRegions::ConstFactorType, Array<OneD, NekDouble> > VarFactorsMap;
static VarFactorsMap NullVarFactorsMap;
/// Describe a linear system.
class GlobalLinSysKey : public GlobalMatrixKey