Commit 2477c58d authored by Dave Moxey's avatar Dave Moxey

Merge branch 'feature/CFS-filters' into 'master'

Extend AeroForces filter to compressible flows

See merge request !815
parents 982d929a 149c6a86
......@@ -16,6 +16,7 @@ v5.0.0
- Add ARPACK thirdparty build capabilities (!828)
- Added native support for csv files in addititon to pts (!760, !835, !906)
- Utilize LAPACK_DIR env variable to find the native blas/lapack install (!827)
- Extend AeroForces filter to compressible flows (!815)
- Remove StdExpansion use from MultiRegion (use Expansions instead). (!831)
- Move steady state check and CFL output from solvers to SolverUtils (!832)
- Remove DG advection implementation from EquationSystem (!832)
......
......@@ -24,7 +24,7 @@ SET(SOLVER_UTILS_SOURCES
Filters/FilterAverageFields.cpp
Filters/FilterCheckpoint.cpp
Filters/FilterEnergy1D.cpp
Filters/FilterEnergyBase.cpp
Filters/FilterEnergy.cpp
Filters/FilterHistoryPoints.cpp
Filters/FilterModalEnergy.cpp
Filters/FilterMovingAverage.cpp
......@@ -68,7 +68,7 @@ SET(SOLVER_UTILS_HEADERS
Filters/FilterAverageFields.h
Filters/FilterCheckpoint.h
Filters/FilterEnergy1D.h
Filters/FilterEnergyBase.h
Filters/FilterEnergy.h
Filters/FilterHistoryPoints.h
Filters/FilterModalEnergy.h
Filters/FilterMovingAverage.h
......
......@@ -1139,7 +1139,7 @@ namespace Nektar
variables[i] = m_boundaryConditions->GetVariable(i);
}
v_ExtraFldOutput(fieldcoeffs, variables);
ExtraFldOutput(fieldcoeffs, variables);
WriteFld(outname, m_fields[0], fieldcoeffs, variables);
}
......
......@@ -128,7 +128,11 @@ class Interpolator;
}
/// Get pressure field if available
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure();
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure();
SOLVER_UTILS_EXPORT inline void ExtraFldOutput(
std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
std::vector<std::string> &variables);
/// Print a summary of parameters and solver characteristics.
SOLVER_UTILS_EXPORT inline void PrintSummary(std::ostream &out);
......@@ -567,7 +571,21 @@ class Interpolator;
{
return v_GetPressure();
}
/**
* Append the coefficients and name of variables with solver specific
* extra variables
*
* @param fieldcoeffs Vector with coefficients
* @param variables Vector with name of variables
*/
inline void EquationSystem::ExtraFldOutput(
std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
std::vector<std::string> &variables)
{
v_ExtraFldOutput(fieldcoeffs, variables);
}
/**
* Prints a summary of variables and problem parameters.
*
......
......@@ -46,8 +46,10 @@ FilterFactory& GetFilterFactory()
}
Filter::Filter(const LibUtilities::SessionReaderSharedPtr& pSession) :
m_session(pSession)
Filter::Filter(const LibUtilities::SessionReaderSharedPtr &pSession,
const std::weak_ptr<EquationSystem> &pEquation) :
m_session(pSession),
m_equ(pEquation)
{
}
......
......@@ -48,6 +48,7 @@ namespace Nektar
namespace SolverUtils
{
class Filter;
class EquationSystem;
/// A shared pointer to a Driver object
typedef std::shared_ptr<Filter> FilterSharedPtr;
......@@ -57,6 +58,7 @@ typedef std::shared_ptr<Filter> FilterSharedPtr;
typedef LibUtilities::NekFactory<
std::string, Filter,
const LibUtilities::SessionReaderSharedPtr&,
const std::weak_ptr<EquationSystem>&,
const std::map<std::string, std::string>&
> FilterFactory;
SOLVER_UTILS_EXPORT FilterFactory& GetFilterFactory();
......@@ -66,7 +68,8 @@ class Filter
public:
typedef std::map<std::string, std::string> ParamMap;
SOLVER_UTILS_EXPORT Filter(
const LibUtilities::SessionReaderSharedPtr& pSession);
const LibUtilities::SessionReaderSharedPtr &pSession,
const std::weak_ptr<EquationSystem> &pEquation);
SOLVER_UTILS_EXPORT virtual ~Filter();
SOLVER_UTILS_EXPORT inline void Initialise(
......@@ -81,7 +84,8 @@ public:
SOLVER_UTILS_EXPORT inline bool IsTimeDependent();
protected:
LibUtilities::SessionReaderSharedPtr m_session;
LibUtilities::SessionReaderSharedPtr m_session;
const std::weak_ptr<EquationSystem> m_equ;
virtual void v_Initialise(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
......
......@@ -43,6 +43,7 @@
#include <MultiRegions/ExpList3D.h>
#include <MultiRegions/ExpList3DHomogeneous1D.h>
#include <SolverUtils/Filters/FilterAeroForces.h>
#include <SolverUtils/Filters/FilterInterfaces.hpp>
#include <LibUtilities/BasicUtils/ParseUtils.h>
using namespace std;
......@@ -60,8 +61,9 @@ std::string FilterAeroForces::className =
*/
FilterAeroForces::FilterAeroForces(
const LibUtilities::SessionReaderSharedPtr &pSession,
const std::weak_ptr<EquationSystem> &pEquation,
const ParamMap &pParams) :
Filter(pSession)
Filter(pSession, pEquation)
{
// OutputFile
auto it = pParams.find("OutputFile");
......@@ -584,6 +586,13 @@ void FilterAeroForces::CalculateForces(
return;
}
// Lock equation system weak pointer
auto equ = m_equ.lock();
ASSERTL0(equ, "Weak pointer expired");
auto fluidEqu = std::dynamic_pointer_cast<FluidInterface>(equ);
ASSERTL0(fluidEqu, "Aero forces filter is incompatible with this solver.");
int i, j, k, n, cnt, elmtid, nq, offset, boundary, plane;
// Get number of quadrature points and dimensions
int physTot = pFields[0]->GetNpoints();
......@@ -600,12 +609,18 @@ void FilterAeroForces::CalculateForces(
Array<OneD, MultiRegions::ExpListSharedPtr>
fields( pFields.num_elements() );
// Arrays of variables in field
Array<OneD, Array<OneD, NekDouble> > physfields(pFields.num_elements());
Array<OneD, Array<OneD, NekDouble> > velocity(nVel);
Array<OneD, NekDouble> pressure;
// Arrays of variables in the element
Array<OneD, Array<OneD, NekDouble> > velocity(expdim);
Array<OneD, NekDouble> P(physTot);
Array<OneD, Array<OneD, NekDouble> > velElmt(expdim);
Array<OneD, NekDouble> pElmt(physTot);
// Velocity gradient
Array<OneD, Array<OneD, NekDouble> > grad( expdim*expdim);
Array<OneD, NekDouble> div;
// Values at the boundary
Array<OneD, NekDouble> Pb;
......@@ -635,10 +650,18 @@ void FilterAeroForces::CalculateForces(
Array<OneD, Array<OneD, NekDouble> > fv( expdim );
// Get viscosity
NekDouble rho = (m_session->DefinesParameter("rho"))
? (m_session->GetParameter("rho"))
: 1;
NekDouble mu = rho*m_session->GetParameter("Kinvis");
NekDouble mu;
if(m_session->DefinesParameter("Kinvis"))
{
NekDouble rho = (m_session->DefinesParameter("rho"))
? (m_session->GetParameter("rho")) : 1;
mu = rho * m_session->GetParameter("Kinvis");
}
else
{
mu = m_session->GetParameter("mu");
}
NekDouble lambda = -2.0/3.0;
// Perform BwdTrans: when we only want the mean force in a 3DH1D
// we work in wavespace, otherwise we use physical space
......@@ -659,12 +682,12 @@ void FilterAeroForces::CalculateForces(
if(m_isHomogeneous1D)
{
BndConds = pFields[0]->GetPlane(0)->GetBndConditions();
BndExp = pFields[0]->GetPlane(0)->GetBndCondExpansions();
BndExp = pFields[0]->GetPlane(0)->GetBndCondExpansions();
}
else
{
BndConds = pFields[0]->GetBndConditions();
BndExp = pFields[0]->GetBndCondExpansions();
BndExp = pFields[0]->GetBndCondExpansions();
}
// For Homogeneous, calculate force on each 2D plane
......@@ -689,7 +712,20 @@ void FilterAeroForces::CalculateForces(
{
fields[n] = pFields[n];
}
}
}
// Get velocity and pressure values
for(n = 0; n < physfields.num_elements(); ++n)
{
physfields[n] = fields[n]->GetPhys();
}
for(n = 0; n < nVel; ++n)
{
velocity[n] = Array<OneD, NekDouble>(fields[n]->GetTotPoints());
}
pressure = Array<OneD, NekDouble>(fields[0]->GetTotPoints());
fluidEqu->GetVelocity(physfields, velocity);
fluidEqu->GetPressure(physfields, pressure);
//Loop all the Boundary Regions
for( cnt = n = 0; n < BndConds.num_elements(); n++)
......@@ -706,21 +742,30 @@ void FilterAeroForces::CalculateForces(
// Extract fields on this element
for( j=0; j<expdim; j++)
{
velocity[j] = fields[j]->GetPhys() + offset;
velElmt[j] = velocity[j] + offset;
}
P = fields[nVel]->GetPhys() + offset;
pElmt = pressure + offset;
// Compute the velocity gradients
div = Array<OneD, NekDouble>(nq,0.0);
for (j=0; j<expdim; j++)
{
for (k=0; k<expdim; k++)
{
grad[j*expdim+k] =
Array<OneD, NekDouble>(nq,0.0);
elmt->PhysDeriv(k,velocity[j],
elmt->PhysDeriv(k,velElmt[j],
grad[j*expdim+k]);
if( j == k)
{
Vmath::Vadd(nq, grad[j*expdim+k], 1,
div, 1, div, 1);
}
}
}
// Scale div by lambda (for compressible flows)
Vmath::Smul(nq, lambda, div, 1, div, 1);
// identify boundary of element
boundary = m_BCtoTraceID[cnt];
......@@ -746,7 +791,7 @@ void FilterAeroForces::CalculateForces(
// Extract values at boundary
Pb = Array<OneD, NekDouble>(nbc,0.0);
elmt->GetEdgePhysVals(boundary,bc,P,Pb);
elmt->GetEdgePhysVals(boundary,bc,pElmt,Pb);
for(int j = 0; j < expdim*expdim; ++j)
{
gradb[j] = Array<OneD, NekDouble>
......@@ -772,14 +817,14 @@ void FilterAeroForces::CalculateForces(
// Extract values at boundary
Pb = Array<OneD, NekDouble>(nbc,0.0);
elmt->GetFacePhysVals(boundary,bc,P,Pb);
elmt->GetFacePhysVals(boundary,bc,pElmt,Pb);
for(int j = 0; j < expdim*expdim; ++j)
{
gradb[j] = Array<OneD, NekDouble>
(nbc,0.0);
elmt->GetFacePhysVals(boundary,
bc,grad[j],gradb[j]);
}
}
}
break;
......@@ -814,7 +859,15 @@ void FilterAeroForces::CalculateForces(
Vmath::Vvtvp (nbc, gradb[j*expdim+k], 1,
normals[k], 1,
fv[j], 1,
fv[j], 1);
fv[j], 1);
}
if(!fluidEqu->HasConstantDensity())
{
// Add gradient term
Vmath::Vvtvp (nbc, div, 1,
normals[j], 1,
fv[j], 1,
fv[j], 1);
}
Vmath::Smul(nbc, -mu, fv[j], 1, fv[j], 1);
}
......@@ -888,7 +941,7 @@ void FilterAeroForces::CalculateForces(
pFields[i]->HomogeneousFwdTrans(pFields[i]->GetPhys(),
pFields[i]->UpdatePhys());
}
}
}
}
/**
......@@ -927,13 +980,13 @@ void FilterAeroForces::CalculateForcesMapping(
Array<OneD, MultiRegions::ExpListSharedPtr> C ( nVel*nVel);
Array<OneD, MultiRegions::ExpListSharedPtr> CPlane( nVel*nVel);
Array<OneD, Array<OneD, NekDouble> > CElmt ( nVel*nVel);
Array<OneD, Array<OneD, NekDouble> > CBnd ( nVel*nVel);
Array<OneD, Array<OneD, NekDouble> > CBnd ( nVel*nVel);
// Jacobian
MultiRegions::ExpListSharedPtr Jac;
MultiRegions::ExpListSharedPtr JacPlane;
Array<OneD, NekDouble> JacElmt;
Array<OneD, NekDouble> JacBnd;
Array<OneD, NekDouble> JacBnd;
// Communicators to exchange results
LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
......@@ -983,12 +1036,12 @@ void FilterAeroForces::CalculateForcesMapping(
if(m_isHomogeneous1D)
{
BndConds = pFields[0]->GetPlane(0)->GetBndConditions();
BndExp = pFields[0]->GetPlane(0)->GetBndCondExpansions();
BndExp = pFields[0]->GetPlane(0)->GetBndCondExpansions();
}
else
{
BndConds = pFields[0]->GetBndConditions();
BndExp = pFields[0]->GetBndCondExpansions();
BndExp = pFields[0]->GetBndCondExpansions();
}
//
......@@ -1018,7 +1071,7 @@ void FilterAeroForces::CalculateForcesMapping(
AllocateSharedPtr(*Exp3DH1);
}
Jac = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::
AllocateSharedPtr(*Exp3DH1);
AllocateSharedPtr(*Exp3DH1);
}
else
{
......@@ -1038,7 +1091,7 @@ void FilterAeroForces::CalculateForcesMapping(
AllocateSharedPtr(*Exp2D);
}
Jac = MemoryManager<MultiRegions::ExpList2D>::
AllocateSharedPtr(*Exp2D);
AllocateSharedPtr(*Exp2D);
}
break;
}
......
......@@ -52,11 +52,11 @@ public:
/// Creates an instance of this class
static FilterSharedPtr create(
const LibUtilities::SessionReaderSharedPtr &pSession,
const std::map<std::string, std::string> &pParams)
const std::weak_ptr<EquationSystem> &pEquation,
const std::map<std::string, std::string> &pParams)
{
FilterSharedPtr p = MemoryManager<FilterAeroForces>::
AllocateSharedPtr(pSession, pParams);
//p->InitObject();
AllocateSharedPtr(pSession, pEquation, pParams);
return p;
}
......@@ -65,14 +65,15 @@ public:
SOLVER_UTILS_EXPORT FilterAeroForces(
const LibUtilities::SessionReaderSharedPtr &pSession,
const std::map<std::string, std::string> &pParams);
const std::weak_ptr<EquationSystem> &pEquation,
const std::map<std::string, std::string> &pParams);
SOLVER_UTILS_EXPORT virtual ~FilterAeroForces();
SOLVER_UTILS_EXPORT void GetForces(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
Array<OneD, NekDouble> &Aeroforces,
const NekDouble &time);
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
Array<OneD, NekDouble> &Aeroforces,
const NekDouble &time);
protected:
virtual void v_Initialise(
......@@ -117,7 +118,7 @@ private:
Array<OneD, Array<OneD, NekDouble> > m_Ftplane;
NekDouble m_lastTime;
GlobalMapping::MappingSharedPtr m_mapping;
GlobalMapping::MappingSharedPtr m_mapping;
void CalculateForces(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
......
......@@ -45,9 +45,22 @@ std::string FilterAverageFields::className =
FilterAverageFields::FilterAverageFields(
const LibUtilities::SessionReaderSharedPtr &pSession,
const std::weak_ptr<EquationSystem> &pEquation,
const ParamMap &pParams)
: FilterFieldConvert(pSession, pParams)
: FilterFieldConvert(pSession, pEquation, pParams)
{
// Load sampling frequency
auto it = pParams.find("SampleFrequency");
if (it == pParams.end())
{
m_sampleFrequency = 1;
}
else
{
LibUtilities::Equation equ(
m_session->GetExpressionEvaluator(), it->second);
m_sampleFrequency = round(equ.Evaluate());
}
}
FilterAverageFields::~FilterAverageFields()
......@@ -56,12 +69,13 @@ FilterAverageFields::~FilterAverageFields()
void FilterAverageFields::v_ProcessSample(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
const NekDouble &time)
{
for(int n = 0; n < pFields.num_elements(); ++n)
for(int n = 0; n < m_outFields.size(); ++n)
{
Vmath::Vadd(m_outFields[n].num_elements(),
pFields[n]->GetCoeffs(),
fieldcoeffs[n],
1,
m_outFields[n],
1,
......
......@@ -50,10 +50,11 @@ public:
/// Creates an instance of this class
static FilterSharedPtr create(
const LibUtilities::SessionReaderSharedPtr &pSession,
const std::map<std::string, std::string> &pParams)
const std::weak_ptr<EquationSystem> &pEquation,
const std::map<std::string, std::string> &pParams)
{
FilterSharedPtr p = MemoryManager<FilterAverageFields>
::AllocateSharedPtr(pSession, pParams);
::AllocateSharedPtr(pSession, pEquation, pParams);
return p;
}
......@@ -62,12 +63,14 @@ public:
SOLVER_UTILS_EXPORT FilterAverageFields(
const LibUtilities::SessionReaderSharedPtr &pSession,
const std::weak_ptr<EquationSystem> &pEquation,
const ParamMap &pParams);
SOLVER_UTILS_EXPORT virtual ~FilterAverageFields();
protected:
virtual void v_ProcessSample(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
const NekDouble &time);
virtual void v_PrepareOutput(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
......
......@@ -45,8 +45,9 @@ std::string FilterCheckpoint::className =
FilterCheckpoint::FilterCheckpoint(
const LibUtilities::SessionReaderSharedPtr &pSession,
const std::weak_ptr<EquationSystem> &pEquation,
const ParamMap &pParams) :
Filter(pSession)
Filter(pSession, pEquation)
{
// OutputFile
auto it = pParams.find("OutputFile");
......
......@@ -50,10 +50,11 @@ public:
/// Creates an instance of this class
static FilterSharedPtr create(
const LibUtilities::SessionReaderSharedPtr &pSession,
const std::map<std::string, std::string> &pParams) {
const std::weak_ptr<EquationSystem> &pEquation,
const std::map<std::string, std::string> &pParams)
{
FilterSharedPtr p = MemoryManager<FilterCheckpoint>
::AllocateSharedPtr(pSession, pParams);
//p->InitObject();
::AllocateSharedPtr(pSession, pEquation, pParams);
return p;
}
......@@ -62,6 +63,7 @@ public:
SOLVER_UTILS_EXPORT FilterCheckpoint(
const LibUtilities::SessionReaderSharedPtr &pSession,
const std::weak_ptr<EquationSystem> &pEquation,
const ParamMap &pParams);
SOLVER_UTILS_EXPORT virtual ~FilterCheckpoint();
......
///////////////////////////////////////////////////////////////////////////////
//
// File FilterEnergyBase.cpp
// File FilterEnergy.cpp
//
// For more information, please see: http://www.nektar.info
//
......@@ -35,7 +35,8 @@
#include <iomanip>
#include <SolverUtils/Filters/FilterEnergyBase.h>
#include <SolverUtils/Filters/FilterEnergy.h>
#include <SolverUtils/Filters/FilterInterfaces.hpp>
using namespace std;
......@@ -43,15 +44,17 @@ namespace Nektar
{
namespace SolverUtils
{
FilterEnergyBase::FilterEnergyBase(
std::string FilterEnergy::className = SolverUtils::GetFilterFactory().
RegisterCreatorFunction("Energy", FilterEnergy::create);
FilterEnergy::FilterEnergy(
const LibUtilities::SessionReaderSharedPtr &pSession,
const ParamMap &pParams,
const bool pConstDensity)
: Filter (pSession),
const std::weak_ptr<EquationSystem> &pEquation,
const ParamMap &pParams)
: Filter (pSession, pEquation),
m_index (-1),
m_homogeneous (false),
m_planes (),
m_constDensity(pConstDensity)
m_planes ()
{
std::string outName;
......@@ -90,12 +93,12 @@ FilterEnergyBase::FilterEnergyBase(
m_outputFrequency = round(equ.Evaluate());
}
FilterEnergyBase::~FilterEnergyBase()
FilterEnergy::~FilterEnergy()
{
}
void FilterEnergyBase::v_Initialise(
void FilterEnergy::v_Initialise(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
const NekDouble &time)
{
......@@ -134,7 +137,7 @@ void FilterEnergyBase::v_Initialise(
v_Update(pFields, time);
}
void FilterEnergyBase::v_Update(
void FilterEnergy::v_Update(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
const NekDouble &time)
{
......@@ -147,16 +150,33 @@ void FilterEnergyBase::v_Update(
return;
}
// Lock equation system pointer
auto equ = m_equ.lock();
ASSERTL0(equ, "Weak pointer expired");
auto fluidEqu = std::dynamic_pointer_cast<FluidInterface>(equ);
ASSERTL0(fluidEqu, "Energy filter is incompatible with this solver.");
// Store physical values in an array
Array<OneD, Array<OneD, NekDouble> > physfields(pFields.num_elements());
for(i = 0; i < pFields.num_elements(); ++i)
{
physfields[i] = pFields[i]->GetPhys();
}
// Calculate kinetic energy.
NekDouble Ek = 0.0;
Array<OneD, NekDouble> tmp(nPoints, 0.0);