#include <ForcingSyntheticEddy.h>

Inheritance diagram for Nektar::SolverUtils::ForcingSyntheticEddy:

Static Public Member Functions
static SOLVER_UTILS_EXPORT ForcingSharedPtr	create (const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
	Creates an instance of this class. More...

Static Public Member Functions inherited from Nektar::SolverUtils::Forcing
static SOLVER_UTILS_EXPORT std::vector< ForcingSharedPtr >	Load (const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields=0)

Static Public Attributes
static SOLVER_UTILS_EXPORT std::string	className
	Name of the class. More...

Protected Member Functions
SOLVER_UTILS_EXPORT void	v_InitObject (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce) override
	Read input from xml file and initialise the class members. The main parameters are the characteristic lengths, Reynolds stresses and the synthetic eddy region (box of eddies). More...

SOLVER_UTILS_EXPORT void	v_Apply (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time) override
	Apply forcing term if an eddy left the box of eddies and update the eddies positions. More...

SOLVER_UTILS_EXPORT void	v_ApplyCoeff (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time) override
	Apply forcing term if an eddy left the box of eddies and update the eddies positions. More...

SOLVER_UTILS_EXPORT void	SetCholeskyReyStresses (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
	Calculates the Cholesky decomposition of the Reynolds Stresses in each degree of freedom of the mesh. More...

SOLVER_UTILS_EXPORT void	SetNumberOfEddies ()
	Set the Number of the Eddies in the box. More...

SOLVER_UTILS_EXPORT void	ComputeRefLenghts ()
	Set reference lengths. More...

SOLVER_UTILS_EXPORT void	ComputeXiMax ()
	Set Xi max. More...

SOLVER_UTILS_EXPORT Array< OneD, Array< OneD, int > >	GenerateRandomOneOrMinusOne ()
	Compute Random 1 or -1 value. More...

SOLVER_UTILS_EXPORT NekDouble	ComputeConstantC (int row, int col)
	Compute Constant C. More...

SOLVER_UTILS_EXPORT NekDouble	ComputeGaussian (NekDouble coord, NekDouble xiMaxVal, NekDouble constC=1.0)
	Compute Gaussian. More...

SOLVER_UTILS_EXPORT bool	InsideBoxOfEddies (NekDouble coord0, NekDouble coord1, NekDouble coord2)
	Check if point is inside the box of eddies. More...

SOLVER_UTILS_EXPORT Array< OneD, Array< OneD, NekDouble > >	ComputeStochasticSignal (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
	Compute Stochastic Signal. More...

SOLVER_UTILS_EXPORT Array< OneD, Array< OneD, NekDouble > >	ComputeVelocityFluctuation (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, Array< OneD, NekDouble > > stochasticSignal)
	Compute Velocity Fluctuation. More...

SOLVER_UTILS_EXPORT Array< OneD, Array< OneD, NekDouble > >	ComputeCharConvTurbTime (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
	Compute Characteristic Convective Turbulent Time. More...

SOLVER_UTILS_EXPORT Array< OneD, Array< OneD, NekDouble > >	ComputeSmoothingFactor (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, Array< OneD, NekDouble > > convTurb)
	Compute Smoothing Factor. More...

SOLVER_UTILS_EXPORT void	SetBoxOfEddiesMask (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
	Set box of eddies mask. More...

SOLVER_UTILS_EXPORT void	InitialiseForcingEddy (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
	Initialise forcing term for each eddy. More...

SOLVER_UTILS_EXPORT void	ComputeInitialRandomLocationOfEddies ()
	Compute the initial position of the eddies inside the box. More...

SOLVER_UTILS_EXPORT void	UpdateEddiesPositions ()
	Update positions of the eddies inside the box. More...

SOLVER_UTILS_EXPORT void	RemoveEddiesFromForcing (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
	Remove eddy from forcing term. More...

SOLVER_UTILS_EXPORT void	ComputeInitialLocationTestCase ()
	Compute the initial location of the eddies for the test case. More...

SOLVER_UTILS_EXPORT	ForcingSyntheticEddy (const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation)

SOLVER_UTILS_EXPORT	~ForcingSyntheticEddy (void) override=default

Protected Member Functions inherited from Nektar::SolverUtils::Forcing
SOLVER_UTILS_EXPORT	Forcing (const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation)
	Constructor. More...

virtual SOLVER_UTILS_EXPORT void	v_InitObject (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)=0

virtual SOLVER_UTILS_EXPORT void	v_Apply (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time)=0

virtual SOLVER_UTILS_EXPORT void	v_PreApply (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time)

virtual SOLVER_UTILS_EXPORT void	v_ApplyCoeff (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time)

SOLVER_UTILS_EXPORT SessionFunctionSharedPtr	GetFunction (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const LibUtilities::SessionReaderSharedPtr &pSession, std::string pName, bool pCache=false)
	Get a SessionFunction by name. More...

SOLVER_UTILS_EXPORT void	EvaluateTimeFunction (LibUtilities::SessionReaderSharedPtr pSession, std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, NekDouble pTime=NekDouble(0))

SOLVER_UTILS_EXPORT void	EvaluateTimeFunction (const NekDouble pTime, const LibUtilities::EquationSharedPtr &pEqn, Array< OneD, NekDouble > &pArray)

Protected Attributes
std::map< int, LibUtilities::EquationSharedPtr >	m_R

Array< OneD, Array< OneD, NekDouble > >	m_Cholesky
	Cholesky decomposition of the Reynolds Stresses for all points in the mesh. More...

NekDouble	m_Ub
	Bulk velocity. More...

NekDouble	m_sigma
	Standard deviation. More...

Array< OneD, NekDouble >	m_lyz
	Inlet length in the y-direction and z-direction. More...

Array< OneD, NekDouble >	m_rc
	Center of the inlet plane. More...

int	m_N
	Number of eddies in the box. More...

Array< OneD, NekDouble >	m_l
	Characteristic lenght scales. More...

Array< OneD, NekDouble >	m_lref
	Reference lenghts. More...

Array< OneD, NekDouble >	m_xiMax
	Xi max. More...

NekDouble	m_xiMaxMin = 0.4

int	m_spacedim
	Space dimension. More...

NekDouble	m_gamma
	Ration of specific heats. More...

Array< OneD, int >	m_mask
	Box of eddies mask. More...

Array< OneD, Array< OneD, NekDouble > >	m_eddyPos
	Eddy position. More...

bool	m_calcForcing {true}
	Check when the forcing should be applied. More...

std::vector< unsigned int >	m_eddiesIDForcing
	Eddies that add to the domain. More...

NekDouble	m_currTime
	Current time. More...

bool	m_tCase
	Check for test case. More...

Array< OneD, Array< OneD, Array< OneD, NekDouble > > >	m_ForcingEddy
	Forcing for each eddy. More...

Protected Attributes inherited from Nektar::SolverUtils::Forcing
LibUtilities::SessionReaderSharedPtr	m_session
	Session reader. More...

const std::weak_ptr< EquationSystem >	m_equ
	Weak pointer to equation system using this forcing. More...

Array< OneD, Array< OneD, NekDouble > >	m_Forcing
	Evaluated forcing function. More...

int	m_NumVariable
	Number of variables. More...

std::map< std::string, SolverUtils::SessionFunctionSharedPtr >	m_sessionFunctions
	Map of known SessionFunctions. More...

Friends
class	MemoryManager< ForcingSyntheticEddy >

Additional Inherited Members
Public Member Functions inherited from Nektar::SolverUtils::Forcing
virtual SOLVER_UTILS_EXPORT	~Forcing ()

SOLVER_UTILS_EXPORT void	InitObject (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
	Initialise the forcing object. More...

SOLVER_UTILS_EXPORT void	Apply (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time)
	Apply the forcing. More...

SOLVER_UTILS_EXPORT void	PreApply (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time)
	Change the advection velocity before applying the forcing. For example, subtracting the frame velocity from the advection velocity in the MovingRefercenceFrame. More...

SOLVER_UTILS_EXPORT void	ApplyCoeff (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time)
	Apply the forcing. More...

SOLVER_UTILS_EXPORT const Array< OneD, const Array< OneD, NekDouble > > &	GetForces ()

SOLVER_UTILS_EXPORT Array< OneD, Array< OneD, NekDouble > > &	UpdateForces ()

Detailed Description

Definition at line 45 of file ForcingSyntheticEddy.h.

Constructor & Destructor Documentation

◆ ForcingSyntheticEddy()

Nektar::SolverUtils::ForcingSyntheticEddy::ForcingSyntheticEddy	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const std::weak_ptr< EquationSystem > &	pEquation
	)

protected

Definition at line 49 of file ForcingSyntheticEddy.cpp.

    : Forcing(pSession, pEquation)
{
}

◆ ~ForcingSyntheticEddy()

SOLVER_UTILS_EXPORT Nektar::SolverUtils::ForcingSyntheticEddy::~ForcingSyntheticEddy ( void )

overrideprotecteddefault

Member Function Documentation

◆ ComputeCharConvTurbTime()

Array< OneD, Array< OneD, NekDouble > > Nektar::SolverUtils::ForcingSyntheticEddy::ComputeCharConvTurbTime ( const Array< OneD, MultiRegions::ExpListSharedPtr > & pFields )

protected

Compute Characteristic Convective Turbulent Time.

Compute characteristic convective turbulent time.

Parameters

pFields Pointer to fields.

Definition at line 285 of file ForcingSyntheticEddy.cpp.

{
    // Total number of quadrature points
    int nqTot = pFields[0]->GetTotPoints();
    // Characteristic convective turbulent time
    Array<OneD, Array<OneD, NekDouble>> convTurbTime(m_spacedim);
 
    for (size_t k = 0; k < m_spacedim; ++k)
    {
        convTurbTime[k] = Array<OneD, NekDouble>(nqTot, 0.0);
        for (size_t i = 0; i < nqTot; ++i)
        {
            NekDouble convTurbLength = m_xiMaxMin * m_lref[0];
            // 3*k because of the structure of the m_l parameter
            // to obtain lxk.
            if ((m_l[3 * k] > m_xiMaxMin * m_lref[0]) && (m_mask[i]))
            {
                convTurbLength = m_l[3 * k];
            }
            convTurbTime[k][i] = convTurbLength / m_Ub;
        }
    }
 
    return convTurbTime;
}

References m_l, m_lref, m_mask, m_spacedim, m_Ub, and m_xiMaxMin.

Referenced by Nektar::SolverUtils::ForcingCFSSyntheticEddy::CalculateForcing(), and Nektar::SolverUtils::ForcingIncNSSyntheticEddy::CalculateForcing().

◆ ComputeConstantC()

NekDouble Nektar::SolverUtils::ForcingSyntheticEddy::ComputeConstantC	(	int	row,
		int	col
	)

protected

Compute Constant C.

Compute constant C for the gaussian funcion.

Parameters

row	index for the rows of the matrix.
col	index for the columns of the matrix.

Returns: Value of C.

Definition at line 614 of file ForcingSyntheticEddy.cpp.

{
    NekDouble sizeLenScale = m_xiMax[col * m_spacedim + row];
 
    // Integration
    NekDouble sum  = 0.0;
    NekDouble step = 0.025;
    NekDouble xi0  = -1;
    NekDouble xif  = 1;
    while (xi0 < xif)
    {
        sum += (ComputeGaussian(xi0 + step, sizeLenScale) *
                    ComputeGaussian(xi0 + step, sizeLenScale) +
                ComputeGaussian(xi0, sizeLenScale) *
                    ComputeGaussian(xi0, sizeLenScale));
        xi0 += step;
    }
 
    return (0.5 * 0.5 * step * sum);
}

References ComputeGaussian(), m_spacedim, and m_xiMax.

Referenced by ComputeStochasticSignal().

◆ ComputeGaussian()

NekDouble Nektar::SolverUtils::ForcingSyntheticEddy::ComputeGaussian	(	NekDouble	coord,
		NekDouble	xiMaxVal,
		NekDouble	constC = `1.0`
	)

protected

Compute Gaussian.

Compute standard Gaussian with zero mean.

Parameters

coord Coordianate.

Returns: Gaussian value for the coordinate.

Definition at line 591 of file ForcingSyntheticEddy.cpp.

{
    NekDouble epsilon = 1e-6;
    if (abs(coord) <= (xiMaxVal + epsilon))
    {
        return ((1.0 / (m_sigma * sqrt(2.0 * M_PI * constC))) *
                exp(-0.5 * (coord / m_sigma) * (coord / m_sigma)));
    }
    else
    {
        return 0.0;
    }
}

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

Referenced by ComputeConstantC(), and ComputeStochasticSignal().

◆ ComputeInitialLocationTestCase()

void Nektar::SolverUtils::ForcingSyntheticEddy::ComputeInitialLocationTestCase ( )

protected

Compute the initial location of the eddies for the test case.

Place eddies in specific locations in the test case geometry for consistency and comparison.

This function was design for a three-dimensional channel flow test case (ChanFlow3d_infTurb). It is only called for testing purposes (unit test).

Definition at line 900 of file ForcingSyntheticEddy.cpp.

{
    m_N       = 3; // Redefine number of eddies
    m_eddyPos = Array<OneD, Array<OneD, NekDouble>>(m_N);
 
    // First eddy (center)
    m_eddyPos[0]    = Array<OneD, NekDouble>(m_spacedim);
    m_eddyPos[0][0] = (m_rc[0] - m_lref[0]) + 0.6 * 2 * m_lref[0];
    m_eddyPos[0][1] = m_rc[1];
    m_eddyPos[0][2] = m_rc[2];
 
    // Second eddy (top)
    m_eddyPos[1]    = Array<OneD, NekDouble>(m_spacedim);
    m_eddyPos[1][0] = (m_rc[0] - m_lref[0]) + 0.7 * 2 * m_lref[0];
    m_eddyPos[1][1] = (m_rc[1] - 0.5 * m_lyz[0]) + 0.9 * m_lyz[0];
    m_eddyPos[1][2] = m_rc[2];
 
    // Third eddy (bottom)
    m_eddyPos[2]    = Array<OneD, NekDouble>(m_spacedim);
    m_eddyPos[2][0] = (m_rc[0] - m_lref[0]) + 0.8 * 2 * m_lref[0];
    m_eddyPos[2][1] = (m_rc[1] - 0.5 * m_lyz[0]) + 0.1 * m_lyz[0];
    m_eddyPos[2][2] = m_rc[2];
}

References m_eddyPos, m_lref, m_lyz, m_N, m_rc, and m_spacedim.

Referenced by v_InitObject().

◆ ComputeInitialRandomLocationOfEddies()

void Nektar::SolverUtils::ForcingSyntheticEddy::ComputeInitialRandomLocationOfEddies ( )

protected

Compute the initial position of the eddies inside the box.

Calculate distribution of eddies in the box.

Definition at line 565 of file ForcingSyntheticEddy.cpp.

{
    m_eddyPos = Array<OneD, Array<OneD, NekDouble>>(m_N);
 
    for (size_t n = 0; n < m_N; ++n)
    {
        m_eddyPos[n] = Array<OneD, NekDouble>(m_spacedim);
        // Generate randomly eddies inside the box
        m_eddyPos[n][0] =
            (m_rc[0] - m_lref[0]) +
            (NekSingle(std::rand()) / NekSingle(RAND_MAX)) * 2 * m_lref[0];
        m_eddyPos[n][1] =
            (m_rc[1] - 0.5 * m_lyz[0]) +
            (NekSingle(std::rand()) / NekSingle(RAND_MAX)) * m_lyz[0];
        m_eddyPos[n][2] =
            (m_rc[2] - 0.5 * m_lyz[1]) +
            (NekSingle(std::rand()) / NekSingle(RAND_MAX)) * m_lyz[1];
    }
}

References m_eddyPos, m_lref, m_lyz, m_N, m_rc, and m_spacedim.

Referenced by v_InitObject().

◆ ComputeRefLenghts()

void Nektar::SolverUtils::ForcingSyntheticEddy::ComputeRefLenghts ( )

protected

Set reference lengths.

Calculates the reference lenghts.

Definition at line 768 of file ForcingSyntheticEddy.cpp.

{
    m_lref    = Array<OneD, NekDouble>(m_spacedim, 0.0);
    m_lref[0] = m_l[0];
    m_lref[1] = m_l[1];
    m_lref[2] = m_l[2];
 
    // The l_{x}^{ref}, l_{y}^{ref} and l_{z}^{ref}
    // are the maximum among the velocity components
    // over the domain.
    for (size_t i = 0; i < m_spacedim; ++i)
    {
        for (size_t j = 1; j < m_spacedim; ++j)
        {
            if (m_l[m_spacedim * j + i] > m_lref[i])
            {
                m_lref[i] = m_l[m_spacedim * j + i];
            }
        }
    }
}

References m_l, m_lref, and m_spacedim.

Referenced by v_InitObject().

◆ ComputeSmoothingFactor()

Array< OneD, Array< OneD, NekDouble > > Nektar::SolverUtils::ForcingSyntheticEddy::ComputeSmoothingFactor	(	const Array< OneD, MultiRegions::ExpListSharedPtr > &	pFields,
		Array< OneD, Array< OneD, NekDouble > >	convTurbTime
	)

protected

Compute Smoothing Factor.

Compute smoothing factor to avoid strong variations of the source term across the domain.

Parameters

pFields	Pointer to fields.
convTurbTime	Convective turbulent time.

Definition at line 320 of file ForcingSyntheticEddy.cpp.

{
    // Total number of quadrature points
    int nqTot = pFields[0]->GetTotPoints();
    // Number of elements
    int nElmts = pFields[0]->GetNumElmts();
    // Total number of quadrature points of each element
    int nqe;
    // Characteristic convective turbulent time
    Array<OneD, Array<OneD, NekDouble>> smoothFac(m_spacedim);
    // Counter
    int count = 0;
    // module
    NekDouble mod;
    // Create Array with zeros
    for (size_t i = 0; i < m_spacedim; ++i)
    {
        smoothFac[i] = Array<OneD, NekDouble>(nqTot, 0.0);
    }
 
    for (size_t e = 0; e < nElmts; ++e)
    {
        nqe = pFields[0]->GetExp(e)->GetTotPoints();
 
        Array<OneD, NekDouble> coords0(nqe, 0.0);
        Array<OneD, NekDouble> coords1(nqe, 0.0);
        Array<OneD, NekDouble> coords2(nqe, 0.0);
 
        pFields[0]->GetExp(e)->GetCoords(coords0, coords1, coords2);
 
        for (size_t i = 0; i < nqe; ++i)
        {
            if (m_mask[count + i])
            {
                mod = (coords0[i] - m_rc[0] + m_lref[0]) *
                      (coords0[i] - m_rc[0] + m_lref[0]);
 
                smoothFac[0][count + i] =
                    exp((-0.5 * M_PI * mod) /
                        (convTurbTime[0][count + i] *
                         convTurbTime[0][count + i] * m_Ub * m_Ub));
                smoothFac[1][count + i] =
                    exp((-0.5 * M_PI * mod) /
                        (convTurbTime[1][count + i] *
                         convTurbTime[1][count + i] * m_Ub * m_Ub));
                smoothFac[2][count + i] =
                    exp((-0.5 * M_PI * mod) /
                        (convTurbTime[2][count + i] *
                         convTurbTime[2][count + i] * m_Ub * m_Ub));
            }
        }
        count += nqe;
    }
 
    return smoothFac;
}

References m_lref, m_mask, m_rc, m_spacedim, and m_Ub.

Referenced by Nektar::SolverUtils::ForcingCFSSyntheticEddy::CalculateForcing(), and Nektar::SolverUtils::ForcingIncNSSyntheticEddy::CalculateForcing().

◆ ComputeStochasticSignal()

Array< OneD, Array< OneD, NekDouble > > Nektar::SolverUtils::ForcingSyntheticEddy::ComputeStochasticSignal ( const Array< OneD, MultiRegions::ExpListSharedPtr > & pFields )

protected

Compute Stochastic Signal.

Compute stochastic signal.

Parameters

pFields Pointer to fields.

Returns: Stochastic signal.

Definition at line 429 of file ForcingSyntheticEddy.cpp.

{
    // Total number of quadrature points
    int nqTot = pFields[0]->GetTotPoints();
    // Number of elements
    int nElmts = pFields[0]->GetNumElmts();
    // Total number of quadrature points of each element
    int nqe;
    // Stochastic Signal vector
    Array<OneD, Array<OneD, NekDouble>> stochasticSignal(m_N);
    // Random numbers: -1 and 1
    Array<OneD, Array<OneD, int>> epsilonSign;
 
    // Generate only for the new eddies after the first time step
    epsilonSign = GenerateRandomOneOrMinusOne();
 
    // Calculate the stochastic signal for the eddies
    for (auto &n : m_eddiesIDForcing)
    {
        stochasticSignal[n] = Array<OneD, NekDouble>(nqTot * m_spacedim, 0.0);
 
        // Evaluate the function at interpolation points for each component
        for (size_t j = 0; j < m_spacedim; ++j)
        {
            // Count the number of quadrature points
            int nqeCount = 0;
 
            for (size_t e = 0; e < nElmts; ++e)
            {
                nqe = pFields[0]->GetExp(e)->GetTotPoints();
 
                Array<OneD, NekDouble> coords0(nqe, 0.0);
                Array<OneD, NekDouble> coords1(nqe, 0.0);
                Array<OneD, NekDouble> coords2(nqe, 0.0);
 
                pFields[0]->GetExp(e)->GetCoords(coords0, coords1, coords2);
 
                // i: degrees of freedom, j: direction, n: eddies
                for (size_t i = 0; i < nqe; ++i)
                {
                    if (m_mask[nqeCount + i])
                    {
                        stochasticSignal[n][j * nqTot + nqeCount + i] =
                            epsilonSign[j][n] *
                            ComputeGaussian((coords0[i] - m_eddyPos[n][0]) /
                                                m_lref[0],
                                            m_xiMax[j * m_spacedim + 0],
                                            ComputeConstantC(0, j)) *
                            ComputeGaussian((coords1[i] - m_eddyPos[n][1]) /
                                                m_lref[1],
                                            m_xiMax[j * m_spacedim + 1],
                                            ComputeConstantC(1, j)) *
                            ComputeGaussian((coords2[i] - m_eddyPos[n][2]) /
                                                m_lref[2],
                                            m_xiMax[j * m_spacedim + 2],
                                            ComputeConstantC(2, j));
                    }
                }
                nqeCount += nqe;
            }
        }
    }
 
    return stochasticSignal;
}

References ComputeConstantC(), ComputeGaussian(), GenerateRandomOneOrMinusOne(), m_eddiesIDForcing, m_eddyPos, m_lref, m_mask, m_N, m_spacedim, and m_xiMax.

Referenced by Nektar::SolverUtils::ForcingCFSSyntheticEddy::CalculateForcing(), and Nektar::SolverUtils::ForcingIncNSSyntheticEddy::CalculateForcing().

◆ ComputeVelocityFluctuation()

Array< OneD, Array< OneD, NekDouble > > Nektar::SolverUtils::ForcingSyntheticEddy::ComputeVelocityFluctuation	(	const Array< OneD, MultiRegions::ExpListSharedPtr > &	pFields,
		Array< OneD, Array< OneD, NekDouble > >	stochasticSignal
	)

protected

Compute Velocity Fluctuation.

Calculate velocity fluctuation for the source term.

Parameters

pFields	Pointer to fields.
stochasticSignal	Stochastic signal.

Returns: Velocity fluctuation.

Definition at line 387 of file ForcingSyntheticEddy.cpp.

{
    // Total number of quadrature points
    int nqTot = pFields[0]->GetTotPoints();
    // Velocity fluctuation vector
    Array<OneD, Array<OneD, NekDouble>> velFluc(m_N);
    // Control loop for the m_Cholesky (Cholesky decomposition matrix)
    int l;
 
    for (auto &n : m_eddiesIDForcing)
    {
        velFluc[n] = Array<OneD, NekDouble>(nqTot * m_spacedim, 0.0);
        for (size_t k = 0; k < m_spacedim; ++k)
        {
            for (size_t j = 0; j < k + 1; ++j)
            {
                for (size_t i = 0; i < nqTot; ++i)
                {
                    if (m_mask[i])
                    {
                        l = k + j * (2 * m_spacedim - j - 1) * 0.5;
                        velFluc[n][k * nqTot + i] +=
                            m_Cholesky[i][l] *
                            stochasticSignal[n][j * nqTot + i];
                    }
                }
            }
        }
    }
 
    return velFluc;
}

References m_Cholesky, m_eddiesIDForcing, m_mask, m_N, and m_spacedim.

Referenced by Nektar::SolverUtils::ForcingCFSSyntheticEddy::CalculateForcing(), and Nektar::SolverUtils::ForcingIncNSSyntheticEddy::CalculateForcing().

◆ ComputeXiMax()

void Nektar::SolverUtils::ForcingSyntheticEddy::ComputeXiMax ( )

protected

Set Xi max.

Calculates the \(\xi_{max}\).

Definition at line 793 of file ForcingSyntheticEddy.cpp.

{
    NekDouble value;
    m_xiMax = Array<OneD, NekDouble>(m_spacedim * m_spacedim, 0.0);
 
    for (size_t i = 0; i < m_spacedim; i++)
    {
        for (size_t j = 0; j < m_spacedim; j++)
        {
            value = (m_l[m_spacedim * j + i] / m_lref[i]);
            if (value > m_xiMaxMin)
            {
                m_xiMax[m_spacedim * j + i] = value;
            }
            else
            {
                m_xiMax[m_spacedim * j + i] = m_xiMaxMin;
            }
        }
    }
}

References m_l, m_lref, m_spacedim, m_xiMax, and m_xiMaxMin.

Referenced by v_InitObject().

◆ create()

static SOLVER_UTILS_EXPORT ForcingSharedPtr Nektar::SolverUtils::ForcingSyntheticEddy::create	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const std::weak_ptr< EquationSystem > &	pEquation,
		const Array< OneD, MultiRegions::ExpListSharedPtr > &	pFields,
		const unsigned int &	pNumForcingFields,
		const TiXmlElement *	pForce
	)

inlinestatic

Creates an instance of this class.

Definition at line 51 of file ForcingSyntheticEddy.h.

    {
        ForcingSharedPtr p =
            MemoryManager<ForcingSyntheticEddy>::AllocateSharedPtr(pSession,
                                                                   pEquation);
        p->InitObject(pFields, pNumForcingFields, pForce);
        return p;
    }

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

◆ GenerateRandomOneOrMinusOne()

Array< OneD, Array< OneD, int > > Nektar::SolverUtils::ForcingSyntheticEddy::GenerateRandomOneOrMinusOne ( )

protected

Compute Random 1 or -1 value.

Generate random 1 or -1 values to be use to compute the stochastic signal.

Returns: ramdom 1 or -1 values.

Definition at line 640 of file ForcingSyntheticEddy.cpp.

{
    Array<OneD, Array<OneD, int>> epsilonSign(m_spacedim);
 
    // j: component of the stochastic signal
    // n: eddy
    for (size_t j = 0; j < m_spacedim; ++j)
    {
        epsilonSign[j] = Array<OneD, int>(m_N, 0.0);
        for (auto &n : m_eddiesIDForcing)
        {
            // Convert to -1 or to 1
            // Check if test case
            if (!m_tCase)
            {
                epsilonSign[j][n] =
                    ((NekSingle(std::rand()) / NekSingle(RAND_MAX)) <= 0.5) ? -1
                                                                            : 1;
            }
            else
            {
                // Positive values only for the test case
                epsilonSign[j][n] = 1;
            }
        }
    }
 
    return epsilonSign;
}

References m_eddiesIDForcing, m_N, m_spacedim, and m_tCase.

Referenced by ComputeStochasticSignal().

◆ InitialiseForcingEddy()

void Nektar::SolverUtils::ForcingSyntheticEddy::InitialiseForcingEddy ( const Array< OneD, MultiRegions::ExpListSharedPtr > & pFields )

protected

Initialise forcing term for each eddy.

Initialise Forcing term for each eddy.

Definition at line 717 of file ForcingSyntheticEddy.cpp.

{
    // Total number of quadrature points
    int nqTot = pFields[0]->GetTotPoints();
    // Number of Variables
    int nVars     = pFields.size();
    m_ForcingEddy = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(m_N);
 
    for (size_t i = 0; i < m_N; ++i)
    {
        m_ForcingEddy[i] = Array<OneD, Array<OneD, NekDouble>>(nVars);
        for (size_t j = 0; j < nVars; ++j)
        {
            m_ForcingEddy[i][j] = Array<OneD, NekDouble>(nqTot);
            for (size_t k = 0; k < nqTot; ++k)
            {
                m_ForcingEddy[i][j][k] = 0.0;
            }
        }
    }
}

References m_ForcingEddy, and m_N.

Referenced by v_InitObject().

◆ InsideBoxOfEddies()

bool Nektar::SolverUtils::ForcingSyntheticEddy::InsideBoxOfEddies	(	NekDouble	coord0,
		NekDouble	coord1,
		NekDouble	coord2
	)

protected

Check if point is inside the box of eddies.

Check it point is inside the box of eddies.

Parameters

coord0	coordinate in the x-direction.
coord1	coordinate in the y-direction.
coord2	coordinate in the z-direction.

Returns: true or false

Definition at line 748 of file ForcingSyntheticEddy.cpp.

{
    NekDouble tol = 1e-6;
    if ((coord0 >= (m_rc[0] - m_lref[0] - m_lref[0])) &&
        (coord0 <= (m_rc[0] + m_lref[0] + tol)) &&
        (coord1 >= (m_rc[1] - 0.5 * m_lyz[0] - tol)) &&
        (coord1 <= (m_rc[1] + 0.5 * m_lyz[0] + tol)) &&
        (coord2 >= (m_rc[2] - 0.5 * m_lyz[1] - tol)) &&
        (coord2 <= (m_rc[2] + 0.5 * m_lyz[1] + tol)))
    {
        return true;
    }
 
    return false;
}

References m_lref, m_lyz, and m_rc.

Referenced by SetBoxOfEddiesMask().

◆ RemoveEddiesFromForcing()

void Nektar::SolverUtils::ForcingSyntheticEddy::RemoveEddiesFromForcing ( const Array< OneD, MultiRegions::ExpListSharedPtr > & pFields )

protected

Remove eddy from forcing term.

Parameters

pFields Pointer to fields.

Definition at line 542 of file ForcingSyntheticEddy.cpp.

{
    // Total number of quadrature points
    int nqTot = pFields[0]->GetTotPoints();
    // Number of Variables
    int nVars = pFields.size();
 
    for (auto &n : m_eddiesIDForcing)
    {
        for (size_t j = 0; j < nVars; ++j)
        {
            for (size_t i = 0; i < nqTot; ++i)
            {
                m_Forcing[j][i] -= m_ForcingEddy[n][j][i];
            }
        }
    }
}

References m_eddiesIDForcing, Nektar::SolverUtils::Forcing::m_Forcing, and m_ForcingEddy.

Referenced by Nektar::SolverUtils::ForcingCFSSyntheticEddy::CalculateForcing(), and Nektar::SolverUtils::ForcingIncNSSyntheticEddy::CalculateForcing().

◆ SetBoxOfEddiesMask()

void Nektar::SolverUtils::ForcingSyntheticEddy::SetBoxOfEddiesMask ( const Array< OneD, MultiRegions::ExpListSharedPtr > & pFields )

protected

Set box of eddies mask.

Set box of eddies mask to be used to seprate the degrees of freedom (coordinates) inside and outside the box of eddies.

Parameters

pFields Pointer to fields.

Definition at line 678 of file ForcingSyntheticEddy.cpp.

{
    // Total number of quadrature points
    int nqTot = pFields[0]->GetTotPoints();
    // Number of elements
    int nElmts = pFields[0]->GetNumElmts();
    // Total number of quadrature points of each element
    int nqe;
    // Mask
    m_mask = Array<OneD, int>(nqTot, 0); // 0 for ouside, 1 for inside
    // Counter
    int count = 0;
 
    for (size_t e = 0; e < nElmts; ++e)
    {
        nqe = pFields[0]->GetExp(e)->GetTotPoints();
 
        Array<OneD, NekDouble> coords0(nqe, 0.0);
        Array<OneD, NekDouble> coords1(nqe, 0.0);
        Array<OneD, NekDouble> coords2(nqe, 0.0);
 
        pFields[0]->GetExp(e)->GetCoords(coords0, coords1, coords2);
 
        for (size_t i = 0; i < nqe; ++i)
        {
            if (InsideBoxOfEddies(coords0[i], coords1[i], coords2[i]))
            {
                // 0 for outside, 1 for inside
                m_mask[count + i] = 1;
            }
        }
        count += nqe;
    }
}

References InsideBoxOfEddies(), and m_mask.

Referenced by v_InitObject().

◆ SetCholeskyReyStresses()

void Nektar::SolverUtils::ForcingSyntheticEddy::SetCholeskyReyStresses ( const Array< OneD, MultiRegions::ExpListSharedPtr > & pFields )

protected

Calculates the Cholesky decomposition of the Reynolds Stresses in each degree of freedom of the mesh.

Parameters

pFields Pointer to fields.

Definition at line 821 of file ForcingSyntheticEddy.cpp.

{
    // Total number of quadrature points
    int nqTot = pFields[0]->GetTotPoints();
    // Total number of quadrature points of each element
    int nqe;
    // Number of elements
    int nElmts = pFields[0]->GetNumElmts();
    // Count the degrees of freedom
    int nqeCount = 0;
    // Block diagonal size
    int diagSize = m_spacedim * (m_spacedim + 1) * 0.5;
    // Cholesky decomposition matrix for the synthetic eddy region (box)
    m_Cholesky = Array<OneD, Array<OneD, NekDouble>>(nqTot);
 
    for (size_t e = 0; e < nElmts; ++e)
    {
        nqe = pFields[0]->GetExp(e)->GetTotPoints();
 
        Array<OneD, NekDouble> coords0(nqe, 0.0);
        Array<OneD, NekDouble> coords1(nqe, 0.0);
        Array<OneD, NekDouble> coords2(nqe, 0.0);
 
        // Coordinates (for each degree of freedom) for the element k.
        pFields[0]->GetExp(e)->GetCoords(coords0, coords1, coords2);
 
        // Evaluate Cholesky decomposition
        for (size_t i = 0; i < nqe; ++i)
        {
            int l = 0;
            // Size of Cholesky decomposition matrix - aux vector
            Array<OneD, NekDouble> A(diagSize, 0.0);
 
            while (l < diagSize)
            {
                // Variable to compute the Cholesky decomposition for each
                // degree of freedom
                A[l] = m_R[l]->Evaluate(coords0[i], coords1[i], coords2[i]);
                l++;
            }
            int info = 0;
            Lapack::Dpptrf('L', m_spacedim, A.get(), info);
            if (info < 0)
            {
                std::string message =
                    "ERROR: The " + std::to_string(-info) +
                    "th parameter had an illegal parameter for dpptrf";
                NEKERROR(ErrorUtil::efatal, message.c_str());
            }
            m_Cholesky[nqeCount + i] = Array<OneD, NekDouble>(diagSize);
            for (size_t l = 0; l < diagSize; ++l)
            {
                m_Cholesky[nqeCount + i][l] = A[l];
            }
        }
        nqeCount += nqe;
    }
}

References Lapack::Dpptrf(), Nektar::ErrorUtil::efatal, m_Cholesky, m_R, m_spacedim, and NEKERROR.

Referenced by v_InitObject().

◆ SetNumberOfEddies()

void Nektar::SolverUtils::ForcingSyntheticEddy::SetNumberOfEddies ( )

protected

Set the Number of the Eddies in the box.

Calculate the number of eddies that are going to be injected in the synthetic eddy region (box).

Definition at line 885 of file ForcingSyntheticEddy.cpp.

{
    m_N = int((m_lyz[0] * m_lyz[1]) /
              (4 * m_lref[m_spacedim - 2] * m_lref[m_spacedim - 1])) +
          1;
}

References m_lref, m_lyz, m_N, and m_spacedim.

Referenced by v_InitObject().

◆ UpdateEddiesPositions()

void Nektar::SolverUtils::ForcingSyntheticEddy::UpdateEddiesPositions ( )

protected

Update positions of the eddies inside the box.

Update the position of the eddies for every time step.

Definition at line 500 of file ForcingSyntheticEddy.cpp.

{
    NekDouble dt = m_session->GetParameter("TimeStep");
 
    for (size_t n = 0; n < m_N; ++n)
    {
        // Check if any eddy went through the outlet plane (box)
        if ((m_eddyPos[n][0] + m_Ub * dt) < (m_rc[0] + m_lref[0]))
        {
            m_eddyPos[n][0] = m_eddyPos[n][0] + m_Ub * dt;
        }
        else
        {
            // Generate a new one in the inlet plane
            m_eddyPos[n][0] = m_rc[0] - m_lref[0];
            // Check if test case
            if (!m_tCase)
            {
                m_eddyPos[n][1] =
                    (m_rc[1] - 0.5 * m_lyz[0]) +
                    (NekSingle(std::rand()) / NekSingle(RAND_MAX)) * (m_lyz[0]);
                m_eddyPos[n][2] =
                    (m_rc[2] - 0.5 * m_lyz[1] + 0.5 * m_lref[2]) +
                    (NekSingle(std::rand()) / NekSingle(RAND_MAX)) * (m_lyz[1]);
            }
            else
            {
                // same place (center) for the test case
                m_eddyPos[n][1] = m_rc[1];
                m_eddyPos[n][2] = m_rc[2];
            }
            m_eddiesIDForcing.push_back(n);
            m_calcForcing = true;
        }
    }
}

References m_calcForcing, m_eddiesIDForcing, m_eddyPos, m_lref, m_lyz, m_N, m_rc, Nektar::SolverUtils::Forcing::m_session, m_tCase, and m_Ub.

Referenced by Nektar::SolverUtils::ForcingCFSSyntheticEddy::v_Apply(), Nektar::SolverUtils::ForcingIncNSSyntheticEddy::v_Apply(), and Nektar::SolverUtils::ForcingCFSSyntheticEddy::v_ApplyCoeff().

◆ v_Apply()

void Nektar::SolverUtils::ForcingSyntheticEddy::v_Apply	(	const Array< OneD, MultiRegions::ExpListSharedPtr > &	pFields,
		const Array< OneD, Array< OneD, NekDouble > > &	inarray,
		Array< OneD, Array< OneD, NekDouble > > &	outarray,
		const NekDouble &	time
	)

overrideprotectedvirtual

Apply forcing term if an eddy left the box of eddies and update the eddies positions.

Parameters

pFields	Pointer to fields.
inarray	Given fields. The fields are in in physical space.
outarray	Calculated solution after forcing term being applied in physical space.
time	time.

Implements Nektar::SolverUtils::Forcing.

Reimplemented in Nektar::SolverUtils::ForcingCFSSyntheticEddy, and Nektar::SolverUtils::ForcingIncNSSyntheticEddy.

Definition at line 254 of file ForcingSyntheticEddy.cpp.

259{

260}

◆ v_ApplyCoeff()

void Nektar::SolverUtils::ForcingSyntheticEddy::v_ApplyCoeff	(	const Array< OneD, MultiRegions::ExpListSharedPtr > &	pFields,
		const Array< OneD, Array< OneD, NekDouble > > &	inarray,
		Array< OneD, Array< OneD, NekDouble > > &	outarray,
		const NekDouble &	time
	)

overrideprotectedvirtual

Apply forcing term if an eddy left the box of eddies and update the eddies positions.

Parameters

pFields	Pointer to fields.
inarray	Given fields. The fields are in in physical space.
outarray	Calculated solution after forcing term being applied in coefficient space.
time	time.

Reimplemented from Nektar::SolverUtils::Forcing.

Reimplemented in Nektar::SolverUtils::ForcingCFSSyntheticEddy.

Definition at line 272 of file ForcingSyntheticEddy.cpp.

277{

278}

◆ v_InitObject()

void Nektar::SolverUtils::ForcingSyntheticEddy::v_InitObject	(	const Array< OneD, MultiRegions::ExpListSharedPtr > &	pFields,
		const unsigned int &	pNumForcingFields,
		const TiXmlElement *	pForce
	)

overrideprotectedvirtual

Read input from xml file and initialise the class members. The main parameters are the characteristic lengths, Reynolds stresses and the synthetic eddy region (box of eddies).

Parameters

pFields	Pointer to fields.
pNumForcingField	Number of forcing fields.
pForce	Xml element describing the mapping.

Implements Nektar::SolverUtils::Forcing.

Definition at line 65 of file ForcingSyntheticEddy.cpp.

{
    // Space dimension
    m_spacedim = pFields[0]->GetGraph()->GetMeshDimension();
 
    if (m_spacedim != 3)
    {
        NEKERROR(Nektar::ErrorUtil::efatal,
                 "Sythetic eddy method "
                 "only supports fully three-dimensional simulations");
    }
 
    // Get gamma parameter
    m_session->LoadParameter("Gamma", m_gamma, 1.4);
 
    const TiXmlElement *elmtInfTurb;
 
    // Reynolds Stresses
    elmtInfTurb = pForce->FirstChildElement("ReynoldsStresses");
    ASSERTL0(elmtInfTurb, "Requires ReynoldsStresses tag specifying function "
                          "name which prescribes the reynodls stresses.");
 
    std::string reyStressesFuncName = elmtInfTurb->GetText();
    ASSERTL0(m_session->DefinesFunction(reyStressesFuncName),
             "Function '" + reyStressesFuncName +
                 "' is not defined in the session.");
 
    // Reynolds stresses variables. Do not change the order of the variables.
    std::vector<std::string> reyStressesVars = {"r00", "r10", "r20",
                                                "r11", "r21", "r22"};
 
    for (size_t i = 0; i < reyStressesVars.size(); ++i)
    {
        std::string varStr = reyStressesVars[i];
        if (m_session->DefinesFunction(reyStressesFuncName, varStr))
        {
            m_R[i] = m_session->GetFunction(reyStressesFuncName, varStr);
        }
        else
        {
            NEKERROR(Nektar::ErrorUtil::efatal,
                     "Missing parameter '" + varStr +
                         "' in the FUNCTION NAME = " + reyStressesFuncName +
                         ".");
        }
    }
 
    // Characteristic length scales
    elmtInfTurb = pForce->FirstChildElement("CharLengthScales");
    ASSERTL0(elmtInfTurb, "Requires CharLengthScales tag specifying function "
                          "name which prescribes the characteristic "
                          "length scales.");
 
    std::string charLenScalesFuncName = elmtInfTurb->GetText();
    ASSERTL0(m_session->DefinesFunction(charLenScalesFuncName),
             "Function '" + charLenScalesFuncName +
                 "' is not defined in the session.");
 
    // Characteristic length scale variables
    // Do not change the order of the variables
    std::vector<std::string> lenScalesVars = {"l00", "l10", "l20", "l01", "l11",
                                              "l21", "l02", "l12", "l22"};
    LibUtilities::EquationSharedPtr clsAux;
    m_l = Array<OneD, NekDouble>(m_spacedim * m_spacedim, 0.0);
 
    for (size_t i = 0; i < lenScalesVars.size(); ++i)
    {
        std::string varStr = lenScalesVars[i];
        if (m_session->DefinesFunction(charLenScalesFuncName, varStr))
        {
            clsAux = m_session->GetFunction(charLenScalesFuncName, varStr);
            m_l[i] = (NekDouble)clsAux->Evaluate();
        }
        else
        {
            NEKERROR(Nektar::ErrorUtil::efatal,
                     "Missing parameter '" + varStr +
                         "' in the FUNCTION NAME = " + charLenScalesFuncName +
                         ".");
        }
    }
 
    // Read box of eddies parameters
    m_rc = Array<OneD, NekDouble>(m_spacedim, 0.0);
    // Array<OneD, NekDouble> m_lyz(m_spacedim - 1, 0.0);
    m_lyz       = Array<OneD, NekDouble>(m_spacedim - 1, 0.0);
    elmtInfTurb = pForce->FirstChildElement("BoxOfEddies");
    ASSERTL0(elmtInfTurb,
             "Unable to find BoxOfEddies tag. in SyntheticTurbulence forcing");
 
    if (elmtInfTurb)
    {
        std::stringstream boxStream;
        std::string boxStr = elmtInfTurb->GetText();
        boxStream.str(boxStr);
        size_t countVar = 0;
        for (size_t i = 0; i < (2 * m_spacedim - 1); ++i)
        {
            boxStream >> boxStr;
            if (i < m_spacedim)
            {
                m_rc[i] = boost::lexical_cast<NekDouble>(boxStr);
            }
            else
            {
                m_lyz[i - m_spacedim] = boost::lexical_cast<NekDouble>(boxStr);
            }
            countVar += 1;
        }
        if (countVar != (2 * m_spacedim - 1))
        {
            NEKERROR(Nektar::ErrorUtil::efatal,
                     "Missing parameter in the BoxOfEddies tag");
        }
    }
 
    // Read sigma
    elmtInfTurb = pForce->FirstChildElement("Sigma");
    ASSERTL0(elmtInfTurb,
             "Unable to find Sigma tag. in SyntheticTurbulence forcing");
    std::string sigmaStr = elmtInfTurb->GetText();
    m_sigma              = boost::lexical_cast<NekDouble>(sigmaStr);
 
    // Read bulk velocity
    elmtInfTurb = pForce->FirstChildElement("BulkVelocity");
    ASSERTL0(elmtInfTurb,
             "Unable to find BulkVelocity tag. in SyntheticTurbulence forcing");
    std::string bVelStr = elmtInfTurb->GetText();
    m_Ub                = boost::lexical_cast<NekDouble>(bVelStr);
 
    // Read flag to check if the run is a test case
    elmtInfTurb          = pForce->FirstChildElement("TestCase");
    const char *tcaseStr = (elmtInfTurb) ? elmtInfTurb->GetText() : "NoName";
    m_tCase              = (strcmp(tcaseStr, "ChanFlow3D") == 0) ? true : false;
 
    // Set Cholesky decomposition of the Reynolds Stresses in the domain
    SetCholeskyReyStresses(pFields);
    // Compute reference lengths
    ComputeRefLenghts();
    // Compute Xi maximum
    ComputeXiMax();
    // Set Number of Eddies
    SetNumberOfEddies();
    // Set mask
    SetBoxOfEddiesMask(pFields);
    // Set Forcing for each eddy
    InitialiseForcingEddy(pFields);
    // Check for test case
    if (!m_tCase)
    {
        // Compute initial location of the eddies in the box
        ComputeInitialRandomLocationOfEddies();
    }
    else
    {
        // Compute initial location of the eddies for the test case
        ComputeInitialLocationTestCase();
    }
 
    // Seed to generate random positions for the eddies
    srand(time(nullptr));
 
    // Initialise member from the base class
    m_Forcing = Array<OneD, Array<OneD, NekDouble>>(pFields.size());
    for (int i = 0; i < pFields.size(); ++i)
    {
        m_Forcing[i] = Array<OneD, NekDouble>(pFields[0]->GetTotPoints(), 0.0);
    }
 
    // Initialise eddies ID vector
    for (size_t n = 0; n < m_N; ++n)
    {
        m_eddiesIDForcing.push_back(n);
    }
}

References ASSERTL0, ComputeInitialLocationTestCase(), ComputeInitialRandomLocationOfEddies(), ComputeRefLenghts(), ComputeXiMax(), Nektar::ErrorUtil::efatal, InitialiseForcingEddy(), m_eddiesIDForcing, Nektar::SolverUtils::Forcing::m_Forcing, m_gamma, m_l, m_lyz, m_N, m_R, m_rc, Nektar::SolverUtils::Forcing::m_session, m_sigma, m_spacedim, m_tCase, m_Ub, NEKERROR, SetBoxOfEddiesMask(), SetCholeskyReyStresses(), and SetNumberOfEddies().

Friends And Related Function Documentation

◆ MemoryManager< ForcingSyntheticEddy >

friend class MemoryManager< ForcingSyntheticEddy >

friend

Definition at line 1 of file ForcingSyntheticEddy.h.

Member Data Documentation

◆ className

std::string Nektar::SolverUtils::ForcingSyntheticEddy::className

static

Initial value:

=
    GetForcingFactory().RegisterCreatorFunction(
        "SyntheticTurbulence", ForcingSyntheticEddy::create,
        "Synthetic Eddy Turbulence Method")

Name of the class.

Definition at line 65 of file ForcingSyntheticEddy.h.

◆ m_calcForcing

bool Nektar::SolverUtils::ForcingSyntheticEddy::m_calcForcing {true}

protected

Check when the forcing should be applied.

Definition at line 170 of file ForcingSyntheticEddy.h.

Referenced by UpdateEddiesPositions(), Nektar::SolverUtils::ForcingCFSSyntheticEddy::v_Apply(), Nektar::SolverUtils::ForcingIncNSSyntheticEddy::v_Apply(), and Nektar::SolverUtils::ForcingCFSSyntheticEddy::v_ApplyCoeff().

◆ m_Cholesky

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::ForcingSyntheticEddy::m_Cholesky

protected

Cholesky decomposition of the Reynolds Stresses for all points in the mesh.

Definition at line 142 of file ForcingSyntheticEddy.h.

Referenced by ComputeVelocityFluctuation(), and SetCholeskyReyStresses().

◆ m_currTime

NekDouble Nektar::SolverUtils::ForcingSyntheticEddy::m_currTime

protected

Current time.

Definition at line 174 of file ForcingSyntheticEddy.h.

Referenced by Nektar::SolverUtils::ForcingCFSSyntheticEddy::v_ApplyCoeff().

◆ m_eddiesIDForcing

std::vector<unsigned int> Nektar::SolverUtils::ForcingSyntheticEddy::m_eddiesIDForcing

protected

Eddies that add to the domain.

Definition at line 172 of file ForcingSyntheticEddy.h.

Referenced by Nektar::SolverUtils::ForcingCFSSyntheticEddy::CalculateForcing(), Nektar::SolverUtils::ForcingIncNSSyntheticEddy::CalculateForcing(), Nektar::SolverUtils::ForcingCFSSyntheticEddy::ComputeDensityFluctuation(), ComputeStochasticSignal(), ComputeVelocityFluctuation(), GenerateRandomOneOrMinusOne(), RemoveEddiesFromForcing(), UpdateEddiesPositions(), and v_InitObject().

◆ m_eddyPos

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::ForcingSyntheticEddy::m_eddyPos

protected

Eddy position.

Definition at line 168 of file ForcingSyntheticEddy.h.

Referenced by ComputeInitialLocationTestCase(), ComputeInitialRandomLocationOfEddies(), ComputeStochasticSignal(), and UpdateEddiesPositions().

◆ m_ForcingEddy

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::ForcingSyntheticEddy::m_ForcingEddy

protected

Forcing for each eddy.

Definition at line 178 of file ForcingSyntheticEddy.h.

Referenced by Nektar::SolverUtils::ForcingCFSSyntheticEddy::CalculateForcing(), Nektar::SolverUtils::ForcingIncNSSyntheticEddy::CalculateForcing(), InitialiseForcingEddy(), and RemoveEddiesFromForcing().

◆ m_gamma

NekDouble Nektar::SolverUtils::ForcingSyntheticEddy::m_gamma

protected

Ration of specific heats.

Definition at line 164 of file ForcingSyntheticEddy.h.

Referenced by Nektar::SolverUtils::ForcingCFSSyntheticEddy::ComputeDensityFluctuation(), and v_InitObject().

◆ m_l

Array<OneD, NekDouble> Nektar::SolverUtils::ForcingSyntheticEddy::m_l

protected

Characteristic lenght scales.

Definition at line 154 of file ForcingSyntheticEddy.h.

Referenced by ComputeCharConvTurbTime(), ComputeRefLenghts(), ComputeXiMax(), and v_InitObject().

◆ m_lref

Array<OneD, NekDouble> Nektar::SolverUtils::ForcingSyntheticEddy::m_lref

protected

Reference lenghts.

Definition at line 156 of file ForcingSyntheticEddy.h.

Referenced by ComputeCharConvTurbTime(), ComputeInitialLocationTestCase(), ComputeInitialRandomLocationOfEddies(), ComputeRefLenghts(), ComputeSmoothingFactor(), ComputeStochasticSignal(), ComputeXiMax(), InsideBoxOfEddies(), SetNumberOfEddies(), and UpdateEddiesPositions().

◆ m_lyz

Array<OneD, NekDouble> Nektar::SolverUtils::ForcingSyntheticEddy::m_lyz

protected

Inlet length in the y-direction and z-direction.

Definition at line 148 of file ForcingSyntheticEddy.h.

Referenced by ComputeInitialLocationTestCase(), ComputeInitialRandomLocationOfEddies(), InsideBoxOfEddies(), SetNumberOfEddies(), UpdateEddiesPositions(), and v_InitObject().

◆ m_mask

Array<OneD, int> Nektar::SolverUtils::ForcingSyntheticEddy::m_mask

protected

Box of eddies mask.

Definition at line 166 of file ForcingSyntheticEddy.h.

Referenced by Nektar::SolverUtils::ForcingCFSSyntheticEddy::CalculateForcing(), Nektar::SolverUtils::ForcingIncNSSyntheticEddy::CalculateForcing(), ComputeCharConvTurbTime(), Nektar::SolverUtils::ForcingCFSSyntheticEddy::ComputeDensityFluctuation(), Nektar::SolverUtils::ForcingCFSSyntheticEddy::ComputeRhoMachMean(), ComputeSmoothingFactor(), ComputeStochasticSignal(), ComputeVelocityFluctuation(), and SetBoxOfEddiesMask().

◆ m_N

int Nektar::SolverUtils::ForcingSyntheticEddy::m_N

protected

Number of eddies in the box.

Definition at line 152 of file ForcingSyntheticEddy.h.

Referenced by Nektar::SolverUtils::ForcingCFSSyntheticEddy::ComputeDensityFluctuation(), ComputeInitialLocationTestCase(), ComputeInitialRandomLocationOfEddies(), ComputeStochasticSignal(), ComputeVelocityFluctuation(), GenerateRandomOneOrMinusOne(), InitialiseForcingEddy(), SetNumberOfEddies(), UpdateEddiesPositions(), and v_InitObject().

◆ m_R

std::map<int, LibUtilities::EquationSharedPtr> Nektar::SolverUtils::ForcingSyntheticEddy::m_R

protected

Definition at line 139 of file ForcingSyntheticEddy.h.

Referenced by SetCholeskyReyStresses(), and v_InitObject().

◆ m_rc

Array<OneD, NekDouble> Nektar::SolverUtils::ForcingSyntheticEddy::m_rc

protected

Center of the inlet plane.

Definition at line 150 of file ForcingSyntheticEddy.h.

Referenced by ComputeInitialLocationTestCase(), ComputeInitialRandomLocationOfEddies(), ComputeSmoothingFactor(), InsideBoxOfEddies(), UpdateEddiesPositions(), and v_InitObject().

◆ m_sigma

NekDouble Nektar::SolverUtils::ForcingSyntheticEddy::m_sigma

protected

Standard deviation.

Definition at line 146 of file ForcingSyntheticEddy.h.

Referenced by ComputeGaussian(), and v_InitObject().

◆ m_spacedim

int Nektar::SolverUtils::ForcingSyntheticEddy::m_spacedim

protected

Space dimension.

Definition at line 162 of file ForcingSyntheticEddy.h.

Referenced by Nektar::SolverUtils::ForcingCFSSyntheticEddy::CalculateForcing(), Nektar::SolverUtils::ForcingIncNSSyntheticEddy::CalculateForcing(), ComputeCharConvTurbTime(), ComputeConstantC(), ComputeInitialLocationTestCase(), ComputeInitialRandomLocationOfEddies(), ComputeRefLenghts(), ComputeSmoothingFactor(), ComputeStochasticSignal(), ComputeVelocityFluctuation(), ComputeXiMax(), GenerateRandomOneOrMinusOne(), SetCholeskyReyStresses(), SetNumberOfEddies(), and v_InitObject().

◆ m_tCase

bool Nektar::SolverUtils::ForcingSyntheticEddy::m_tCase

protected

Check for test case.

Definition at line 176 of file ForcingSyntheticEddy.h.

Referenced by GenerateRandomOneOrMinusOne(), UpdateEddiesPositions(), and v_InitObject().

◆ m_Ub

NekDouble Nektar::SolverUtils::ForcingSyntheticEddy::m_Ub

protected

Bulk velocity.

Definition at line 144 of file ForcingSyntheticEddy.h.

Referenced by ComputeCharConvTurbTime(), Nektar::SolverUtils::ForcingCFSSyntheticEddy::ComputeDensityFluctuation(), ComputeSmoothingFactor(), UpdateEddiesPositions(), and v_InitObject().

◆ m_xiMax

Array<OneD, NekDouble> Nektar::SolverUtils::ForcingSyntheticEddy::m_xiMax

protected

Xi max.

Definition at line 158 of file ForcingSyntheticEddy.h.

Referenced by ComputeConstantC(), ComputeStochasticSignal(), and ComputeXiMax().

◆ m_xiMaxMin

NekDouble Nektar::SolverUtils::ForcingSyntheticEddy::m_xiMaxMin = 0.4

protected

Definition at line 160 of file ForcingSyntheticEddy.h.

Referenced by ComputeCharConvTurbTime(), and ComputeXiMax().

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Protected Attributes

Friends

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ ForcingSyntheticEddy()

◆ ~ForcingSyntheticEddy()

Member Function Documentation

◆ ComputeCharConvTurbTime()

◆ ComputeConstantC()

◆ ComputeGaussian()

◆ ComputeInitialLocationTestCase()

◆ ComputeInitialRandomLocationOfEddies()

◆ ComputeRefLenghts()

◆ ComputeSmoothingFactor()

◆ ComputeStochasticSignal()

◆ ComputeVelocityFluctuation()

◆ ComputeXiMax()

◆ create()

◆ GenerateRandomOneOrMinusOne()

◆ InitialiseForcingEddy()

◆ InsideBoxOfEddies()

◆ RemoveEddiesFromForcing()

◆ SetBoxOfEddiesMask()

◆ SetCholeskyReyStresses()

◆ SetNumberOfEddies()

◆ UpdateEddiesPositions()

◆ v_Apply()

◆ v_ApplyCoeff()

◆ v_InitObject()

Friends And Related Function Documentation

◆ MemoryManager< ForcingSyntheticEddy >

Member Data Documentation

◆ className

◆ m_calcForcing

◆ m_Cholesky

◆ m_currTime

◆ m_eddiesIDForcing

◆ m_eddyPos

◆ m_ForcingEddy

◆ m_gamma

◆ m_l

◆ m_lref

◆ m_lyz

◆ m_mask

◆ m_N

◆ m_R

◆ m_rc

◆ m_sigma

◆ m_spacedim

◆ m_tCase

◆ m_Ub

◆ m_xiMax

◆ m_xiMaxMin