#include <DriverSteadyState.h>

Inheritance diagram for Nektar::SolverUtils::DriverSteadyState:

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...

SOLVER_UTILS_EXPORT	~DriverSteadyState () override=default
	Destructor. More...

SOLVER_UTILS_EXPORT void	v_InitObject (std::ostream &out=std::cout) override
	Initialises EquationSystem class members. More...

SOLVER_UTILS_EXPORT void	v_Execute (std::ostream &out=std::cout) override
	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...

	~DriverModifiedArnoldi () override=default
	Destructor. More...

void	v_InitObject (std::ostream &out=std::cout) override
	Virtual function for initialisation implementation. More...

void	v_Execute (std::ostream &out=std::cout) override
	Virtual function for solve implementation. More...

Protected Member Functions inherited from Nektar::SolverUtils::DriverArnoldi
	DriverArnoldi (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
	Constructor. More...

	~DriverArnoldi () override=default
	Destructor. More...

void	v_InitObject (std::ostream &out=std::cout) override
	Virtual function for initialisation implementation. More...

void	v_Execute (std::ostream &out=std::cout) override
	Virtual function for solve implementation. 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	GetMaskInfo (std::vector< std::vector< LibUtilities::EquationSharedPtr > > &selectedDomains, std::set< int > &unselectedVariables)

SOLVER_UTILS_EXPORT void	ArnoldiSummary (std::ostream &out)

Protected Member Functions inherited from Nektar::SolverUtils::Driver
	Driver (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
	Initialises EquationSystem class members. More...

virtual SOLVER_UTILS_EXPORT void	v_InitObject (std::ostream &out=std::cout)
	Virtual function for initialisation implementation. More...

virtual SOLVER_UTILS_EXPORT void	v_Execute (std::ostream &out=std::cout)=0
	Virtual function for solve implementation. More...

Private 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)

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

Static Private Attributes
static std::string	driverLookupId

Friends
class	MemoryManager< DriverSteadyState >

Additional Inherited Members
Public Member Functions inherited from Nektar::SolverUtils::DriverArnoldi
Array< OneD, NekDouble >	GetRealEvl (void)

Array< OneD, NekDouble >	GetImagEvl (void)

Public Member Functions inherited from Nektar::SolverUtils::Driver
virtual	~Driver ()=default
	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 ()

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 iterations. 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
	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...

Static Protected Attributes inherited from Nektar::SolverUtils::DriverModifiedArnoldi
static std::string	driverLookupId

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

static std::string	evolutionOperatorDef

static std::string	driverDefault

Detailed Description

Definition at line 44 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 56 of file DriverSteadyState.cpp.

    : DriverModifiedArnoldi(pSession, pGraph)
{
}

◆ ~DriverSteadyState()

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

overrideprotecteddefault

Destructor.

Member Function Documentation

◆ ComputeOptimization()

void Nektar::SolverUtils::DriverSteadyState::ComputeOptimization ( )

private

Definition at line 368 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
	)

private

Encapsulated SFD method

Definition at line 348 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
	)

private

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

Definition at line 592 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]->GetTime()
             << ";\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]->GetTime()
               << "\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 50 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
	)

private

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 481 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(), Nektar::UnitTests::d(), and tinysimd::sqrt().

Referenced by GradientDescentMethod().

◆ GradientDescentMethod()

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

private

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 398 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 ( )

private

Definition at line 644 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
	)

private

Definition at line 522 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
	)

private

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 333 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 )

overrideprotectedvirtual

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 74 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" ||
        m_session->GetSolverInfo("EqType") == "NavierStokesImplicitCFE")
    {
        // 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 )

overrideprotectedvirtual

Initialises EquationSystem class members.

Reimplemented from Nektar::SolverUtils::DriverModifiedArnoldi.

Definition at line 66 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 111 of file DriverSteadyState.h.

Referenced by PrintSummarySFD(), and v_Execute().

◆ AdaptiveTOL

NekDouble Nektar::SolverUtils::DriverSteadyState::AdaptiveTOL

private

Definition at line 110 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 61 of file DriverSteadyState.h.

◆ cpuTime

NekDouble Nektar::SolverUtils::DriverSteadyState::cpuTime

private

Definition at line 89 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

◆ Diff_q1_q0

NekDouble Nektar::SolverUtils::DriverSteadyState::Diff_q1_q0

private

Definition at line 106 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ Diff_q_qBar

NekDouble Nektar::SolverUtils::DriverSteadyState::Diff_q_qBar

private

Definition at line 105 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ driverLookupId

string Nektar::SolverUtils::DriverSteadyState::driverLookupId

staticprivate

Initial value:

=

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

Nektar::LibUtilities::SessionReader::RegisterEnumValue

static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)

Registers an enumeration value.

Definition: SessionReader.h:560

Definition at line 144 of file DriverSteadyState.h.

◆ elapsed

NekDouble Nektar::SolverUtils::DriverSteadyState::elapsed

private

Definition at line 91 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

◆ FlowPartiallyConverged

bool Nektar::SolverUtils::DriverSteadyState::FlowPartiallyConverged

private

For adaptive SFD method.

Definition at line 109 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ FrequencyEV

NekDouble Nektar::SolverUtils::DriverSteadyState::FrequencyEV

private

Definition at line 114 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ GrowthRateEV

NekDouble Nektar::SolverUtils::DriverSteadyState::GrowthRateEV

private

Definition at line 113 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ M11

NekDouble Nektar::SolverUtils::DriverSteadyState::M11

private

Definition of the SFD operator.

Definition at line 100 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

◆ M12

NekDouble Nektar::SolverUtils::DriverSteadyState::M12

private

Definition at line 101 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

◆ M21

NekDouble Nektar::SolverUtils::DriverSteadyState::M21

private

Definition at line 102 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

◆ M22

NekDouble Nektar::SolverUtils::DriverSteadyState::M22

private

Definition at line 103 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

◆ m_Check

int Nektar::SolverUtils::DriverSteadyState::m_Check

private

Definition at line 81 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ m_Check_BaseFlow

int Nektar::SolverUtils::DriverSteadyState::m_Check_BaseFlow

private

Definition at line 82 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ m_checksteps

int Nektar::SolverUtils::DriverSteadyState::m_checksteps

private

Definition at line 85 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 95 of file DriverSteadyState.h.

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

◆ m_dt

NekDouble Nektar::SolverUtils::DriverSteadyState::m_dt

private

Definition at line 97 of file DriverSteadyState.h.

Referenced by SetSFDOperator(), and v_Execute().

◆ m_file

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

private

Definition at line 92 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

◆ m_infosteps

int Nektar::SolverUtils::DriverSteadyState::m_infosteps

private

Definition at line 84 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ m_NonConvergingStepsCounter

int Nektar::SolverUtils::DriverSteadyState::m_NonConvergingStepsCounter

private

Definition at line 112 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ m_stepCounter

int Nektar::SolverUtils::DriverSteadyState::m_stepCounter

private

Definition at line 80 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

◆ m_X

NekDouble Nektar::SolverUtils::DriverSteadyState::m_X

private

Definition at line 96 of file DriverSteadyState.h.

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

◆ NumVar_SFD

int Nektar::SolverUtils::DriverSteadyState::NumVar_SFD

private

Definition at line 86 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

◆ timer

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

private

Definition at line 88 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

◆ TOL

NekDouble Nektar::SolverUtils::DriverSteadyState::TOL

private

Definition at line 83 of file DriverSteadyState.h.

Referenced by PrintSummarySFD(), and v_Execute().

◆ totalTime

NekDouble Nektar::SolverUtils::DriverSteadyState::totalTime

private

Definition at line 90 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Private Member Functions

Private Attributes

Static Private Attributes

Friends

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ DriverSteadyState()

◆ ~DriverSteadyState()

Member Function Documentation

◆ ComputeOptimization()

◆ ComputeSFD()

◆ ConvergenceHistory()

◆ create()

◆ EvalEV_ScalarSFD()

◆ GradientDescentMethod()

◆ PrintSummarySFD()

◆ ReadEVfile()

◆ SetSFDOperator()

◆ v_Execute()

◆ v_InitObject()

Friends And Related Function Documentation

◆ MemoryManager< DriverSteadyState >

Member Data Documentation

◆ AdaptiveTime

◆ AdaptiveTOL

◆ className

◆ cpuTime

◆ Diff_q1_q0

◆ Diff_q_qBar

◆ driverLookupId

◆ elapsed

◆ FlowPartiallyConverged

◆ FrequencyEV

◆ GrowthRateEV

◆ M11

◆ M12

◆ M21

◆ M22

◆ m_Check

◆ m_Check_BaseFlow

◆ m_checksteps

◆ m_Delta

◆ m_dt

◆ m_file

◆ m_infosteps

◆ m_NonConvergingStepsCounter

◆ m_stepCounter

◆ m_X

◆ NumVar_SFD

◆ timer

◆ TOL

◆ totalTime