Nektar::SolverUtils::DriverSteadyState Class Reference

#include <DriverSteadyState.h>

Inheritance diagram for Nektar::SolverUtils::DriverSteadyState:

Collaboration 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)
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)
	Creates an instance of this class. More...
Static Public Member Functions inherited from Nektar::SolverUtils::DriverModifiedArnoldi
static DriverSharedPtr	create (const LibUtilities::SessionReaderSharedPtr &pSession)
	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)
	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)
	Constructor. More...
virtual	~DriverModifiedArnoldi ()
	Destructor. More...
Protected Member Functions inherited from Nektar::SolverUtils::DriverArnoldi
	DriverArnoldi (const LibUtilities::SessionReaderSharedPtr pSession)
	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)
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)
	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
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
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...
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 49 of file DriverSteadyState.h.

Constructor & Destructor Documentation

Nektar::SolverUtils::DriverSteadyState::DriverSteadyState ( const LibUtilities::SessionReaderSharedPtr  pSession )

protected

Constructor.

Definition at line 56 of file DriverSteadyState.cpp.

    : DriverModifiedArnoldi(pSession)
{
}

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

protectedvirtual

Destructor.

Definition at line 66 of file DriverSteadyState.cpp.

{
}

Member Function Documentation

void Nektar::SolverUtils::DriverSteadyState::ComputeOptimization

(

)

Definition at line 357 of file DriverSteadyState.cpp.

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

Referenced by v_Execute().

{
    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);
}

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 337 of file DriverSteadyState.cpp.

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

Referenced by v_Execute().

{
    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].num_elements(), M11, q0[i],    1, q1[i], 1, q1[i], 1 );
    Vmath::Svtvp(q1[i].num_elements(), M12, qBar0[i], 1, q1[i], 1, q1[i], 1 );

    Vmath::Svtvp(qBar1[i].num_elements(), M21, q0[i],    1, qBar1[i], 1,
                                                            qBar1[i], 1 );
    Vmath::Svtvp(qBar1[i].num_elements(), M22, qBar0[i], 1, qBar1[i], 1,
                                                            qBar1[i], 1 );
}

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 581 of file DriverSteadyState.cpp.

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::Timer::Start(), Nektar::Timer::Stop(), Nektar::Timer::TimePerTest(), timer, and totalTime.

Referenced by v_Execute().

{
    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();
}

static DriverSharedPtr Nektar::SolverUtils::DriverSteadyState::create ( const LibUtilities::SessionReaderSharedPtr &  pSession )

inlinestatic

Creates an instance of this class.

Definition at line 55 of file DriverSteadyState.h.

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

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

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 473 of file DriverSteadyState.cpp.

Referenced by GradientDescentMethod().

{
    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);
}

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) = *u(n).

Definition at line 387 of file DriverSteadyState.cpp.

References EvalEV_ScalarSFD().

Referenced by ComputeOptimization(), and v_Execute().

{
    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;
        }

    }
}

void Nektar::SolverUtils::DriverSteadyState::PrintSummarySFD

(

)

Definition at line 635 of file DriverSteadyState.cpp.

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

Referenced by v_Execute().

{
    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;
}

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

Definition at line 509 of file DriverSteadyState.cpp.

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

Referenced by ComputeOptimization().

{
    // 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();
}

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 324 of file DriverSteadyState.cpp.

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

Referenced by ComputeOptimization(), and v_Execute().

{
    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)));
}

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 80 of file DriverSteadyState.cpp.

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::Timer::Start(), timer, TOL, totalTime, and Nektar::SolverUtils::DriverModifiedArnoldi::v_Execute().

{
    // With a loop over "DoSolve", this Driver implements the
    // "encaplulated" Selective Frequency Damping method
    // 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>();

    // Condition necessary to run SFD for the compressible case
    NumVar_SFD = m_equ[m_nequ - 1]->UpdateFields()[0]->GetCoordim(0);
    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;
        }
    }
}

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

protectedvirtual

Second-stage initialisation.

Reimplemented from Nektar::SolverUtils::DriverModifiedArnoldi.

Definition at line 74 of file DriverSteadyState.cpp.

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

{
    DriverModifiedArnoldi::v_InitObject(out);
}

Friends And Related Function Documentation

friend class MemoryManager< DriverSteadyState >

friend

Definition at line 52 of file DriverSteadyState.h.

Member Data Documentation

NekDouble Nektar::SolverUtils::DriverSteadyState::AdaptiveTime

private

Definition at line 152 of file DriverSteadyState.h.

Referenced by PrintSummarySFD(), and v_Execute().

NekDouble Nektar::SolverUtils::DriverSteadyState::AdaptiveTOL

private

Definition at line 151 of file DriverSteadyState.h.

Referenced by PrintSummarySFD(), and v_Execute().

string Nektar::SolverUtils::DriverSteadyState::className

static

Initial value:

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

Name of the class.

Definition at line 64 of file DriverSteadyState.h.

NekDouble Nektar::SolverUtils::DriverSteadyState::cpuTime

private

Definition at line 130 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

NekDouble Nektar::SolverUtils::DriverSteadyState::Diff_q1_q0

private

Definition at line 147 of file DriverSteadyState.h.

Referenced by v_Execute().

NekDouble Nektar::SolverUtils::DriverSteadyState::Diff_q_qBar

private

Definition at line 146 of file DriverSteadyState.h.

Referenced by v_Execute().

string Nektar::SolverUtils::DriverSteadyState::driverLookupId

staticprotected

Initial value:

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

Definition at line 118 of file DriverSteadyState.h.

NekDouble Nektar::SolverUtils::DriverSteadyState::elapsed

private

Definition at line 132 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

bool Nektar::SolverUtils::DriverSteadyState::FlowPartiallyConverged

private

For adaptive SFD method.

Definition at line 150 of file DriverSteadyState.h.

Referenced by v_Execute().

NekDouble Nektar::SolverUtils::DriverSteadyState::FrequencyEV

private

Definition at line 155 of file DriverSteadyState.h.

Referenced by v_Execute().

NekDouble Nektar::SolverUtils::DriverSteadyState::GrowthRateEV

private

Definition at line 154 of file DriverSteadyState.h.

Referenced by v_Execute().

NekDouble Nektar::SolverUtils::DriverSteadyState::M11

private

Definition of the SFD operator.

Definition at line 141 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

NekDouble Nektar::SolverUtils::DriverSteadyState::M12

private

Definition at line 142 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

NekDouble Nektar::SolverUtils::DriverSteadyState::M21

private

Definition at line 143 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

NekDouble Nektar::SolverUtils::DriverSteadyState::M22

private

Definition at line 144 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

int Nektar::SolverUtils::DriverSteadyState::m_Check

private

Definition at line 122 of file DriverSteadyState.h.

Referenced by v_Execute().

int Nektar::SolverUtils::DriverSteadyState::m_Check_BaseFlow

private

Definition at line 123 of file DriverSteadyState.h.

Referenced by v_Execute().

int Nektar::SolverUtils::DriverSteadyState::m_checksteps

private

Definition at line 126 of file DriverSteadyState.h.

Referenced by v_Execute().

NekDouble Nektar::SolverUtils::DriverSteadyState::m_Delta

private

Definition of the SFD parameters:

Definition at line 136 of file DriverSteadyState.h.

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

NekDouble Nektar::SolverUtils::DriverSteadyState::m_dt

private

Definition at line 138 of file DriverSteadyState.h.

Referenced by SetSFDOperator(), and v_Execute().

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

private

Definition at line 133 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

int Nektar::SolverUtils::DriverSteadyState::m_infosteps

private

Definition at line 125 of file DriverSteadyState.h.

Referenced by v_Execute().

int Nektar::SolverUtils::DriverSteadyState::m_NonConvergingStepsCounter

private

Definition at line 153 of file DriverSteadyState.h.

Referenced by v_Execute().

int Nektar::SolverUtils::DriverSteadyState::m_stepCounter

private

Definition at line 121 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

NekDouble Nektar::SolverUtils::DriverSteadyState::m_X

private

Definition at line 137 of file DriverSteadyState.h.

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

int Nektar::SolverUtils::DriverSteadyState::NumVar_SFD

private

Definition at line 127 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

Timer Nektar::SolverUtils::DriverSteadyState::timer

private

Definition at line 129 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

NekDouble Nektar::SolverUtils::DriverSteadyState::TOL

private

Definition at line 124 of file DriverSteadyState.h.

Referenced by PrintSummarySFD(), and v_Execute().

NekDouble Nektar::SolverUtils::DriverSteadyState::totalTime

private

Definition at line 131 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().