Nektar::SolverUtils::DriverSteadyState Class Reference

#include <DriverSteadyState.h>

Inheritance diagram for Nektar::SolverUtils::DriverSteadyState:

Public Member Functions
void	PrintSummarySFD ()
void	ConvergenceHistory (const Array< OneD, const Array< OneD, NekDouble > > &qBar1, const Array< OneD, const Array< OneD, NekDouble > > &q0, NekDouble &MaxNormDiff_q_qBar, NekDouble &MaxNormDiff_q1_q0)
void	ComputeSFD (const int i, const Array< OneD, const Array< OneD, NekDouble > > &q0, const Array< OneD, const Array< OneD, NekDouble > > &qBar0, Array< OneD, Array< OneD, NekDouble > > &q1, Array< OneD, Array< OneD, NekDouble > > &qBar1)
void	EvalEV_ScalarSFD (const NekDouble &X_input, const NekDouble &Delta_input, const std::complex< NekDouble > &alpha, NekDouble &MaxEV)
void	GradientDescentMethod (const std::complex< NekDouble > &alpha, NekDouble &X_output, NekDouble &Delta_output)
void	ReadEVfile (int &KrylovSubspaceDim, NekDouble &growthEV, NekDouble &frequencyEV)
void	ComputeOptimization ()
void	SetSFDOperator (const NekDouble X_input, const NekDouble Delta_input)
Public Member Functions inherited from Nektar::SolverUtils::DriverArnoldi
SOLVER_UTILS_EXPORT void	ArnoldiSummary (std::ostream &out)
SOLVER_UTILS_EXPORT const Array< OneD, const NekDouble > &	GetMaskCoeff () const
SOLVER_UTILS_EXPORT const Array< OneD, const NekDouble > &	GetMaskPhys () const
Public Member Functions inherited from Nektar::SolverUtils::Driver
virtual	~Driver ()
	Destructor. More...
SOLVER_UTILS_EXPORT void	InitObject (std::ostream &out=std::cout)
	Initialise Object. More...
SOLVER_UTILS_EXPORT void	Execute (std::ostream &out=std::cout)
	Execute driver. More...
SOLVER_UTILS_EXPORT Array< OneD, EquationSystemSharedPtr >	GetEqu ()
SOLVER_UTILS_EXPORT Array< OneD, NekDouble >	GetRealEvl (void)
SOLVER_UTILS_EXPORT Array< OneD, NekDouble >	GetImagEvl (void)

Static Public Member Functions
static DriverSharedPtr	create (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
	Creates an instance of this class. More...
Static Public Member Functions inherited from Nektar::SolverUtils::DriverModifiedArnoldi
static DriverSharedPtr	create (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
	Creates an instance of this class. More...

Static Public Attributes
static std::string	className
	Name of the class. More...
Static Public Attributes inherited from Nektar::SolverUtils::DriverModifiedArnoldi
static std::string	className
	Name of the class. More...

Protected Member Functions
SOLVER_UTILS_EXPORT	DriverSteadyState (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
	Constructor. More...
virtual SOLVER_UTILS_EXPORT	~DriverSteadyState ()
	Destructor. More...
virtual SOLVER_UTILS_EXPORT void	v_InitObject (std::ostream &out=std::cout)
	Second-stage initialisation. More...
virtual SOLVER_UTILS_EXPORT void	v_Execute (std::ostream &out=std::cout)
	Virtual function for solve implementation. More...
Protected Member Functions inherited from Nektar::SolverUtils::DriverModifiedArnoldi
	DriverModifiedArnoldi (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
	Constructor. More...
virtual	~DriverModifiedArnoldi ()
	Destructor. More...
Protected Member Functions inherited from Nektar::SolverUtils::DriverArnoldi
	DriverArnoldi (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
	Constructor. More...
virtual	~DriverArnoldi ()
	Destructor. More...
void	CopyArnoldiArrayToField (Array< OneD, NekDouble > &array)
	Copy Arnoldi storage to fields. More...
void	CopyFieldToArnoldiArray (Array< OneD, NekDouble > &array)
	Copy fields to Arnoldi storage. More...
void	CopyFwdToAdj ()
	Copy the forward field to the adjoint system in transient growth calculations. More...
void	WriteFld (std::string file, std::vector< Array< OneD, NekDouble > > coeffs)
	Write coefficients to file. More...
void	WriteFld (std::string file, Array< OneD, NekDouble > coeffs)
void	WriteEvs (std::ostream &evlout, const int k, const NekDouble real, const NekDouble imag, NekDouble resid=NekConstants::kNekUnsetDouble, bool DumpInverse=true)
void	MaskInit ()
	init mask More...
void	GetUnmaskFunction (std::vector< std::vector< LibUtilities::EquationSharedPtr > > &unmaskfun)
virtual Array< OneD, NekDouble >	v_GetRealEvl (void)
virtual Array< OneD, NekDouble >	v_GetImagEvl (void)
Protected Member Functions inherited from Nektar::SolverUtils::Driver
	Driver (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
	Initialises EquationSystem class members. More...

Static Protected Attributes
static std::string	driverLookupId
Static Protected Attributes inherited from Nektar::SolverUtils::Driver
static std::string	evolutionOperatorLookupIds []
static std::string	evolutionOperatorDef
static std::string	driverDefault

Private Attributes
int	m_stepCounter
int	m_Check
int	m_Check_BaseFlow
NekDouble	TOL
int	m_infosteps
int	m_checksteps
int	NumVar_SFD
LibUtilities::Timer	timer
NekDouble	cpuTime
NekDouble	totalTime
NekDouble	elapsed
std::ofstream	m_file
NekDouble	m_Delta
	Definition of the SFD parameters: More...
NekDouble	m_X
NekDouble	m_dt
NekDouble	M11
	Definition of the SFD operator. More...
NekDouble	M12
NekDouble	M21
NekDouble	M22
NekDouble	Diff_q_qBar
NekDouble	Diff_q1_q0
bool	FlowPartiallyConverged
	For adaptive SFD method. More...
NekDouble	AdaptiveTOL
NekDouble	AdaptiveTime
int	m_NonConvergingStepsCounter
NekDouble	GrowthRateEV
NekDouble	FrequencyEV

Friends
class	MemoryManager< DriverSteadyState >

Additional Inherited Members
Protected Attributes inherited from Nektar::SolverUtils::DriverArnoldi
int	m_kdim
int	m_nvec
	Dimension of Krylov subspace. More...
int	m_nits
	Number of vectors to test. More...
NekDouble	m_evtol
	Maxmum number of iterations. More...
NekDouble	m_period
	Tolerance of iteratiosn. More...
bool	m_timeSteppingAlgorithm
	Period of time stepping algorithm. More...
int	m_infosteps
	underlying operator is time stepping More...
int	m_nfields
	interval to dump information if required. More...
NekDouble	m_realShift
NekDouble	m_imagShift
int	m_negatedOp
Array< OneD, NekDouble >	m_real_evl
	Operator in solve call is negated. More...
Array< OneD, NekDouble >	m_imag_evl
bool	m_useMask
Array< OneD, NekDouble >	m_maskCoeffs
Array< OneD, NekDouble >	m_maskPhys
Protected Attributes inherited from Nektar::SolverUtils::Driver
LibUtilities::CommSharedPtr	m_comm
	Communication object. More...
LibUtilities::SessionReaderSharedPtr	m_session
	Session reader object. More...
LibUtilities::SessionReaderSharedPtr	session_LinNS
	I the Coupling between SFD and arnoldi. More...
SpatialDomains::MeshGraphSharedPtr	m_graph
	MeshGraph object. More...
Array< OneD, EquationSystemSharedPtr >	m_equ
	Equation system to solve. More...
int	m_nequ
	number of equations More...
enum EvolutionOperatorType	m_EvolutionOperator
	Evolution Operator. More...

Detailed Description

Definition at line 48 of file DriverSteadyState.h.

Constructor & Destructor Documentation

DriverSteadyState()

Nektar::SolverUtils::DriverSteadyState::DriverSteadyState ( const LibUtilities::SessionReaderSharedPtr  pSession, const SpatialDomains::MeshGraphSharedPtr  pGraph  )

protected

Constructor.

Definition at line 59 of file DriverSteadyState.cpp.

    : DriverModifiedArnoldi(pSession, pGraph)
{
}

~DriverSteadyState()

Nektar::SolverUtils::DriverSteadyState::~DriverSteadyState ( )

protectedvirtual

Destructor.

Definition at line 70 of file DriverSteadyState.cpp.

{
}

Member Function Documentation

ComputeOptimization()

void Nektar::SolverUtils::DriverSteadyState::ComputeOptimization

(

)

Definition at line 372 of file DriverSteadyState.cpp.

{
    NekDouble growthEV(0.0);
    NekDouble frequencyEV(0.0);
 
    // m_kdim is the dimension of Krylov subspace (defined in the xml file and
    // used in DriverArnoldi.cpp)
    ReadEVfile(m_kdim, growthEV, frequencyEV);
 
    cout << "\n\tgrowthEV = " << growthEV << endl;
    cout << "\tfrequencyEV = " << frequencyEV << endl;
 
    complex<NekDouble> ApproxEV = polar(exp(growthEV), frequencyEV);
 
    NekDouble X_new = m_X;
    NekDouble Delta_new = m_Delta;
 
    GradientDescentMethod(ApproxEV, X_new, Delta_new);
 
    m_X = X_new;
    m_Delta = Delta_new;
 
    SetSFDOperator(m_X, m_Delta);
}

References GradientDescentMethod(), m_Delta, Nektar::SolverUtils::DriverArnoldi::m_kdim, m_X, ReadEVfile(), and SetSFDOperator().

Referenced by v_Execute().

ComputeSFD()

void Nektar::SolverUtils::DriverSteadyState::ComputeSFD	(	const int	i,
		const Array< OneD, const Array< OneD, NekDouble > > &	q0,
		const Array< OneD, const Array< OneD, NekDouble > > &	qBar0,
		Array< OneD, Array< OneD, NekDouble > > &	q1,
		Array< OneD, Array< OneD, NekDouble > > &	qBar1
	)

Encapsulated SFD method

Definition at line 352 of file DriverSteadyState.cpp.

{
    q1[i]    = Array<OneD, NekDouble> (m_equ[m_nequ - 1]->GetTotPoints(),0.0);
    qBar1[i] = Array<OneD, NekDouble> (m_equ[m_nequ - 1]->GetTotPoints(),0.0);
 
    ///Encapsulated SFD method
    Vmath::Svtvp(q1[i].size(), M11, q0[i],    1, q1[i], 1, q1[i], 1 );
    Vmath::Svtvp(q1[i].size(), M12, qBar0[i], 1, q1[i], 1, q1[i], 1 );
 
    Vmath::Svtvp(qBar1[i].size(), M21, q0[i],    1, qBar1[i], 1,
                                                            qBar1[i], 1 );
    Vmath::Svtvp(qBar1[i].size(), M22, qBar0[i], 1, qBar1[i], 1,
                                                            qBar1[i], 1 );
}

References M11, M12, M21, M22, Nektar::SolverUtils::Driver::m_equ, Nektar::SolverUtils::Driver::m_nequ, and Vmath::Svtvp().

Referenced by v_Execute().

ConvergenceHistory()

void Nektar::SolverUtils::DriverSteadyState::ConvergenceHistory	(	const Array< OneD, const Array< OneD, NekDouble > > &	qBar1,
		const Array< OneD, const Array< OneD, NekDouble > > &	q0,
		NekDouble &	MaxNormDiff_q_qBar,
		NekDouble &	MaxNormDiff_q1_q0
	)

This routine evaluates |q-qBar|_inf (and |q1-q0|_inf) and writes the values in "ConvergenceHistory.txt"

Definition at line 596 of file DriverSteadyState.cpp.

{
    Array<OneD, NekDouble > NormDiff_q_qBar(NumVar_SFD, 1.0);
    Array<OneD, NekDouble > NormDiff_q1_q0(NumVar_SFD, 1.0);
    MaxNormDiff_q_qBar=0.0;
    MaxNormDiff_q1_q0=0.0;
 
    for(int i = 0; i < NumVar_SFD; ++i)
    {
        NormDiff_q_qBar[i] = m_equ[m_nequ - 1]->LinfError(i, qBar1[i]);
        NormDiff_q1_q0[i] = m_equ[m_nequ - 1]->LinfError(i, q0[i]);
 
        if (MaxNormDiff_q_qBar < NormDiff_q_qBar[i])
        {
            MaxNormDiff_q_qBar = NormDiff_q_qBar[i];
        }
        if (MaxNormDiff_q1_q0 < NormDiff_q1_q0[i])
        {
            MaxNormDiff_q1_q0 = NormDiff_q1_q0[i];
        }
    }
 
    timer.Stop();
    elapsed  = timer.TimePerTest(1);
    cpuTime += elapsed;
    totalTime += elapsed;
 
    if (m_comm->GetRank() == 0)
    {
        cout << "SFD - Step: " <<  left <<  m_stepCounter+1
             << ";\tTime: " << left << m_equ[m_nequ - 1]->GetFinalTime()
             << ";\tCPU time = " << left << cpuTime << " s"
             << ";\tTot time = " << left << totalTime << " s"
             << ";\tX = " << left << m_X
             << ";\tDelta = " << left << m_Delta
             << ";\t|q-qBar|inf = " << left << MaxNormDiff_q_qBar << endl;
        std::ofstream m_file( "ConvergenceHistory.txt", std::ios_base::app);
        m_file << m_stepCounter+1 << "\t"
               << m_equ[m_nequ - 1]->GetFinalTime() << "\t"
               << totalTime << "\t"
               << MaxNormDiff_q_qBar << "\t"
               << MaxNormDiff_q1_q0 << endl;
        m_file.close();
    }
 
    cpuTime = 0.0;
    timer.Start();
}

References cpuTime, elapsed, Nektar::SolverUtils::Driver::m_comm, m_Delta, Nektar::SolverUtils::Driver::m_equ, m_file, Nektar::SolverUtils::Driver::m_nequ, m_stepCounter, m_X, NumVar_SFD, Nektar::LibUtilities::Timer::Start(), Nektar::LibUtilities::Timer::Stop(), Nektar::LibUtilities::Timer::TimePerTest(), timer, and totalTime.

Referenced by v_Execute().

create()

static DriverSharedPtr Nektar::SolverUtils::DriverSteadyState::create ( const LibUtilities::SessionReaderSharedPtr &  pSession, const SpatialDomains::MeshGraphSharedPtr &  pGraph  )

inlinestatic

Creates an instance of this class.

Definition at line 54 of file DriverSteadyState.h.

    {
        DriverSharedPtr p = MemoryManager<DriverSteadyState>
            ::AllocateSharedPtr(pSession, pGraph);
        p->InitObject();
        return p;
    }

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.

EvalEV_ScalarSFD()

void Nektar::SolverUtils::DriverSteadyState::EvalEV_ScalarSFD	(	const NekDouble &	X_input,
		const NekDouble &	Delta_input,
		const std::complex< NekDouble > &	alpha,
		NekDouble &	MaxEV
	)

This routine evaluates the maximum eigenvalue of the SFD system when applied to the 1D model u(n+1) = alpha*u(n)

Definition at line 488 of file DriverSteadyState.cpp.

{
    NekDouble A11 = ( 1.0 + X_input * Delta_input *
            exp(-(X_input + 1.0/Delta_input)) )/(1.0 + X_input*Delta_input);
    NekDouble A12 = ( X_input*Delta_input - X_input * Delta_input *
            exp(-(X_input + 1.0/Delta_input)) )/(1.0 + X_input*Delta_input);
    NekDouble A21 = ( 1.0 - 1.0 *
            exp(-(X_input + 1.0/Delta_input)) )/(1 + X_input*Delta_input);
    NekDouble A22 = ( X_input*Delta_input + 1.0 *
            exp(-(X_input + 1.0/Delta_input)) )/(1.0 + X_input*Delta_input);
 
    complex<NekDouble> B11 = alpha;
    NekDouble B12 = 0.0;
    NekDouble B21 = 0.0;
    NekDouble B22 = 1.0;
 
    complex<NekDouble> a = A11*B11 + A12*B21;
    NekDouble b = A11*B12 + A12*B22;
    complex<NekDouble> c = A21*B11 + A22*B21;
    NekDouble d = A21*B12 + A22*B22;
 
    complex<NekDouble> delt = sqrt((a-d)*(a-d) + 4.0*b*c);
    complex<NekDouble> lambda_1 = 0.5*(a+d + delt);
    complex<NekDouble> lambda_2 = 0.5*(a+d - delt);
 
    NekDouble NORM_1 = abs(lambda_1);
    NekDouble NORM_2 = abs(lambda_2);
 
    MaxEV = max(NORM_1, NORM_2);
}

References tinysimd::abs(), and tinysimd::sqrt().

Referenced by GradientDescentMethod().

GradientDescentMethod()

void Nektar::SolverUtils::DriverSteadyState::GradientDescentMethod	(	const std::complex< NekDouble > &	alpha,
		NekDouble &	X_output,
		NekDouble &	Delta_output
	)

This routine implements a gradient descent method to find the parameters X end Delta which give the minimum eigenlavue of the SFD problem applied to the scalar case u(n+1) = \alpha*u(n).

Definition at line 402 of file DriverSteadyState.cpp.

{
    cout << "\n\tWe enter the Gradient Descent Method [...]" << endl;
    bool OptParmFound = false;
    bool Descending = true;
    NekDouble X_input = X_output;
    NekDouble Delta_input = Delta_output;
 
    NekDouble X_init = X_output;
    NekDouble Delta_init = Delta_output;
    int stepCounter(0);
 
    NekDouble F0(0.0);
    NekDouble F0x(0.0);
    NekDouble F0y(0.0);
    NekDouble F1(0.0);
    NekDouble dx = 0.00000001;
    NekDouble dirx(0.0);
    NekDouble diry(0.0);
    NekDouble s(0.0);
    NekDouble CurrentMin = 1.0;
 
    while (OptParmFound == false)
    {
        Descending = true;
        EvalEV_ScalarSFD(X_input, Delta_input, alpha, F0);
        EvalEV_ScalarSFD(X_input + dx, Delta_input, alpha, F0x);
        EvalEV_ScalarSFD(X_input, Delta_input + dx, alpha, F0y);
        dirx =  (F0 - F0x);
        diry =  (F0 - F0y);
        s = abs(0.000001/dirx);
        X_output = X_input + s*dirx;
        Delta_output = Delta_input + s*diry;
        F1 = F0;
 
        while (Descending == true)
        {
            CurrentMin = F1;
            X_input = X_output;
            Delta_input = Delta_output;
            EvalEV_ScalarSFD(X_output, Delta_output, alpha, F1);
 
            if (F1 > CurrentMin)
            {
                Descending = false;
            }
            else
            {
                s = s + s*0.01;
                X_output = X_input + s*dirx;
                Delta_output = Delta_input + s*diry;
            }
 
            if(stepCounter > 9999999)
            {
                //We are stuck in this loop..
                //Then we restart it with different initail conditions
                Descending = false;
                X_input = X_init;
                Delta_init = Delta_init + Delta_init*0.1;
                Delta_input = Delta_init;
                stepCounter = 0;
            }
            stepCounter++;
        }
 
        if (abs(F0-F1) < dx)
        {
            cout << "\tThe Gradient Descent Method has converged!" << endl;
            EvalEV_ScalarSFD(X_output, Delta_output, alpha, F1);
            cout << "\n\tThe updated parameters are: X_tilde = " << X_output
                 << " and Delta_tilde = " << Delta_output << endl;
            OptParmFound = true;
        }
 
    }
}

References tinysimd::abs(), and EvalEV_ScalarSFD().

Referenced by ComputeOptimization(), and v_Execute().

PrintSummarySFD()

void Nektar::SolverUtils::DriverSteadyState::PrintSummarySFD

(

)

Definition at line 650 of file DriverSteadyState.cpp.

{
    cout << "\n====================================="
            "=================================="
         << endl;
    cout << "Parameters for the SFD method:" << endl;
    cout << "\tControl Coefficient: X = " << m_X << endl;
    cout << "\tFilter Width: Delta = " << m_Delta << endl;
    cout << "The simulation is stopped when |q-qBar|inf < " << TOL << endl;
    if (m_EvolutionOperator == eAdaptiveSFD)
    {
        cout << "\nWe use the adaptive SFD method:" << endl;
        cout << "  The parameters are updated every " << AdaptiveTime
             << " time units;" << endl;
        cout << "  until |q-qBar|inf becomes smaller than " << AdaptiveTOL
             << endl;
    }
    cout << "====================================="
            "==================================\n" << endl;
}

References AdaptiveTime, AdaptiveTOL, Nektar::SolverUtils::eAdaptiveSFD, m_Delta, Nektar::SolverUtils::Driver::m_EvolutionOperator, m_X, and TOL.

Referenced by v_Execute().

ReadEVfile()

void Nektar::SolverUtils::DriverSteadyState::ReadEVfile	(	int &	KrylovSubspaceDim,
		NekDouble &	growthEV,
		NekDouble &	frequencyEV
	)

Definition at line 524 of file DriverSteadyState.cpp.

{
    // This routine reads the .evl file written by the Arnoldi algorithm
    // (written in September 2014)
    std::string   line;
    int           NumLinesInFile = 0;
    std::string   EVfileName = m_session->GetSessionName() +  ".evl";
    std::ifstream EVfile(EVfileName.c_str());
    ASSERTL0(EVfile.good(), "Cannot open .evl file.");
 
    if(EVfile)
    {
        // This block counts the total number of lines of the .evl file
        // We keep going util we reach the end of the file
        while(getline(EVfile, line))
        {
            ++NumLinesInFile;
        }
 
        // It may happen that the Stability method that have produced the .elv
        // file converges in less than m_kdim iterations. In this case,
        // KrylovSubspaceDim has to be changed here
        if(NumLinesInFile < KrylovSubspaceDim*2.0
                                + KrylovSubspaceDim*(KrylovSubspaceDim+1.0)/2.0)
        {
            for(int i = 1; i <= KrylovSubspaceDim; ++i)
            {
                if(NumLinesInFile == i*2.0 + i*(i+1.0)/2.0)
                {
                    KrylovSubspaceDim = i;
                }
            }
        }
 
        // go back to the beginning of the file
        EVfile.clear();
        EVfile.seekg(0, ios::beg);
 
        // We now want to go to the line where the most unstable eigenlavue was
        // written
        for(int i = 0; i < (NumLinesInFile - KrylovSubspaceDim); ++i)
        {
            std::getline(EVfile, line);
            cout << "Discard line: " << line << endl;
        }
 
        std::vector<std::string> tokens;
        std::getline(EVfile, line);
        boost::algorithm::split(tokens, line,
                boost::is_any_of("\t "),boost::token_compress_on);
 
        ASSERTL0(tokens.size() >= 5,
                 "Unexpected formatting of .evl file while reading line:\n"
                 + line);
        growthEV = boost::lexical_cast<NekDouble>(tokens[4]);
        frequencyEV = boost::lexical_cast<NekDouble>(tokens[5]);
    }
    else
    {
        cout << "An error occurred when opening the .evl file" << endl;
    }
    EVfile.close();
}

References ASSERTL0, and Nektar::SolverUtils::Driver::m_session.

Referenced by ComputeOptimization().

SetSFDOperator()

void Nektar::SolverUtils::DriverSteadyState::SetSFDOperator	(	const NekDouble	X_input,
		const NekDouble	Delta_input
	)

This routine defines the encapsulated SFD operator with first-order splitting and exact resolution of the second subproblem. (See http://scitation.aip.org/content/aip/journal/pof2/26/3/10.1063/1.4867482 for details)

Definition at line 339 of file DriverSteadyState.cpp.

{
    NekDouble X = X_input*m_dt;
    NekDouble Delta = Delta_input/m_dt;
    NekDouble coeff = 1.0/(1.0 + X*Delta);
    M11 = coeff*(1.0 + X*Delta*exp(-(X + 1.0/Delta)));
    M12 = coeff*(X*Delta*(1.0 - exp(-(X + 1.0/Delta))));
    M21 = coeff*(1.0 - exp(-(X + 1.0/Delta)));
    M22 = coeff*(X*Delta + exp(-(X + 1.0/Delta)));
}

References M11, M12, M21, M22, and m_dt.

Referenced by ComputeOptimization(), and v_Execute().

v_Execute()

void Nektar::SolverUtils::DriverSteadyState::v_Execute ( std::ostream &  out = std::cout )

protectedvirtual

Virtual function for solve implementation.

Definition of variables used in this algorithm

Call the Navier-Stokes solver for one time step

Copy the current flow field into q0

Apply the linear operator to the outcome of the solver

Update qBar

Copy the output of the SFD method into the current flow field

Loop for the adaptive SFD method

We save the final solution into a .fld file

Reimplemented from Nektar::SolverUtils::DriverModifiedArnoldi.

Definition at line 84 of file DriverSteadyState.cpp.

{
    // With a loop over "DoSolve", this Driver implements the
    // "encapsulated" Selective Frequency Damping method(SFD)
    // to find the steady state of a flow above the critical Reynolds
    // number.
    m_equ[m_nequ - 1]->PrintSummary(out);
 
    m_session->LoadParameter("IO_InfoSteps", m_infosteps, 1000);
    m_session->LoadParameter("IO_CheckSteps", m_checksteps, 100000);
    m_session->LoadParameter("ControlCoeff",m_X, 1);
    m_session->LoadParameter("FilterWidth", m_Delta, 2);
 
    // To evaluate optimum SFD parameters if growth rate provided in the
    // xml file
    m_session->LoadParameter("GrowthRateEV", GrowthRateEV, 0.0);
 
    // To evaluate optimum SFD parameters if frequency provided in the xml
    // file
    m_session->LoadParameter("FrequencyEV", FrequencyEV, 0.0);
 
    // Determine when SFD method is converged
    m_session->LoadParameter("TOL", TOL, 1.0e-08);
 
    // Used only for the Adaptive SFD method
    m_session->LoadParameter("AdaptiveTOL", AdaptiveTOL, 1.0e-02);
 
    // Used only for the Adaptive SFD method
    m_session->LoadParameter("AdaptiveTime", AdaptiveTime, 50.0*m_Delta);
 
    if (m_comm->GetRank() == 0)
    {
        PrintSummarySFD();
    }
 
    timer.Start();
 
    // Definition of shared pointer (used only for the Adaptive SFD method)
    AdvectionSystemSharedPtr A = m_equ[0]->as<AdvectionSystem>();
 
    NumVar_SFD = m_equ[m_nequ - 1]->UpdateFields()[0]->GetCoordim(0);
    // SFD to run for incompressible case with scalar field
    if (m_session->GetSolverInfo("EqType")=="UnsteadyNavierStokes"||
        m_session->GetSolverInfo("EqType")=="SteadyNavierStokes")
    {
        int nConvectiveFields = m_session->GetVariables().size();
        if (boost::iequals(m_session->GetVariable(nConvectiveFields-1),"p"))
        {
            nConvectiveFields -= 1;
        }
        NumVar_SFD = nConvectiveFields;
    }
    // Condition necessary to run SFD for the compressible case
    if (m_session->GetSolverInfo("EqType") == "EulerCFE" ||
        m_session->GetSolverInfo("EqType") == "NavierStokesCFE")
    {
        // Number of variables for the compressible equations
        NumVar_SFD += 2;
    }
    if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
    {
        if (m_session->GetSolverInfo("HOMOGENEOUS") == "1D")
        {
            NumVar_SFD += 1;
        }
    }
 
    // We store the time step
    m_dt = m_equ[m_nequ - 1]->GetTimeStep();
 
    // Evaluate optimum SFD parameters if dominent EV given by xml file
    if (GrowthRateEV != 0.0 && FrequencyEV != 0.0)
    {
        cout << "Besed on the dominant EV given in the xml file,"
             << "a 1D model is used to evaluate the optumum parameters"
             << "of the SFD method:" << endl;
        complex<NekDouble> EV = polar(exp(GrowthRateEV), FrequencyEV);
        GradientDescentMethod(EV, m_X, m_Delta);
    }
 
 
    // We set up the elements of the operator of the encapsulated
    // formulation of the selective frequencive damping method
    SetSFDOperator(m_X, m_Delta);
 
    // m_steps is set to 1. Then "m_equ[m_nequ - 1]->DoSolve()" will run
    // for only one time step
    m_equ[m_nequ - 1]->SetSteps(1);
    ofstream m_file("ConvergenceHistory.txt", ios::out | ios::trunc);
 
    Array<OneD, Array<OneD, NekDouble> > q0(NumVar_SFD);
    Array<OneD, Array<OneD, NekDouble> > q1(NumVar_SFD);
    Array<OneD, Array<OneD, NekDouble> > qBar0(NumVar_SFD);
    Array<OneD, Array<OneD, NekDouble> > qBar1(NumVar_SFD);
 
    for(int i = 0; i < NumVar_SFD; ++i)
    {
        q0[i] = Array<OneD, NekDouble> (m_equ[m_nequ-1]->GetTotPoints(),
                                        0.0); //q0 is initialised
        qBar0[i] = Array<OneD, NekDouble> (m_equ[m_nequ-1]->GetTotPoints(),
                                        0.0);
        m_equ[m_nequ - 1]->CopyFromPhysField(i, qBar0[i]);
    }
 
    ///Definition of variables used in this algorithm
    m_stepCounter               = 0;
    m_Check                     = 0;
    m_Check_BaseFlow            = 1;
    m_NonConvergingStepsCounter = 0;
    Diff_q_qBar                 = 1.0;
    Diff_q1_q0                  = 1.0;
    cpuTime                     = 0.0;
    elapsed                     = 0.0;
    totalTime                   = 0.0;
    FlowPartiallyConverged      = false;
 
    while (max(Diff_q_qBar, Diff_q1_q0) > TOL)
    {
        ///Call the Navier-Stokes solver for one time step
        m_equ[m_nequ - 1]->DoSolve();
 
        for(int i = 0; i < NumVar_SFD; ++i)
        {
            ///Copy the current flow field into q0
            m_equ[m_nequ - 1]->CopyFromPhysField(i, q0[i]);
 
            ///Apply the linear operator to the outcome of the solver
            ComputeSFD(i, q0, qBar0, q1, qBar1);
 
            ///Update qBar
            qBar0[i] = qBar1[i];
 
            ///Copy the output of the SFD method into the current flow field
            m_equ[m_nequ - 1]->CopyToPhysField(i, q1[i]);
        }
 
        if(m_infosteps && !((m_stepCounter+1)%m_infosteps))
        {
            ConvergenceHistory(qBar1, q0, Diff_q_qBar, Diff_q1_q0);
 
            ///Loop for the adaptive SFD method
            if (m_EvolutionOperator    == eAdaptiveSFD &&
                FlowPartiallyConverged == false)
            {
 
                m_NonConvergingStepsCounter++;
 
                if (Diff_q_qBar < AdaptiveTOL)
                {
                    if (m_comm->GetRank() == 0)
                    {
                    cout << "\n\t The SFD method is converging: we compute "
                         << "stability analysis using the 'partially "
                         << "converged' steady state as base flow:\n" << endl;
                    }
 
                    m_equ[m_nequ - 1]->Checkpoint_BaseFlow(m_Check_BaseFlow);
                    m_Check_BaseFlow++;
 
                    A->GetAdvObject()->SetBaseFlow(q0,m_equ[0]->UpdateFields());
                    DriverModifiedArnoldi::v_Execute(out);
 
                    if (m_comm->GetRank() == 0)
                    {
                        // Compute the update of the SFD parameters only on
                        // one processor
                        ComputeOptimization();
                    }
                    else
                    {
                        // On all the other processors, the parameters are set
                        // to 0
                        m_X = 0;
                        m_Delta = 0;
                    }
                    // The we give to all the processors the value of X and
                    // Delta of the first processor
                    m_comm->AllReduce(m_X,     Nektar::LibUtilities::ReduceSum);
                    m_comm->AllReduce(m_Delta, Nektar::LibUtilities::ReduceSum);
 
                    FlowPartiallyConverged = true;
                }
                else if (m_NonConvergingStepsCounter * m_dt * m_infosteps
                            >= AdaptiveTime)
                {
                    if (m_comm->GetRank() == 0)
                    {
                        cout << "\n\t We compute stability analysis using"
                             << " the current flow field as base flow:\n"
                             << endl;
                    }
 
                    m_equ[m_nequ - 1]->Checkpoint_BaseFlow(m_Check_BaseFlow);
                    m_Check_BaseFlow++;
 
                    A->GetAdvObject()->SetBaseFlow(q0,m_equ[0]->UpdateFields());
                    DriverModifiedArnoldi::v_Execute(out);
 
                    if (m_comm->GetRank() == 0)
                    {
                        // Compute the update of the SFD parameters only on
                        // one processor
                        ComputeOptimization();
                    }
                    else
                    {
                        // On all the other processors, the parameters are set
                        // to 0
                        m_X = 0;
                        m_Delta = 0;
                    }
                    // The we give to all the processors the value of X and
                    // Delta of the first processor
                    m_comm->AllReduce(m_X, Nektar::LibUtilities::ReduceSum);
                    m_comm->AllReduce(m_Delta, Nektar::LibUtilities::ReduceSum);
 
                    m_NonConvergingStepsCounter = 0;
                }
            }
        }
 
        if(m_checksteps && m_stepCounter&&(!((m_stepCounter+1)%m_checksteps)))
        {
            m_Check++;
            m_equ[m_nequ - 1]->Checkpoint_Output(m_Check);
        }
        m_stepCounter++;
    }
 
    m_file.close();
 
    ///We save the final solution into a .fld file
    m_equ[m_nequ - 1]->Output();
 
    for(int j = 0; j < m_equ[m_nequ - 1]->GetNvariables(); ++j)
    {
        NekDouble vL2Error = m_equ[m_nequ - 1]->L2Error(j,false);
        NekDouble vLinfError = m_equ[m_nequ - 1]->LinfError(j);
        if (m_comm->GetRank() == 0)
        {
            out << "L 2 error (variable " << m_equ[m_nequ - 1]->GetVariable(j)
                << ") : " << vL2Error << endl;
            out << "L inf error (variable " << m_equ[m_nequ - 1]->GetVariable(j)
                << ") : " << vLinfError << endl;
        }
    }
}

References AdaptiveTime, AdaptiveTOL, ComputeOptimization(), ComputeSFD(), ConvergenceHistory(), cpuTime, Diff_q1_q0, Diff_q_qBar, Nektar::SolverUtils::eAdaptiveSFD, elapsed, FlowPartiallyConverged, FrequencyEV, GradientDescentMethod(), GrowthRateEV, m_Check, m_Check_BaseFlow, m_checksteps, Nektar::SolverUtils::Driver::m_comm, m_Delta, m_dt, Nektar::SolverUtils::Driver::m_equ, Nektar::SolverUtils::Driver::m_EvolutionOperator, m_file, m_infosteps, Nektar::SolverUtils::Driver::m_nequ, m_NonConvergingStepsCounter, Nektar::SolverUtils::Driver::m_session, m_stepCounter, m_X, NumVar_SFD, PrintSummarySFD(), Nektar::LibUtilities::ReduceSum, SetSFDOperator(), Nektar::LibUtilities::Timer::Start(), timer, TOL, totalTime, and Nektar::SolverUtils::DriverModifiedArnoldi::v_Execute().

v_InitObject()

void Nektar::SolverUtils::DriverSteadyState::v_InitObject ( std::ostream &  out = std::cout )

protectedvirtual

Second-stage initialisation.

Reimplemented from Nektar::SolverUtils::DriverModifiedArnoldi.

Definition at line 78 of file DriverSteadyState.cpp.

{
    DriverModifiedArnoldi::v_InitObject(out);
}

References Nektar::SolverUtils::DriverModifiedArnoldi::v_InitObject().

Friends And Related Function Documentation

MemoryManager< DriverSteadyState >

friend class MemoryManager< DriverSteadyState >

friend

Definition at line 1 of file DriverSteadyState.h.

Member Data Documentation

AdaptiveTime

NekDouble Nektar::SolverUtils::DriverSteadyState::AdaptiveTime

private

Definition at line 154 of file DriverSteadyState.h.

Referenced by PrintSummarySFD(), and v_Execute().

AdaptiveTOL

NekDouble Nektar::SolverUtils::DriverSteadyState::AdaptiveTOL

private

Definition at line 153 of file DriverSteadyState.h.

Referenced by PrintSummarySFD(), and v_Execute().

className

string Nektar::SolverUtils::DriverSteadyState::className

static

Initial value:

= GetDriverFactory().RegisterCreatorFunction(
                "SteadyState", DriverSteadyState::create)

Name of the class.

Definition at line 65 of file DriverSteadyState.h.

cpuTime

NekDouble Nektar::SolverUtils::DriverSteadyState::cpuTime

private

Definition at line 132 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

Diff_q1_q0

NekDouble Nektar::SolverUtils::DriverSteadyState::Diff_q1_q0

private

Definition at line 149 of file DriverSteadyState.h.

Referenced by v_Execute().

Diff_q_qBar

NekDouble Nektar::SolverUtils::DriverSteadyState::Diff_q_qBar

private

Definition at line 148 of file DriverSteadyState.h.

Referenced by v_Execute().

driverLookupId

string Nektar::SolverUtils::DriverSteadyState::driverLookupId

staticprotected

Initial value:

= LibUtilities::SessionReader::RegisterEnumValue(
                "Driver","SteadyState",0)

Definition at line 120 of file DriverSteadyState.h.

elapsed

NekDouble Nektar::SolverUtils::DriverSteadyState::elapsed

private

Definition at line 134 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

FlowPartiallyConverged

bool Nektar::SolverUtils::DriverSteadyState::FlowPartiallyConverged

private

For adaptive SFD method.

Definition at line 152 of file DriverSteadyState.h.

Referenced by v_Execute().

FrequencyEV

NekDouble Nektar::SolverUtils::DriverSteadyState::FrequencyEV

private

Definition at line 157 of file DriverSteadyState.h.

Referenced by v_Execute().

GrowthRateEV

NekDouble Nektar::SolverUtils::DriverSteadyState::GrowthRateEV

private

Definition at line 156 of file DriverSteadyState.h.

Referenced by v_Execute().

M11

NekDouble Nektar::SolverUtils::DriverSteadyState::M11

private

Definition of the SFD operator.

Definition at line 143 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

M12

NekDouble Nektar::SolverUtils::DriverSteadyState::M12

private

Definition at line 144 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

M21

NekDouble Nektar::SolverUtils::DriverSteadyState::M21

private

Definition at line 145 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

M22

NekDouble Nektar::SolverUtils::DriverSteadyState::M22

private

Definition at line 146 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

m_Check

int Nektar::SolverUtils::DriverSteadyState::m_Check

private

Definition at line 124 of file DriverSteadyState.h.

Referenced by v_Execute().

m_Check_BaseFlow

int Nektar::SolverUtils::DriverSteadyState::m_Check_BaseFlow

private

Definition at line 125 of file DriverSteadyState.h.

Referenced by v_Execute().

m_checksteps

int Nektar::SolverUtils::DriverSteadyState::m_checksteps

private

Definition at line 128 of file DriverSteadyState.h.

Referenced by v_Execute().

m_Delta

NekDouble Nektar::SolverUtils::DriverSteadyState::m_Delta

private

Definition of the SFD parameters:

Definition at line 138 of file DriverSteadyState.h.

Referenced by ComputeOptimization(), ConvergenceHistory(), PrintSummarySFD(), and v_Execute().

m_dt

NekDouble Nektar::SolverUtils::DriverSteadyState::m_dt

private

Definition at line 140 of file DriverSteadyState.h.

Referenced by SetSFDOperator(), and v_Execute().

m_file

std::ofstream Nektar::SolverUtils::DriverSteadyState::m_file

private

Definition at line 135 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

m_infosteps

int Nektar::SolverUtils::DriverSteadyState::m_infosteps

private

Definition at line 127 of file DriverSteadyState.h.

Referenced by v_Execute().

m_NonConvergingStepsCounter

int Nektar::SolverUtils::DriverSteadyState::m_NonConvergingStepsCounter

private

Definition at line 155 of file DriverSteadyState.h.

Referenced by v_Execute().

m_stepCounter

int Nektar::SolverUtils::DriverSteadyState::m_stepCounter

private

Definition at line 123 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

m_X

NekDouble Nektar::SolverUtils::DriverSteadyState::m_X

private

Definition at line 139 of file DriverSteadyState.h.

Referenced by ComputeOptimization(), ConvergenceHistory(), PrintSummarySFD(), and v_Execute().

NumVar_SFD

int Nektar::SolverUtils::DriverSteadyState::NumVar_SFD

private

Definition at line 129 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

timer

LibUtilities::Timer Nektar::SolverUtils::DriverSteadyState::timer

private

Definition at line 131 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

TOL

NekDouble Nektar::SolverUtils::DriverSteadyState::TOL

private

Definition at line 126 of file DriverSteadyState.h.

Referenced by PrintSummarySFD(), and v_Execute().

totalTime

NekDouble Nektar::SolverUtils::DriverSteadyState::totalTime

private

Definition at line 133 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().