Commit 9ea8a971 authored by Kilian Lackhove's avatar Kilian Lackhove

Merge branch 'feature/tidySolvers' into 'master'

Tidy EquationSystem and solvers

See merge request !832
parents 6ec5b72b 3b8e103c
......@@ -17,6 +17,8 @@ v5.0.0
- Added native support for csv files in addititon to pts (!760 !835)
- Utilize LAPACK_DIR env variable to find the native blas/lapack install (!827)
- 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)
- Simplify RawType typedefs (!840)
- Remove unused files from BasicUtils (!841)
- Remove checks for old boost versions which are no longer supported (!841)
......@@ -49,6 +51,10 @@ v5.0.0
- Enable output to multiple files (!844)
- Allow using xml file without expansion tag in FieldConvert (!849)
**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)
**CompressibleFlowSolver**
- Add 3D regression tests (!567)
......
......@@ -63,6 +63,92 @@ AdvectionSystem::~AdvectionSystem()
void AdvectionSystem::v_InitObject()
{
UnsteadySystem::v_InitObject();
m_session->LoadParameter("IO_CFLSteps", m_cflsteps, 0);
}
/**
*
*/
bool AdvectionSystem::v_PostIntegrate(int step)
{
bool result = UnsteadySystem::v_PostIntegrate(step);
if(m_cflsteps && !((step+1)%m_cflsteps))
{
int elmtid;
NekDouble cfl = GetCFLEstimate(elmtid);
if(m_comm->GetRank() == 0)
{
if( m_HomogeneousType == eNotHomogeneous)
{
cout << "CFL: ";
}
else
{
cout << "CFL (zero plane): ";
}
cout << cfl << " (in elmt " << elmtid << ")" << endl;
}
}
return result;
}
/**
*
*/
Array<OneD, NekDouble> AdvectionSystem::GetElmtCFLVals(void)
{
int nelmt = m_fields[0]->GetExpSize();
const Array<OneD, int> expOrder = GetNumExpModesPerExp();
const NekDouble cLambda = 0.2; // Spencer book pag. 317
Array<OneD, NekDouble> stdVelocity(nelmt, 0.0);
stdVelocity = v_GetMaxStdVelocity();
Array<OneD, NekDouble> cfl(nelmt, 0.0);
for(int el = 0; el < nelmt; ++el)
{
cfl[el] = m_timestep*(stdVelocity[el] * cLambda *
(expOrder[el]-1) * (expOrder[el]-1));
}
return cfl;
}
/**
*
*/
NekDouble AdvectionSystem::GetCFLEstimate(int &elmtid)
{
int n_element = m_fields[0]->GetExpSize();
Array<OneD, NekDouble> cfl = GetElmtCFLVals();
elmtid = Vmath::Imax(n_element,cfl,1);
NekDouble CFL,CFL_loc;
CFL = CFL_loc = cfl[elmtid];
m_comm->AllReduce(CFL,LibUtilities::ReduceMax);
// unshuffle elmt id if data is not stored in consecutive order.
elmtid = m_fields[0]->GetExp(elmtid)->GetGeom()->GetGlobalID();
if(CFL != CFL_loc)
{
elmtid = -1;
}
m_comm->AllReduce(elmtid,LibUtilities::ReduceMax);
// express element id with respect to plane
if(m_HomogeneousType == eHomogeneous1D)
{
elmtid = elmtid%m_fields[0]->GetPlane(0)->GetExpSize();
}
return CFL;
}
}
......
......@@ -54,14 +54,32 @@ public:
SOLVER_UTILS_EXPORT virtual void v_InitObject();
/// Returns the advection object held by this instance.
AdvectionSharedPtr GetAdvObject()
SOLVER_UTILS_EXPORT AdvectionSharedPtr GetAdvObject()
{
return m_advObject;
}
SOLVER_UTILS_EXPORT Array<OneD, NekDouble> GetElmtCFLVals(void);
SOLVER_UTILS_EXPORT NekDouble GetCFLEstimate(int &elmtid);
protected:
/// Advection term
SolverUtils::AdvectionSharedPtr m_advObject;
SOLVER_UTILS_EXPORT virtual bool v_PostIntegrate(int step);
SOLVER_UTILS_EXPORT virtual Array<OneD, NekDouble> v_GetMaxStdVelocity()
{
ASSERTL0(false,
"v_GetMaxStdVelocity is not implemented by the base class.");
Array<OneD, NekDouble> dummy(1);
return dummy;
}
private:
/// dump cfl estimate
int m_cflsteps;
};
/// Shared pointer to an AdvectionSystem class
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -54,10 +54,10 @@ namespace Nektar
/// Calculate the larger time-step mantaining the problem stable.
SOLVER_UTILS_EXPORT NekDouble GetTimeStep(
const Array<OneD, const Array<OneD, NekDouble> > &inarray);
/// CFL safety factor (comprise between 0 to 1).
NekDouble m_cflSafetyFactor;
protected:
/// Number of time steps between outputting status information.
int m_infosteps;
......@@ -81,6 +81,15 @@ namespace Nektar
/// forward transformed state.
bool m_homoInitialFwd;
/// Tolerance to which steady state should be evaluated at
NekDouble m_steadyStateTol;
/// Check for steady state at step interval
int m_steadyStateSteps;
/// Storage for previous solution for steady-state check
Array<OneD, Array<OneD, NekDouble> > m_previousSolution;
// Steady-state residual file
std::ofstream m_errFile;
std::vector<int> m_intVariables;
std::vector<FilterSharedPtr> m_filters;
......@@ -108,34 +117,11 @@ namespace Nektar
SOLVER_UTILS_EXPORT virtual void v_AppendOutput1D(
Array<OneD, Array<OneD, NekDouble> > &solution1D);
///
SOLVER_UTILS_EXPORT virtual void v_NumericalFlux(
Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &numflux);
///
SOLVER_UTILS_EXPORT virtual void v_NumericalFlux(
Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &numfluxX,
Array<OneD, Array<OneD, NekDouble> > &numfluxY );
///
SOLVER_UTILS_EXPORT virtual void v_NumFluxforScalar(
const Array<OneD, Array<OneD, NekDouble> > &ufield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &uflux);
///
SOLVER_UTILS_EXPORT virtual void v_NumFluxforVector(
const Array<OneD, Array<OneD, NekDouble> > &ufield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &qfield,
Array<OneD, Array<OneD, NekDouble> > &qflux);
SOLVER_UTILS_EXPORT virtual NekDouble v_GetTimeStep(
const Array<OneD, const Array<OneD, NekDouble> > &inarray);
SOLVER_UTILS_EXPORT virtual bool v_PreIntegrate(int step);
SOLVER_UTILS_EXPORT virtual bool v_PostIntegrate(int step);
SOLVER_UTILS_EXPORT virtual bool v_SteadyStateCheck(int step);
SOLVER_UTILS_EXPORT virtual bool v_RequireFwdTrans()
{
......@@ -153,21 +139,10 @@ namespace Nektar
private:
///
void WeakPenaltyforScalar(
const int var,
const Array<OneD, const NekDouble> &physfield,
Array<OneD, NekDouble> &penaltyflux,
NekDouble time=0.0);
///
void WeakPenaltyforVector(
const int var,
const int dir,
const Array<OneD, const NekDouble> &physfield,
Array<OneD, NekDouble> &penaltyflux,
NekDouble C11,
NekDouble time=0.0);
void InitializeSteadyState();
bool CheckSteadyState(int step);
};
}
......
......@@ -47,13 +47,14 @@ namespace Nektar
CFLtester::CFLtester(
const LibUtilities::SessionReaderSharedPtr& pSession)
: UnsteadySystem(pSession)
: UnsteadySystem(pSession),
AdvectionSystem(pSession)
{
}
void CFLtester::v_InitObject()
{
UnsteadySystem::v_InitObject();
AdvectionSystem::v_InitObject();
m_velocity = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
std::vector<std::string> vel;
......@@ -64,6 +65,60 @@ namespace Nektar
GetFunction( "AdvectionVelocity")->Evaluate(vel, m_velocity);
// Type of advection class to be used
switch(m_projectionType)
{
// Discontinuous field
case MultiRegions::eDiscontinuous:
{
// Do not forwards transform initial condition
m_homoInitialFwd = false;
// Define the normal velocity fields
if (m_fields[0]->GetTrace())
{
m_traceVn = Array<OneD, NekDouble>(GetTraceNpoints());
}
string advName;
string riemName;
m_session->LoadSolverInfo(
"AdvectionType", advName, "WeakDG");
m_advObject = SolverUtils::
GetAdvectionFactory().CreateInstance(advName, advName);
m_advObject->SetFluxVector(
&CFLtester::GetFluxVector, this);
m_session->LoadSolverInfo(
"UpwindType", riemName, "Upwind");
m_riemannSolver = SolverUtils::
GetRiemannSolverFactory().CreateInstance(riemName);
m_riemannSolver->SetScalar(
"Vn", &CFLtester::GetNormalVelocity, this);
m_advObject->SetRiemannSolver(m_riemannSolver);
m_advObject->InitObject(m_session, m_fields);
break;
}
// Continuous field
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
string advName;
m_session->LoadSolverInfo(
"AdvectionType", advName, "NonConservative");
m_advObject = SolverUtils::
GetAdvectionFactory().CreateInstance(advName, advName);
m_advObject->SetFluxVector(
&CFLtester::GetFluxVector, this);
break;
}
default:
{
ASSERTL0(false, "Unsupported projection type.");
break;
}
}
if (m_explicitAdvection)
{
m_ode.DefineOdeRhs (&CFLtester::DoOdeRhs, this);
......@@ -84,47 +139,14 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> >&outarray,
const NekDouble time)
{
int j;
int nvariables = inarray.num_elements();
int npoints = GetNpoints();
switch (m_projectionType)
m_advObject->Advect(nvariables, m_fields, m_velocity, inarray,
outarray, time);
for(int j = 0; j < nvariables; ++j)
{
case MultiRegions::eDiscontinuous:
{
int ncoeffs = inarray[0].num_elements();
Array<OneD, Array<OneD, NekDouble> > WeakAdv(nvariables);
WeakAdv[0] = Array<OneD, NekDouble>(ncoeffs*nvariables);
for(j = 1; j < nvariables; ++j)
{
WeakAdv[j] = WeakAdv[j-1] + ncoeffs;
}
WeakDGAdvection(inarray, WeakAdv, true, true);
for(j = 0; j < nvariables; ++j)
{
m_fields[j]->MultiplyByElmtInvMass(WeakAdv[j], WeakAdv[j]);
m_fields[j]->BwdTrans(WeakAdv[j],outarray[j]);
Vmath::Neg(npoints,outarray[j],1);
}
break;
}
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
// Calculate -V\cdot Grad(u);
for(j = 0; j < nvariables; ++j)
{
AdvectionNonConservativeForm(m_velocity,
inarray[j],
outarray[j]);
Vmath::Neg(npoints, outarray[j], 1);
}
break;
}
Vmath::Neg(npoints,outarray[j],1);
}
}
......@@ -133,7 +155,6 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> >&outarray,
const NekDouble time)
{
int j;
int nvariables = inarray.num_elements();
SetBoundaryConditions(time);
......@@ -144,7 +165,7 @@ namespace Nektar
// Just copy over array
int npoints = GetNpoints();
for(j = 0; j < nvariables; ++j)
for(int j = 0; j < nvariables; ++j)
{
Vmath::Vcopy(npoints,inarray[j],1,outarray[j],1);
}
......@@ -155,7 +176,7 @@ namespace Nektar
{
Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs());
for(j = 0; j < nvariables; ++j)
for(int j = 0; j < nvariables; ++j)
{
m_fields[j]->FwdTrans(inarray[j],coeffs);
m_fields[j]->BwdTrans_IterPerExp(coeffs,outarray[j]);
......@@ -168,60 +189,56 @@ namespace Nektar
}
}
void CFLtester::v_GetFluxVector(
const int i,
Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &flux)
/**
* @brief Get the normal velocity
*/
Array<OneD, NekDouble> &CFLtester::GetNormalVelocity()
{
ASSERTL1(flux.num_elements() == m_velocity.num_elements(),
"Dimension of flux array and velocity array do not match");
// Number of trace (interface) points
int nTracePts = GetTraceNpoints();
for (int j = 0; j < flux.num_elements(); ++j)
// Auxiliary variable to compute the normal velocity
Array<OneD, NekDouble> tmp(nTracePts);
// Reset the normal velocity
Vmath::Zero(nTracePts, m_traceVn, 1);
for (int i = 0; i < m_velocity.num_elements(); ++i)
{
Vmath::Vmul(GetNpoints(),
physfield[i], 1,
m_velocity[j], 1,
flux[j], 1);
m_fields[0]->ExtractTracePhys(m_velocity[i], tmp);
Vmath::Vvtvp(nTracePts,
m_traceNormals[i], 1,
tmp, 1,
m_traceVn, 1,
m_traceVn, 1);
}
return m_traceVn;
}
void CFLtester::v_NumericalFlux(
Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &numflux)
void CFLtester::GetFluxVector(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux)
{
int i;
int nTraceNumPoints = GetTraceNpoints();
int nvel = m_spacedim;
Array<OneD, NekDouble > Fwd(nTraceNumPoints);
Array<OneD, NekDouble > Bwd(nTraceNumPoints);
Array<OneD, NekDouble > Vn (nTraceNumPoints,0.0);
ASSERTL1(flux[0].num_elements() == m_velocity.num_elements(),
"Dimension of flux array and velocity array do not match");
// Get Edge Velocity - Could be stored if time independent
for (i = 0; i < nvel; ++i)
{
m_fields[0]->ExtractTracePhys(m_velocity[i], Fwd);
Vmath::Vvtvp(nTraceNumPoints,
m_traceNormals[i], 1,
Fwd, 1, Vn, 1, Vn, 1);
}
int nq = physfield[0].num_elements();
for (i = 0; i < numflux.num_elements(); ++i)
for (int i = 0; i < flux.num_elements(); ++i)
{
m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
//evaulate upwinded m_fields[i]
m_fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux[i]);
// calculate m_fields[i]*Vn
Vmath::Vmul(nTraceNumPoints, numflux[i], 1, Vn, 1, numflux[i], 1);
for (int j = 0; j < flux[0].num_elements(); ++j)
{
Vmath::Vmul(nq, physfield[i], 1, m_velocity[j], 1,
flux[i][j], 1);
}
}
}
void CFLtester::v_GenerateSummary(SummaryList& s)
{
UnsteadySystem::v_GenerateSummary(s);
AdvectionSystem::v_GenerateSummary(s);
}
......
......@@ -36,7 +36,8 @@
#ifndef NEKTAR_SOLVERS_ADRSOLVER_EQUATIONSYSTEMS_CFLTESTER_H
#define NEKTAR_SOLVERS_ADRSOLVER_EQUATIONSYSTEMS_CFLTESTER_H
#include <SolverUtils/UnsteadySystem.h>
#include <SolverUtils/AdvectionSystem.h>
#include <SolverUtils/RiemannSolvers/RiemannSolver.h>
using namespace Nektar::SolverUtils;
......@@ -78,7 +79,7 @@ namespace Nektar
{724.129308, 1297.665239, 1890.995890, 2477.394447, 3025.267650, 3552.997300, 4093.052410, 4609.237151, 5136.747279, 5661.792457, 6397.687000, 7383.760000, 8431.678000, 9540.598000},
{732.175780, 1305.179484, 1897.615616, 2483.199183, 3030.256925, 3557.347822, 4096.823380, 4612.552660, 5139.604517, 5664.300842, 6771.812000, 7388.760000, 8436.678000, 9546.598000}};
class CFLtester : public UnsteadySystem
class CFLtester : public SolverUtils::AdvectionSystem
{
public:
friend class MemoryManager<CFLtester>;
......@@ -93,8 +94,9 @@ namespace Nektar
virtual ~CFLtester();
protected:
Array<OneD, Array<OneD, NekDouble> > m_velocity;
SolverUtils::RiemannSolverSharedPtr m_riemannSolver;
Array<OneD, Array<OneD, NekDouble> > m_velocity;
Array<OneD, NekDouble> m_traceVn;
CFLtester(const LibUtilities::SessionReaderSharedPtr& pSession);
......@@ -106,11 +108,14 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble time);
virtual void v_InitObject();
/// Get the normal velocity
Array<OneD, NekDouble> &GetNormalVelocity();
virtual void v_GetFluxVector(const int i, Array<OneD, Array<OneD, NekDouble> > &physfield, Array<OneD, Array<OneD, NekDouble> > &flux);
virtual void v_InitObject();
virtual void v_NumericalFlux(Array<OneD, Array<OneD, NekDouble> > &physfield, Array<OneD, Array<OneD, NekDouble> > &numflux);
void GetFluxVector(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux);
virtual void v_GenerateSummary(SummaryList& s);
private:
......
......@@ -62,6 +62,85 @@ namespace Nektar
vel.resize(m_spacedim);
GetFunction( "AdvectionVelocity")->Evaluate(vel, m_velocity);
// Type of advection class to be used
switch(m_projectionType)
{
// Discontinuous field
case MultiRegions::eDiscontinuous:
{
// Define the normal velocity fields
if (m_fields[0]->GetTrace())
{
m_traceVn = Array<OneD, NekDouble>(GetTraceNpoints());
}
string advName;
string riemName;
m_session->LoadSolverInfo(
"AdvectionType", advName, "WeakDG");
m_advObject = SolverUtils::
GetAdvectionFactory().CreateInstance(advName, advName);
m_advObject->SetFluxVector(
&EigenValuesAdvection::GetFluxVector, this);
m_session->LoadSolverInfo(
"UpwindType", riemName, "Upwind");
m_riemannSolver = SolverUtils::
GetRiemannSolverFactory().CreateInstance(riemName);
m_riemannSolver->SetScalar(
"Vn", &EigenValuesAdvection::GetNormalVelocity, this);
m_advObject->SetRiemannSolver(m_riemannSolver);
m_advObject->InitObject(m_session, m_fields);
break;
}
// Continuous field
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
string advName;
m_session->LoadSolverInfo(
"AdvectionType", advName, "NonConservative");
m_advObject = SolverUtils::
GetAdvectionFactory().CreateInstance(advName, advName);
m_advObject->SetFluxVector(
&EigenValuesAdvection::GetFluxVector, this);
break;
}
default:
{
ASSERTL0(false, "Unsupported projection type.");
break;
}
}
}
/**
* @brief Get the normal velocity
*/
Array<OneD, NekDouble> &EigenValuesAdvection::GetNormalVelocity()
{
// Number of trace (interface) points
int nTracePts = GetTraceNpoints();
// Auxiliary variable to compute the normal velocity
Array<OneD, NekDouble> tmp(nTracePts);
// Reset the normal velocity
Vmath::Zero(nTracePts, m_traceVn, 1);
for (int i = 0; i < m_velocity.num_elements(); ++i)
{
m_fields[0]->ExtractTracePhys(m_velocity[i], tmp);
Vmath::Vvtvp(nTracePts,
m_traceNormals[i], 1,
tmp, 1,
m_traceVn, 1,
m_traceVn, 1);
}
return m_traceVn;
}
void EigenValuesAdvection::v_DoInitialise()
......@@ -77,7 +156,7 @@ namespace Nektar
void EigenValuesAdvection::v_DoSolve()
{
int nvariables = 1;
int i,dofs = GetNcoeffs();
int dofs = GetNcoeffs();
//bool UseContCoeffs = false;
Array<OneD, Array<OneD, NekDouble> > inarray(nvariables);
......@@ -126,36 +205,20 @@ namespace Nektar
/// we simulate the multiplication by the identity matrix.
/// The results stored in outarray is one of the columns of the weak advection oprators
/// which are then stored in MATRIX for the futher eigenvalues calculation.
m_advObject->Advect(nvariables, m_fields, m_velocity, inarray,
outarray, 0.0);
Vmath::Neg(npoints,outarray[0],1);
switch (m_projectionType)
{
case MultiRegions::eDiscontinuous:
{
WeakDGAdvection(inarray, WeakAdv,true,true,1);
m_fields[0]->MultiplyByElmtInvMass(WeakAdv[0],WeakAdv[0]);
m_fields[0]->BwdTrans(WeakAdv[0],outarray[0]);
Vmath::Neg(npoints,outarray[0],1);
break;
}
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
// Calculate -V\cdot Grad(u);
for(i = 0; i < nvariables; ++i)
for(int i = 0; i < nvariables; ++i)
{
//Projection
m_fields[i]->FwdTrans(inarray[i],WeakAdv[i]);
m_fields[i]->BwdTrans_IterPerExp(WeakAdv[i],tmp[i]);
//Advection operator
AdvectionNonConservativeForm(m_velocity,tmp[i],outarray[i]);