Commit 6b6938c6 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'feature/AutomatedSFD' into 'master'

Merge AutomatedSFD + new Advection class into master

See merge request !371
parents 1b4cb847 3b53fa85
......@@ -37,120 +37,168 @@
namespace Nektar
{
namespace SolverUtils
namespace SolverUtils
{
std::string Driver::evolutionOperatorLookupIds[6] = {
LibUtilities::SessionReader::RegisterEnumValue(
"EvolutionOperator","Nonlinear" ,eNonlinear),
LibUtilities::SessionReader::RegisterEnumValue(
"EvolutionOperator","Direct" ,eDirect),
LibUtilities::SessionReader::RegisterEnumValue(
"EvolutionOperator","Adjoint" ,eAdjoint),
LibUtilities::SessionReader::RegisterEnumValue(
"EvolutionOperator","TransientGrowth",eTransientGrowth),
LibUtilities::SessionReader::RegisterEnumValue(
"EvolutionOperator","SkewSymmetric" ,eSkewSymmetric),
LibUtilities::SessionReader::RegisterEnumValue(
"EvolutionOperator","AdaptiveSFD" ,eAdaptiveSFD)
};
std::string Driver::evolutionOperatorDef =
LibUtilities::SessionReader::RegisterDefaultSolverInfo(
"EvolutionOperator","Nonlinear");
std::string Driver::driverDefault =
LibUtilities::SessionReader::RegisterDefaultSolverInfo(
"Driver", "Standard");
DriverFactory& GetDriverFactory()
{
typedef Loki::SingletonHolder<DriverFactory,
Loki::CreateUsingNew,
Loki::NoDestroy > Type;
return Type::Instance();
}
/**
*
*/
Driver::Driver(const LibUtilities::SessionReaderSharedPtr pSession)
: m_comm(pSession->GetComm()),
m_session(pSession)
{
}
Driver::~Driver()
{
}
/**
*
*/
void Driver::v_InitObject(ostream &out)
{
try
{
std::string Driver::evolutionOperatorLookupIds[5] = {
LibUtilities::SessionReader::RegisterEnumValue("EvolutionOperator","Nonlinear" ,eNonlinear),
LibUtilities::SessionReader::RegisterEnumValue("EvolutionOperator","Direct" ,eDirect),
LibUtilities::SessionReader::RegisterEnumValue("EvolutionOperator","Adjoint" ,eAdjoint),
LibUtilities::SessionReader::RegisterEnumValue("EvolutionOperator","TransientGrowth",eTransientGrowth),
LibUtilities::SessionReader::RegisterEnumValue("EvolutionOperator","SkewSymmetric" ,eSkewSymmetric)
};
std::string Driver::evolutionOperatorDef = LibUtilities::SessionReader::RegisterDefaultSolverInfo("EvolutionOperator","Nonlinear");
std::string Driver::driverDefault = LibUtilities::SessionReader::RegisterDefaultSolverInfo("Driver","Standard");
DriverFactory& GetDriverFactory()
// Retrieve the equation system to solve.
ASSERTL0(m_session->DefinesSolverInfo("EqType"),
"EqType SolverInfo tag must be defined.");
std::string vEquation = m_session->GetSolverInfo("EqType");
if (m_session->DefinesSolverInfo("SolverType"))
{
typedef Loki::SingletonHolder<DriverFactory,
Loki::CreateUsingNew,
Loki::NoDestroy > Type;
return Type::Instance();
vEquation = m_session->GetSolverInfo("SolverType");
}
// Check such a module exists for this equation.
ASSERTL0(GetEquationSystemFactory().ModuleExists(vEquation),
"EquationSystem '" + vEquation + "' is not defined.\n"
"Ensure equation name is correct and module is compiled.\n");
/**
*
*/
Driver::Driver(const LibUtilities::SessionReaderSharedPtr pSession)
: m_comm(pSession->GetComm()),
m_session(pSession)
{
}
Driver::~Driver()
{
}
/**
*
*/
void Driver::v_InitObject(ostream &out)
// Retrieve the type of evolution operator to use
/// @todo At the moment this is Navier-Stokes specific - generalise?
m_EvolutionOperator =
m_session->GetSolverInfoAsEnum<EvolutionOperatorType>(
"EvolutionOperator");
m_nequ = ((m_EvolutionOperator == eTransientGrowth ||
m_EvolutionOperator == eAdaptiveSFD) ? 2 : 1);
m_equ = Array<OneD, EquationSystemSharedPtr>(m_nequ);
// Set the AdvectiveType tag and create EquationSystem objects.
switch (m_EvolutionOperator)
{
try
{
// Retrieve the equation system to solve.
ASSERTL0(m_session->DefinesSolverInfo("EqType"),
"EqType SolverInfo tag must be defined.");
std::string vEquation = m_session->GetSolverInfo("EqType");
if (m_session->DefinesSolverInfo("SolverType"))
{
vEquation = m_session->GetSolverInfo("SolverType");
}
// Check such a module exists for this equation.
ASSERTL0(GetEquationSystemFactory().ModuleExists(vEquation),
"EquationSystem '" + vEquation + "' is not defined.\n"
"Ensure equation name is correct and module is compiled.\n");
// Retrieve the type of evolution operator to use
/// @todo At the moment this is Navier-Stokes specific - generalise?
m_EvolutionOperator = m_session->GetSolverInfoAsEnum<EvolutionOperatorType>("EvolutionOperator");
m_nequ = (m_EvolutionOperator == eTransientGrowth ? 2 : 1);
m_equ = Array<OneD, EquationSystemSharedPtr>(m_nequ);
// Set the AdvectiveType tag and create EquationSystem objects.
switch (m_EvolutionOperator)
{
case eNonlinear:
m_session->SetTag("AdvectiveType","Convective");
m_equ[0] = GetEquationSystemFactory().CreateInstance(vEquation, m_session);
break;
case eDirect:
m_session->SetTag("AdvectiveType","Linearised");
m_equ[0] = GetEquationSystemFactory().CreateInstance(vEquation, m_session);
break;
case eAdjoint:
m_session->SetTag("AdvectiveType","Adjoint");
m_equ[0] = GetEquationSystemFactory().CreateInstance(vEquation, m_session);
break;
case eTransientGrowth:
//forward timestepping
m_session->SetTag("AdvectiveType","Linearised");
m_equ[0] = GetEquationSystemFactory().CreateInstance(vEquation, m_session);
//backward timestepping
m_session->SetTag("AdvectiveType","Adjoint");
m_equ[1] = GetEquationSystemFactory().CreateInstance(vEquation, m_session);
break;
case eSkewSymmetric:
m_session->SetTag("AdvectiveType","SkewSymmetric");
m_equ[0] = GetEquationSystemFactory().CreateInstance(vEquation, m_session);
break;
default:
ASSERTL0(false, "Unrecognised evolution operator.");
}
}
catch (int e)
case eNonlinear:
m_session->SetTag("AdvectiveType","Convective");
m_equ[0] = GetEquationSystemFactory().CreateInstance(
vEquation, m_session);
break;
case eDirect:
m_session->SetTag("AdvectiveType","Linearised");
m_equ[0] = GetEquationSystemFactory().CreateInstance(
vEquation, m_session);
break;
case eAdjoint:
m_session->SetTag("AdvectiveType","Adjoint");
m_equ[0] = GetEquationSystemFactory().CreateInstance(
vEquation, m_session);
break;
case eTransientGrowth:
//forward timestepping
m_session->SetTag("AdvectiveType","Linearised");
m_equ[0] = GetEquationSystemFactory().CreateInstance(
vEquation, m_session);
//backward timestepping
m_session->SetTag("AdvectiveType","Adjoint");
m_equ[1] = GetEquationSystemFactory().CreateInstance(
vEquation, m_session);
break;
case eSkewSymmetric:
m_session->SetTag("AdvectiveType","SkewSymmetric");
m_equ[0] = GetEquationSystemFactory().CreateInstance(
vEquation, m_session);
break;
case eAdaptiveSFD:
{
ASSERTL0(e == -1, "No such class class defined.");
out << "An error occurred during driver initialisation." << endl;
// Coupling SFD method and Arnoldi algorithm
// For having 2 equation systems defined into 2 different
// session files (with the mesh into a file named 'session'.gz)
string meshfile;
string LinNSCondFile;
vector<string> LinNSFilename;
meshfile = m_session->GetSessionName() + ".gz";
LinNSCondFile = m_session->GetSessionName();
LinNSCondFile += "_LinNS.xml";
LinNSFilename.push_back(meshfile);
LinNSFilename.push_back(LinNSCondFile);
session_LinNS = LibUtilities::SessionReader::CreateInstance(
0, NULL, LinNSFilename, m_session->GetComm());
//For running stability analysis
session_LinNS->SetTag("AdvectiveType","Linearised");
m_equ[0] = GetEquationSystemFactory().CreateInstance(
vEquation, session_LinNS);
//For running the SFD method on the nonlinear problem
m_session->SetTag("AdvectiveType","Convective");
m_equ[1] = GetEquationSystemFactory().CreateInstance(
vEquation, m_session);
}
}
break;
default:
ASSERTL0(false, "Unrecognised evolution operator.");
Array<OneD, NekDouble> Driver::v_GetRealEvl(void)
{
ASSERTL0(false,"This routine is not valid in this class");
return NullNekDouble1DArray;
}
Array<OneD, NekDouble> Driver::v_GetImagEvl(void)
{
ASSERTL0(false,"This routine is not valid in this class");
return NullNekDouble1DArray;
}
}
}
}
catch (int e)
{
ASSERTL0(e == -1, "No such class class defined.");
out << "An error occurred during driver initialisation." << endl;
}
}
Array<OneD, NekDouble> Driver::v_GetRealEvl(void)
{
ASSERTL0(false,"This routine is not valid in this class");
return NullNekDouble1DArray;
}
Array<OneD, NekDouble> Driver::v_GetImagEvl(void)
{
ASSERTL0(false,"This routine is not valid in this class");
return NullNekDouble1DArray;
}
}
}
......@@ -46,100 +46,104 @@
namespace Nektar
{
namespace SolverUtils
{
class Driver;
/// A shared pointer to a Driver object
typedef boost::shared_ptr<Driver> DriverSharedPtr;
/// Datatype of the NekFactory used to instantiate classes derived from
/// the Driver class.
typedef LibUtilities::NekFactory<
std::string, Driver,
const LibUtilities::SessionReaderSharedPtr&
> DriverFactory;
SOLVER_UTILS_EXPORT DriverFactory& GetDriverFactory();
/// Base class for the development of solvers.
class Driver
{
public:
/// Destructor
virtual ~Driver();
/// Initialise Object
SOLVER_UTILS_EXPORT inline void InitObject(ostream &out = cout);
/// Execute driver
SOLVER_UTILS_EXPORT inline void Execute(ostream &out = cout);
SOLVER_UTILS_EXPORT inline Array<OneD, EquationSystemSharedPtr> GetEqu();
SOLVER_UTILS_EXPORT Array<OneD, NekDouble> GetRealEvl(void);
SOLVER_UTILS_EXPORT Array<OneD, NekDouble> GetImagEvl(void);
protected:
/// Communication object
LibUtilities::CommSharedPtr m_comm;
/// Session reader object
LibUtilities::SessionReaderSharedPtr m_session;
/// Equation system to solve
Array<OneD, EquationSystemSharedPtr> m_equ;
///number of equations
int m_nequ;
///Evolution Operator
enum EvolutionOperatorType m_EvolutionOperator;
/// Initialises EquationSystem class members.
Driver(const LibUtilities::SessionReaderSharedPtr pSession);
SOLVER_UTILS_EXPORT virtual void v_InitObject(ostream &out = cout);
/// Virtual function for solve implementation.
SOLVER_UTILS_EXPORT virtual void v_Execute(ostream &out = cout) = 0;
SOLVER_UTILS_EXPORT virtual Array<OneD, NekDouble> v_GetRealEvl(void);
SOLVER_UTILS_EXPORT virtual Array<OneD, NekDouble> v_GetImagEvl(void);
static std::string evolutionOperatorLookupIds[];
static std::string evolutionOperatorDef;
static std::string driverDefault;
};
inline void Driver::InitObject(ostream &out)
{
v_InitObject(out);
}
inline void Driver::Execute(ostream &out)
{
v_Execute(out);
}
inline Array<OneD, EquationSystemSharedPtr> Driver::GetEqu()
{
return m_equ;
}
inline Array<OneD, NekDouble> Driver::GetRealEvl()
{
return v_GetRealEvl();
}
inline Array<OneD, NekDouble> Driver::GetImagEvl()
{
return v_GetImagEvl();
}
}
namespace SolverUtils
{
class Driver;
/// A shared pointer to a Driver object
typedef boost::shared_ptr<Driver> DriverSharedPtr;
/// Datatype of the NekFactory used to instantiate classes derived from
/// the Driver class.
typedef LibUtilities::NekFactory<std::string, Driver,
const LibUtilities::SessionReaderSharedPtr&> DriverFactory;
SOLVER_UTILS_EXPORT DriverFactory& GetDriverFactory();
/// Base class for the development of solvers.
class Driver
{
public:
/// Destructor
virtual ~Driver();
/// Initialise Object
SOLVER_UTILS_EXPORT inline void InitObject(ostream &out = cout);
/// Execute driver
SOLVER_UTILS_EXPORT inline void Execute(ostream &out = cout);
SOLVER_UTILS_EXPORT inline Array<OneD, EquationSystemSharedPtr> GetEqu();
SOLVER_UTILS_EXPORT Array<OneD, NekDouble> GetRealEvl(void);
SOLVER_UTILS_EXPORT Array<OneD, NekDouble> GetImagEvl(void);
protected:
/// Communication object
LibUtilities::CommSharedPtr m_comm;
/// Session reader object
LibUtilities::SessionReaderSharedPtr m_session;
/// I the Coupling between SFD and arnoldi
LibUtilities::SessionReaderSharedPtr session_LinNS;
/// Equation system to solve
Array<OneD, EquationSystemSharedPtr> m_equ;
///number of equations
int m_nequ;
///Evolution Operator
enum EvolutionOperatorType m_EvolutionOperator;
/// Initialises EquationSystem class members.
Driver(const LibUtilities::SessionReaderSharedPtr pSession);
SOLVER_UTILS_EXPORT virtual void v_InitObject(ostream &out = cout);
/// Virtual function for solve implementation.
SOLVER_UTILS_EXPORT virtual void v_Execute(ostream &out = cout) = 0;
SOLVER_UTILS_EXPORT virtual Array<OneD, NekDouble> v_GetRealEvl(void);
SOLVER_UTILS_EXPORT virtual Array<OneD, NekDouble> v_GetImagEvl(void);
static std::string evolutionOperatorLookupIds[];
static std::string evolutionOperatorDef;
static std::string driverDefault;
};
inline void Driver::InitObject(ostream &out)
{
v_InitObject(out);
}
inline void Driver::Execute(ostream &out)
{
v_Execute(out);
}
inline Array<OneD, EquationSystemSharedPtr> Driver::GetEqu()
{
return m_equ;
}
inline Array<OneD, NekDouble> Driver::GetRealEvl()
{
return v_GetRealEvl();
}
inline Array<OneD, NekDouble> Driver::GetImagEvl()
{
return v_GetImagEvl();
}
}
} //end of namespace
#endif //NEKTAR_SOLVERS_AUXILIARY_ADRBASE_H
......
......@@ -167,15 +167,26 @@ namespace Nektar
*/
void DriverArnoldi::CopyFieldToArnoldiArray(Array<OneD, NekDouble> &array)
{
Array<OneD, MultiRegions::ExpListSharedPtr> fields;
fields = m_equ[m_nequ-1]->UpdateFields();
if (m_EvolutionOperator == eAdaptiveSFD)
{
//This matters for the Adaptive SFD method
//because m_equ[1] is the nonlinear problem with non homogeneous BCs.
fields = m_equ[0]->UpdateFields();
}
else
{
fields = m_equ[m_nequ-1]->UpdateFields();
}
for (int k = 0; k < m_nfields; ++k)
{
int nq = fields[0]->GetNcoeffs();
Vmath::Vcopy(nq, &fields[k]->GetCoeffs()[0], 1, &array[k*nq], 1);
fields[k]->SetPhysState(false);
}
};
......
......@@ -41,98 +41,98 @@
namespace Nektar
{
namespace SolverUtils
{
class DriverModifiedArnoldi: public DriverArnoldi
{
public:
friend class MemoryManager<DriverModifiedArnoldi>;
/// Creates an instance of this class
static DriverSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession)
namespace SolverUtils
{
DriverSharedPtr p = MemoryManager<DriverModifiedArnoldi>
::AllocateSharedPtr(pSession);
p->InitObject();
return p;
class DriverModifiedArnoldi: public DriverArnoldi
{
public:
friend class MemoryManager<DriverModifiedArnoldi>;
/// Creates an instance of this class
static DriverSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession)
{
DriverSharedPtr p = MemoryManager<DriverModifiedArnoldi>
::AllocateSharedPtr(pSession);
p->InitObject();
return p;
}
///Name of the class
static std::string className;
protected:
/// Constructor
DriverModifiedArnoldi(const LibUtilities::SessionReaderSharedPtr pSession);
/// Destructor
virtual ~DriverModifiedArnoldi();
/// Virtual function for initialisation implementation.
virtual void v_InitObject(ostream &out = cout );
/// Virtual function for solve implementation.
virtual void v_Execute(ostream &out = cout);
private:
/// Generates a new vector in the sequence by applying the linear operator.
void EV_update(Array<OneD, NekDouble> &src,
Array<OneD, NekDouble> &tgt);
/// Generates the upper Hessenberg matrix H and computes its eigenvalues.
void EV_small(Array<OneD, Array<OneD, NekDouble> > &Kseq,
const int ntot,
const Array<OneD, NekDouble> &alpha,
const int kdim,
Array<OneD, NekDouble> &zvec,
Array<OneD, NekDouble> &wr,
Array<OneD, NekDouble> &wi,
NekDouble &resnorm);
/// Tests for convergence of eigenvalues of H.
int EV_test(const int itrn,
const int kdim,
Array<OneD, NekDouble> &zvec,
Array<OneD, NekDouble> &wr,
Array<OneD, NekDouble> &wi,
const NekDouble resnorm,
const int nvec,
ofstream &evlout,
NekDouble &resid0);
/// Sorts a sequence of eigenvectors/eigenvalues by magnitude.
void EV_sort(Array<OneD, NekDouble> &evec,
Array<OneD, NekDouble> &wr,
Array<OneD, NekDouble> &wi,
Array<OneD, NekDouble> &test,
const int dim);
void EV_post(Array<OneD, Array<OneD, NekDouble> > &Tseq,
Array<OneD, Array<OneD, NekDouble> > &Kseq,
const int ntot,
const int kdim,
const int nvec,
Array<OneD, NekDouble> &zvec,
Array<OneD, NekDouble> &wr,
Array<OneD, NekDouble> &wi,
const int icon);
void EV_big(Array<OneD, Array<OneD, NekDouble> > &bvecs,