#include <EulerExponentialTimeIntegrationSchemes.h>

Inheritance diagram for Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme:

Public Member Functions
	EulerExponentialTimeIntegrationScheme (std::string variant, size_t order, std::vector< NekDouble > freeParams)

	~EulerExponentialTimeIntegrationScheme () override

LUE void	SetExponentialCoefficients (Array< OneD, std::complex< NekDouble > > &Lambda)

Public Member Functions inherited from Nektar::LibUtilities::TimeIntegrationSchemeGLM
LUE void	InitializeSecondaryData (TimeIntegrationAlgorithmGLM *phase, NekDouble deltaT) const

Public Member Functions inherited from Nektar::LibUtilities::TimeIntegrationScheme
LUE std::string	GetFullName () const

LUE std::string	GetName () const

LUE std::string	GetVariant () const

LUE size_t	GetOrder () const

LUE std::vector< NekDouble >	GetFreeParams ()

LUE TimeIntegrationSchemeType	GetIntegrationSchemeType ()

LUE NekDouble	GetTimeStability () const

LUE size_t	GetNumIntegrationPhases ()

LUE const TripleArray &	GetSolutionVector () const
	Gets the solution vector of the ODE. More...

LUE TripleArray &	UpdateSolutionVector ()

LUE void	SetSolutionVector (const size_t Offset, const DoubleArray &y)
	Sets the solution vector of the ODE. More...

LUE void	InitializeScheme (const NekDouble deltaT, ConstDoubleArray &y_0, const NekDouble time, const TimeIntegrationSchemeOperators &op)
	Explicit integration of an ODE. More...

LUE ConstDoubleArray &	TimeIntegrate (const size_t timestep, const NekDouble delta_t)

LUE void	print (std::ostream &os) const

LUE void	printFull (std::ostream &os) const

Static Public Member Functions
static TimeIntegrationSchemeSharedPtr	create (std::string variant, size_t order, std::vector< NekDouble > freeParams)

static LUE void	SetupSchemeData (TimeIntegrationAlgorithmGLMSharedPtr &phase, std::string variant, size_t order)

Static Public Attributes
static std::string	className

Protected Member Functions
LUE std::string	v_GetName () const override

LUE std::string	v_GetFullName () const override

LUE NekDouble	v_GetTimeStability () const override

void	v_InitializeSecondaryData (TimeIntegrationAlgorithmGLM *phase, NekDouble deltaT) const override

Protected Member Functions inherited from Nektar::LibUtilities::TimeIntegrationSchemeGLM
LUE std::string	v_GetVariant () const override

LUE size_t	v_GetOrder () const override

LUE std::vector< NekDouble >	v_GetFreeParams () const override

LUE TimeIntegrationSchemeType	v_GetIntegrationSchemeType () const override

LUE size_t	v_GetNumIntegrationPhases () const override

LUE const TripleArray &	v_GetSolutionVector () const override

LUE TripleArray &	v_UpdateSolutionVector () override

LUE void	v_SetSolutionVector (const size_t Offset, const DoubleArray &y) override

LUE void	v_InitializeScheme (const NekDouble deltaT, ConstDoubleArray &y_0, const NekDouble time, const TimeIntegrationSchemeOperators &op) override
	Worker method to initialize the integration scheme. More...

LUE ConstDoubleArray &	v_TimeIntegrate (const size_t timestep, const NekDouble delta_t) override
	Worker method that actually does the time integration. More...

virtual LUE void	v_InitializeSecondaryData (TimeIntegrationAlgorithmGLM *phase, NekDouble deltaT) const

LUE void	v_print (std::ostream &os) const override
	Worker method to print details on the integration scheme. More...

LUE void	v_printFull (std::ostream &os) const override

LUE	TimeIntegrationSchemeGLM (std::string variant, size_t order, std::vector< NekDouble > freeParams)

	~TimeIntegrationSchemeGLM () override

Protected Member Functions inherited from Nektar::LibUtilities::TimeIntegrationScheme
virtual LUE std::string	v_GetFullName () const

virtual LUE std::string	v_GetName () const =0

virtual LUE std::string	v_GetVariant () const =0

virtual LUE size_t	v_GetOrder () const =0

virtual LUE std::vector< NekDouble >	v_GetFreeParams () const =0

virtual LUE TimeIntegrationSchemeType	v_GetIntegrationSchemeType () const =0

virtual LUE NekDouble	v_GetTimeStability () const =0

virtual LUE size_t	v_GetNumIntegrationPhases () const =0

virtual LUE const TripleArray &	v_GetSolutionVector () const =0

virtual LUE TripleArray &	v_UpdateSolutionVector ()=0

virtual LUE void	v_SetSolutionVector (const size_t Offset, const DoubleArray &y)=0

virtual LUE void	v_InitializeScheme (const NekDouble deltaT, ConstDoubleArray &y_0, const NekDouble time, const TimeIntegrationSchemeOperators &op)=0

virtual LUE ConstDoubleArray &	v_TimeIntegrate (const size_t timestep, const NekDouble delta_t)=0

virtual LUE void	v_print (std::ostream &os) const =0

virtual LUE void	v_printFull (std::ostream &os) const =0

LUE	TimeIntegrationScheme (std::string variant, size_t order, std::vector< NekDouble > freeParams)

LUE	TimeIntegrationScheme (const TimeIntegrationScheme &in)=delete

virtual	~TimeIntegrationScheme ()=default

Private Member Functions
NekDouble	factorial (size_t n) const

std::complex< NekDouble >	phi_function (const size_t order, const std::complex< NekDouble > z) const

Additional Inherited Members
Protected Attributes inherited from Nektar::LibUtilities::TimeIntegrationSchemeGLM
TimeIntegrationAlgorithmGLMVector	m_integration_phases

TimeIntegrationSolutionGLMSharedPtr	m_solVector

Detailed Description

Definition at line 61 of file EulerExponentialTimeIntegrationSchemes.h.

Constructor & Destructor Documentation

◆ EulerExponentialTimeIntegrationScheme()

Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::EulerExponentialTimeIntegrationScheme	(	std::string	variant,
		size_t	order,
		std::vector< NekDouble >	freeParams
	)

inline

Definition at line 64 of file EulerExponentialTimeIntegrationSchemes.h.

        : TimeIntegrationSchemeGLM(variant, order, freeParams)
    {
        // Currently up to 4th order is implemented because the number
        // of steps is the same as the order.
        ASSERTL1(variant == "Lawson" || variant == "Norsett",
                 "EulerExponential Time integration scheme unknown variant: " +
                     variant + ". Must be 'Lawson' or 'Norsett'.");
 
        ASSERTL1(((variant == "Lawson" && 1 <= order && order <= 1) ||
                  (variant == "Norsett" && 1 <= order && order <= 4)),
                 "EulerExponential Time integration scheme bad order, "
                 "Lawson (1) or Norsett (1-4): " +
                     std::to_string(order));
 
        m_integration_phases = TimeIntegrationAlgorithmGLMVector(order);
 
        // Currently the next lowest order is used to seed the current
        // order. This is not correct but is an okay approximation.
        for (size_t n = 0; n < order; ++n)
        {
            m_integration_phases[n] = TimeIntegrationAlgorithmGLMSharedPtr(
                new TimeIntegrationAlgorithmGLM(this));
 
            EulerExponentialTimeIntegrationScheme::SetupSchemeData(
                m_integration_phases[n], variant, n + 1);
        }
    }

References ASSERTL1, Nektar::LibUtilities::TimeIntegrationSchemeGLM::m_integration_phases, and SetupSchemeData().

◆ ~EulerExponentialTimeIntegrationScheme()

Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::~EulerExponentialTimeIntegrationScheme ( )

inlineoverride

Definition at line 94 of file EulerExponentialTimeIntegrationSchemes.h.

95 {

96 }

Member Function Documentation

◆ create()

static TimeIntegrationSchemeSharedPtr Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::create	(	std::string	variant,
		size_t	order,
		std::vector< NekDouble >	freeParams
	)

inlinestatic

Definition at line 98 of file EulerExponentialTimeIntegrationSchemes.h.

    {
        TimeIntegrationSchemeSharedPtr p =
            MemoryManager<EulerExponentialTimeIntegrationScheme>::
                AllocateSharedPtr(variant, order, freeParams);
 
        return p;
    }

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

◆ factorial()

NekDouble Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::factorial ( size_t n ) const

inlineprivate

Definition at line 425 of file EulerExponentialTimeIntegrationSchemes.h.

    {
        return (n == 1 || n == 0) ? 1 : n * factorial(n - 1);
    }

References factorial().

Referenced by factorial(), and phi_function().

◆ phi_function()

std::complex< NekDouble > Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::phi_function	(	const size_t	order,
		const std::complex< NekDouble >	z
	)		const

inlineprivate

Phi function

Central to the implementation of exponential integrators is the evaluation of exponential-like functions, commonly denoted by phi (φ) functions. It is convenient to define then as φ0(z) = e^z, in which case the functions obey the recurrence relation.

0: exp(z); 1: (exp(z) - 1.0) / (z); 2: (exp(z) - z - 1.0) / (z * z);

Definition at line 430 of file EulerExponentialTimeIntegrationSchemes.h.

    {
        /**
         * \brief Phi function
         *
         * Central to the implementation of exponential integrators is
         * the evaluation of exponential-like functions, commonly
         * denoted by phi (φ) functions. It is convenient to define
         * then as φ0(z) = e^z, in which case the functions obey the
         * recurrence relation.
 
         * 0:  exp(z);
         * 1: (exp(z)     - 1.0) / (z);
         * 2: (exp(z) - z - 1.0) / (z * z);
         */
 
        if (z == 0.0)
        {
            return 1.0 / factorial(order);
        }
 
        if (order == 0)
        {
            return exp(z);
        }
        else
        {
            return (phi_function(order - 1, z) - 1.0 / factorial(order - 1)) /
                   z;
        }
 
        return 0;
    }

References factorial(), phi_function(), and Nektar::UnitTests::z().

Referenced by phi_function(), and v_InitializeSecondaryData().

◆ SetExponentialCoefficients()

LUE void Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::SetExponentialCoefficients ( Array< OneD, std::complex< NekDouble > > & Lambda )

inline

Lambda Matrix Assumption, parameter Lambda.

The one-dimensional Lambda matrix is a diagonal matrix thus values are non zero if and only i=j. As such, the diagonal Lambda values are stored in an array of complex numbers.

Definition at line 180 of file EulerExponentialTimeIntegrationSchemes.h.

    {
        ASSERTL0(!m_integration_phases.empty(), "No scheme")
 
        /**
         * \brief Lambda Matrix Assumption, parameter Lambda.
         *
         * The one-dimensional Lambda matrix is a diagonal
         * matrix thus values are non zero if and only i=j. As such,
         * the diagonal Lambda values are stored in an array of
         * complex numbers.
         */
 
        // Assume that each phase is an exponential integrator.
        for (size_t i = 0; i < m_integration_phases.size(); i++)
        {
            m_integration_phases[i]->m_L = Lambda;
 
            // Anytime the coefficents are updated reset the nVars to
            // be assured that the exponential matrices are
            // recalculated (e.g. the number of variables may remain
            // the same but the coefficients have changed).
            m_integration_phases[i]->m_lastNVars = 0;
        }
    }

References ASSERTL0, and Nektar::LibUtilities::TimeIntegrationSchemeGLM::m_integration_phases.

◆ SetupSchemeData()

static LUE void Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::SetupSchemeData	(	TimeIntegrationAlgorithmGLMSharedPtr &	phase,
		std::string	variant,
		size_t	order
	)

inlinestatic

Definition at line 110 of file EulerExponentialTimeIntegrationSchemes.h.

    {
        phase->m_schemeType = eExponential;
        phase->m_variant    = variant;
        phase->m_order      = order;
        phase->m_name       = variant + std::string("EulerExponentialOrder") +
                        std::to_string(phase->m_order);
 
        // Parameters for the compact 1 step implementation.
        phase->m_numstages = 1;
        phase->m_numsteps  = phase->m_order; // Okay up to 4th order.
 
        phase->m_A = Array<OneD, Array<TwoD, NekDouble>>(1);
        phase->m_B = Array<OneD, Array<TwoD, NekDouble>>(1);
 
        phase->m_A[0] =
            Array<TwoD, NekDouble>(phase->m_numstages, phase->m_numstages, 0.0);
        phase->m_B[0] =
            Array<TwoD, NekDouble>(phase->m_numsteps, phase->m_numstages, 0.0);
 
        phase->m_U =
            Array<TwoD, NekDouble>(phase->m_numstages, phase->m_numsteps, 0.0);
        phase->m_V =
            Array<TwoD, NekDouble>(phase->m_numsteps, phase->m_numsteps, 0.0);
 
        // Coefficients
 
        // B Phi function for first row first column
        phase->m_B[0][0][0] = 1.0 / phase->m_order; // phi_func(0 or 1)
 
        // B evaluation value shuffling second row first column
        if (phase->m_order > 1)
        {
            phase->m_B[0][1][0] = 1.0; // constant 1
        }
 
        // U Curent time step evaluation first row first column
        phase->m_U[0][0] = 1.0; // constant 1
 
        // V Phi function for first row first column
        phase->m_V[0][0] = 1.0; // phi_func(0)
 
        // V Phi function for first row additional columns
        for (size_t n = 1; n < phase->m_order; ++n)
        {
            phase->m_V[0][n] = 1.0 / phase->m_order; // phi_func(n+1)
        }
 
        // V evaluation value shuffling row n column n-1
        for (size_t n = 2; n < phase->m_order; ++n)
        {
            phase->m_V[n][n - 1] = 1.0; // constant 1
        }
 
        phase->m_numMultiStepValues         = 1;
        phase->m_numMultiStepImplicitDerivs = 0;
        phase->m_numMultiStepExplicitDerivs = phase->m_order - 1;
        phase->m_timeLevelOffset    = Array<OneD, size_t>(phase->m_numsteps);
        phase->m_timeLevelOffset[0] = 0;
 
        // For order > 1 derivatives are needed.
        for (size_t n = 1; n < phase->m_order; ++n)
        {
            phase->m_timeLevelOffset[n] = n;
        }
 
        phase->CheckAndVerify();
    }

References Nektar::LibUtilities::eExponential.

Referenced by EulerExponentialTimeIntegrationScheme().

◆ v_GetFullName()

LUE std::string Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::v_GetFullName ( ) const

inlineoverrideprotectedvirtual

Reimplemented from Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 213 of file EulerExponentialTimeIntegrationSchemes.h.

    {
        return GetVariant() + GetName() + "Order" + std::to_string(GetOrder());
    }

References Nektar::LibUtilities::TimeIntegrationScheme::GetName(), Nektar::LibUtilities::TimeIntegrationScheme::GetOrder(), and Nektar::LibUtilities::TimeIntegrationScheme::GetVariant().

◆ v_GetName()

LUE std::string Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::v_GetName ( ) const

inlineoverrideprotectedvirtual

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 208 of file EulerExponentialTimeIntegrationSchemes.h.

    {
        return std::string("EulerExponential");
    }

◆ v_GetTimeStability()

LUE NekDouble Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::v_GetTimeStability ( ) const

inlineoverrideprotectedvirtual

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 218 of file EulerExponentialTimeIntegrationSchemes.h.

    {
        return 1.0;
    }

◆ v_InitializeSecondaryData()

void Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::v_InitializeSecondaryData	(	TimeIntegrationAlgorithmGLM *	phase,
		NekDouble	deltaT
	)		const

inlineoverrideprotectedvirtual

Lambda Matrix Assumption, member variable phase->m_L

The one-dimensional Lambda matrix is a diagonal matrix thus values are non zero if and only i=j. As such, the diagonal Lambda values are stored in an array of complex numbers.

Reimplemented from Nektar::LibUtilities::TimeIntegrationSchemeGLM.

Definition at line 223 of file EulerExponentialTimeIntegrationSchemes.h.

    {
        /**
         * \brief Lambda Matrix Assumption, member variable phase->m_L
         *
         * The one-dimensional Lambda matrix is a diagonal
         * matrix thus values are non zero if and only i=j. As such,
         * the diagonal Lambda values are stored in an array of
         * complex numbers.
         */
 
        ASSERTL1(phase->m_nvars == phase->m_L.size(),
                 "The number of variables does not match "
                 "the number of exponential coefficents.");
 
        phase->m_A_phi = Array<OneD, Array<TwoD, NekDouble>>(phase->m_nvars);
        phase->m_B_phi = Array<OneD, Array<TwoD, NekDouble>>(phase->m_nvars);
        phase->m_U_phi = Array<OneD, Array<TwoD, NekDouble>>(phase->m_nvars);
        phase->m_V_phi = Array<OneD, Array<TwoD, NekDouble>>(phase->m_nvars);
 
        Array<OneD, NekDouble> phi = Array<OneD, NekDouble>(phase->m_order);
 
        for (size_t k = 0; k < phase->m_nvars; ++k)
        {
            // B Phi function for first row first column
            if (phase->m_variant == "Lawson")
            {
                phi[0] = phi_function(0, deltaT * phase->m_L[k]).real();
            }
            else if (phase->m_variant == "Norsett")
            {
                phi[0] = phi_function(1, deltaT * phase->m_L[k]).real();
            }
            else
            {
                ASSERTL1(false, "Cannot call EulerExponential directly "
                                "use the variant 'Lawson' or 'Norsett'.");
            }
 
            // Set up for multiple steps. For multiple steps one needs
            // to weight the phi functions in much the same way there
            // are weights for multi-step Adams-Bashfort methods.
 
            // For order N the weights are an N x N matrix with
            // values: W[j][i] = std::pow(i, j) and phi_func = W phi.
            // There are other possible wieghting schemes.
            if (phase->m_order == 1)
            {
                // Nothing to do as the value is set above and the
                // weight is just 1.
            }
            else if (phase->m_order == 2)
            {
                // For Order 2 the weights are simply : 1 1
                //                                      0 1
 
                // If one were to solve the system of equations it
                // simply results in subtracting the second order
                // value from the first order value.
                phi[1] = phi_function(2, deltaT * phase->m_L[k]).real();
 
                phi[0] -= phi[1];
            }
            else if (phase->m_order == 3)
            {
                Array<OneD, NekDouble> phi_func =
                    Array<OneD, NekDouble>(phase->m_order);
 
                phi_func[0] = phi[0];
 
                for (size_t m = 1; m < phase->m_order; ++m)
                {
                    phi_func[m] =
                        phi_function(m + 1, deltaT * phase->m_L[k]).real();
                }
 
                NekDouble W[3][3];
 
                // Set up the wieghts and calculate the determinant.
                for (size_t j = 0; j < phase->m_order; ++j)
                {
                    for (size_t i = 0; i < phase->m_order; ++i)
                    {
                        W[j][i] = std::pow(i, j);
                    }
                }
 
                NekDouble W_det = Determinant<3>(W);
 
                // Solve the series of equations using Cramer's rule.
                for (size_t m = 0; m < phase->m_order; ++m)
                {
                    // Assemble the working matrix for this solution.
                    for (size_t j = 0; j < phase->m_order; ++j)
                    {
                        for (size_t i = 0; i < phase->m_order; ++i)
                        {
                            // Fill in the mth column for the mth
                            // solution using the phi function value
                            // otherwise utilize the weights.
                            W[i][j] = (j == m) ? phi_func[i] : std::pow(j, i);
                        }
                    }
 
                    // Get the mth solutiion.
                    phi[m] = Determinant<3>(W) / W_det;
                }
            }
            else if (phase->m_order == 4)
            {
                Array<OneD, NekDouble> phi_func =
                    Array<OneD, NekDouble>(phase->m_order);
 
                phi_func[0] = phi[0];
 
                for (size_t m = 1; m < phase->m_order; ++m)
                {
                    phi_func[m] =
                        phi_function(m + 1, deltaT * phase->m_L[k]).real();
                }
 
                NekDouble W[4][4];
 
                // Set up the weights and calculate the determinant.
                for (size_t j = 0; j < phase->m_order; ++j)
                {
                    for (size_t i = 0; i < phase->m_order; ++i)
                    {
                        W[j][i] = std::pow(i, j);
                    }
                }
 
                NekDouble W_det = Determinant<4>(W);
 
                // Solve the series of equations using Cramer's rule.
                for (size_t m = 0; m < phase->m_order; ++m)
                {
                    // Assemble the working matrix for this solution.
                    for (size_t j = 0; j < phase->m_order; ++j)
                    {
                        for (size_t i = 0; i < phase->m_order; ++i)
                        {
                            // Fill in the mth column for the mth
                            // solution using the phi function value
                            // otherwise utilize the weights.
                            W[i][j] = (j == m) ? phi_func[i] : std::pow(j, i);
                        }
                    }
 
                    // Get the mth solutiion.
                    phi[m] = Determinant<4>(W) / W_det;
                }
            }
            else
            {
                ASSERTL1(false, "Not set up for more than 4th Order.");
            }
 
            // Create the phi based Butcher tableau matrices. Note
            // these matrices are set up using a general formational
            // based on the number of steps (i.e. the order).
            phase->m_A_phi[k] = Array<TwoD, NekDouble>(phase->m_numstages,
                                                       phase->m_numstages, 0.0);
            phase->m_B_phi[k] = Array<TwoD, NekDouble>(phase->m_numsteps,
                                                       phase->m_numstages, 0.0);
            phase->m_U_phi[k] = Array<TwoD, NekDouble>(phase->m_numstages,
                                                       phase->m_numsteps, 0.0);
            phase->m_V_phi[k] = Array<TwoD, NekDouble>(phase->m_numsteps,
                                                       phase->m_numsteps, 0.0);
 
            // B Phi function for first row first column.
            phase->m_B_phi[k][0][0] = phi[0];
 
            // B evaluation value shuffling second row first column.
            if (phase->m_order > 1)
            {
                phase->m_B_phi[k][1][0] = 1.0; // constant 1
            }
 
            // U Curent time step evaluation first row first column.
            phase->m_U_phi[k][0][0] = 1.0; // constant 1
 
            // V Phi function for first row first column.
            phase->m_V_phi[k][0][0] =
                phi_function(0, deltaT * phase->m_L[k]).real();
 
            // V Phi function for first row additional columns.
            for (size_t n = 1; n < phase->m_order; ++n)
            {
                phase->m_V_phi[k][0][n] = phi[n];
            }
 
            // V evaluation value shuffling row n column n-1.
            for (size_t n = 2; n < phase->m_order; ++n)
            {
                phase->m_V_phi[k][n][n - 1] = 1.0; // constant 1
            }
        }
    }

Member Data Documentation

◆ className

std::string Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::className

static

Definition at line 108 of file EulerExponentialTimeIntegrationSchemes.h.

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Private Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ EulerExponentialTimeIntegrationScheme()

◆ ~EulerExponentialTimeIntegrationScheme()

Member Function Documentation

◆ create()

◆ factorial()

◆ phi_function()

◆ SetExponentialCoefficients()

◆ SetupSchemeData()

◆ v_GetFullName()

◆ v_GetName()

◆ v_GetTimeStability()

◆ v_InitializeSecondaryData()

Member Data Documentation

◆ className