Nektar::AdjointAdvection Class Reference

Advection for the adjoint form of the linearised Navier-Stokes equations. More...

#include <AdjointAdvection.h>

Inheritance diagram for Nektar::AdjointAdvection:

Collaboration diagram for Nektar::AdjointAdvection:

Static Public Member Functions
static SolverUtils::AdvectionSharedPtr	create (std::string)
	Creates an instance of this class.

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

Protected Member Functions
DNekBlkMatSharedPtr	GetFloquetBlockMatrix (FloquetMatType mattype, bool UseContCoeffs=false) const
DNekBlkMatSharedPtr	GenFloquetBlockMatrix (FloquetMatType mattype, bool UseContCoeffs=false) const
	AdjointAdvection ()
virtual	~AdjointAdvection ()
virtual void	v_InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
	Initialises the advection object.
virtual void	v_Advect (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &advVel, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time)
	Advects a vector field.
virtual void	v_SetBaseFlow (const Array< OneD, Array< OneD, NekDouble > > &inarray)
	Overrides the base flow used during linearised advection.
void	UpdateBase (const NekDouble m_slices, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble m_time, const NekDouble m_period)
void	DFT (const string file, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const NekDouble m_slices)
void	ImportFldBase (std::string pInfile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, int slice)
	Import Base flow.

Protected Attributes
LibUtilities::SessionReaderSharedPtr	m_session
MultiRegions::ProjectionType	m_projectionType
int	m_spacedim
int	m_expdim
Array< OneD, Array< OneD, NekDouble > >	m_baseflow
	Storage for base flow.
int	m_slices
NekDouble	m_period
Array< OneD, Array< OneD, NekDouble > >	m_interp
LibUtilities::NektarFFTSharedPtr	m_FFT
Array< OneD, NekDouble >	m_tmpIN
Array< OneD, NekDouble >	m_tmpOUT
bool	m_useFFTW
bool	m_SingleMode
	flag to determine if use single mode or not
bool	m_HalfMode
	flag to determine if use half mode or not
bool	m_MultipleModes
	flag to determine if use multiple mode or not
bool	m_homogen_dealiasing
MultiRegions::CoeffState	m_CoeffState
FloquetBlockMatrixMapShPtr	m_FloquetBlockMat
Protected Attributes inherited from Nektar::SolverUtils::Advection
AdvectionFluxVecCB	m_fluxVector
	Callback function to the flux vector (set when advection is in conservative form).
RiemannSolverSharedPtr	m_riemann
	Riemann solver for DG-type schemes.
int	m_spaceDim
	Storage for space dimension. Used for homogeneous extension.

Private Types
enum	FloquetMatType { eForwardsCoeff, eForwardsPhys }
enum	HomogeneousType { eHomogeneous1D, eHomogeneous2D, eHomogeneous3D, eNotHomogeneous }
	Parameter for homogeneous expansions. More...
typedef map< FloquetMatType, DNekBlkMatSharedPtr >	FloquetBlockMatrixMap
	A map between matrix keys and their associated block matrices.
typedef boost::shared_ptr < FloquetBlockMatrixMap >	FloquetBlockMatrixMapShPtr
	A shared pointer to a BlockMatrixMap.

Private Attributes
bool	m_useFFT
	flag to determine if use or not the FFT for transformations
enum HomogeneousType	m_HomogeneousType
NekDouble	m_LhomX
	physical length in X direction (if homogeneous)
NekDouble	m_LhomY
	physical length in Y direction (if homogeneous)
NekDouble	m_LhomZ
	physical length in Z direction (if homogeneous)
int	m_npointsX
	number of points in X direction (if homogeneous)
int	m_npointsY
	number of points in Y direction (if homogeneous)
int	m_npointsZ
	number of points in Z direction (if homogeneous)
int	m_HomoDirec
	number of homogenous directions
int	m_NumMode
	Mode to use in case of single mode analysis.
SpatialDomains::BoundaryConditionsSharedPtr	m_boundaryConditions

Friends
class	MemoryManager< AdjointAdvection >

Additional Inherited Members
Public Member Functions inherited from Nektar::SolverUtils::Advection
SOLVER_UTILS_EXPORT void	InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
	Interface function to initialise the advection object.
SOLVER_UTILS_EXPORT void	Advect (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &advVel, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time)
	Interface function to advect the vector field.
template<typename FuncPointerT , typename ObjectPointerT >
void	SetFluxVector (FuncPointerT func, ObjectPointerT obj)
	Set the flux vector callback function.
void	SetRiemannSolver (RiemannSolverSharedPtr riemann)
	Set a Riemann solver object for this advection object.
void	SetFluxVector (AdvectionFluxVecCB fluxVector)
	Set the flux vector callback function.
void	SetBaseFlow (const Array< OneD, Array< OneD, NekDouble > > &inarray)
	Set the base flow used for linearised advection objects.

Detailed Description

Advection for the adjoint form of the linearised Navier-Stokes equations.

Definition at line 47 of file AdjointAdvection.h.

Member Typedef Documentation

typedef map< FloquetMatType, DNekBlkMatSharedPtr> Nektar::AdjointAdvection::FloquetBlockMatrixMap

private

A map between matrix keys and their associated block matrices.

Definition at line 56 of file AdjointAdvection.h.

typedef boost::shared_ptr<FloquetBlockMatrixMap> Nektar::AdjointAdvection::FloquetBlockMatrixMapShPtr

private

A shared pointer to a BlockMatrixMap.

Definition at line 58 of file AdjointAdvection.h.

Member Enumeration Documentation

enum Nektar::AdjointAdvection::FloquetMatType

private

Enumerator:

eForwardsCoeff
eForwardsPhys

Definition at line 49 of file AdjointAdvection.h.

        {
            eForwardsCoeff,
            eForwardsPhys
        };

enum Nektar::AdjointAdvection::HomogeneousType

private

Parameter for homogeneous expansions.

Enumerator:

eHomogeneous1D
eHomogeneous2D
eHomogeneous3D
eNotHomogeneous

Definition at line 153 of file AdjointAdvection.h.

        {
            eHomogeneous1D,
            eHomogeneous2D,
            eHomogeneous3D,
            eNotHomogeneous
        };

Constructor & Destructor Documentation

Nektar::AdjointAdvection::AdjointAdvection ( )

protected

Definition at line 60 of file AdjointAdvection.cpp.

    : Advection()
{
}

Nektar::AdjointAdvection::~AdjointAdvection ( )

protectedvirtual

Definition at line 258 of file AdjointAdvection.cpp.

{
}

Member Function Documentation

static SolverUtils::AdvectionSharedPtr Nektar::AdjointAdvection::create ( std::string  )

inlinestatic

Creates an instance of this class.

Definition at line 64 of file AdjointAdvection.h.

                                                               {
            return MemoryManager<AdjointAdvection>::AllocateSharedPtr();
        }

void Nektar::AdjointAdvection::DFT ( const string  file, Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields, const NekDouble  m_slices  )

protected

Definition at line 666 of file AdjointAdvection.cpp.

References Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), eForwardsPhys, Nektar::eWrapper, GetFloquetBlockMatrix(), Nektar::LibUtilities::GetNektarFFTFactory(), ImportFldBase(), m_baseflow, m_FFT, m_interp, m_period, m_session, m_slices, m_tmpIN, m_tmpOUT, Vmath::Smul(), Vmath::Vcopy(), and Vmath::Zero().

Referenced by v_InitObject().

{
    int npoints=m_baseflow[0].num_elements();

    //Convected fields
    int ConvectedFields=m_baseflow.num_elements()-1;

    m_interp= Array<OneD, Array<OneD, NekDouble> > (ConvectedFields);
    for(int i=0; i<ConvectedFields;++i)
    {
        m_interp[i]=Array<OneD,NekDouble>(npoints*m_slices);
    }

    //Import the slides into auxiliary vector
    //The base flow should be stored in the form filename_i.bse
    for (int i=0; i< m_slices; ++i)
    {
        char chkout[16] = "";
        sprintf(chkout, "%d", i);
        ImportFldBase(file+"_"+chkout+".bse",pFields,i);
    }


    // Discrete Fourier Transform of the fields
    for(int k=0; k< ConvectedFields;++k)
    {
#ifdef NEKTAR_USING_FFTW

        //Discrete Fourier Transform using FFTW


        Array<OneD, NekDouble> fft_in(npoints*m_slices);
        Array<OneD, NekDouble> fft_out(npoints*m_slices);

        Array<OneD, NekDouble> m_tmpIN(m_slices);
        Array<OneD, NekDouble> m_tmpOUT(m_slices);

        //Shuffle the data
        for(int j= 0; j < m_slices; ++j)
        {
            Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(fft_in[j]),m_slices);
        }

        m_FFT = LibUtilities::GetNektarFFTFactory().CreateInstance("NekFFTW", m_slices);

        //FFT Transform
        for(int i=0; i<npoints; i++)
        {
            m_FFT->FFTFwdTrans(m_tmpIN =fft_in + i*m_slices, m_tmpOUT =fft_out + i*m_slices);

        }

        //Reshuffle data
        for(int s = 0; s < m_slices; ++s)
        {
            Vmath::Vcopy(npoints,&fft_out[s],m_slices,&m_interp[k][s*npoints],1);

        }

        Vmath::Zero(fft_in.num_elements(),&fft_in[0],1);
        Vmath::Zero(fft_out.num_elements(),&fft_out[0],1);
#else
        //Discrete Fourier Transform using MVM


        DNekBlkMatSharedPtr blkmat;
        blkmat = GetFloquetBlockMatrix(eForwardsPhys);

        int nrows = blkmat->GetRows();
        int ncols = blkmat->GetColumns();

        Array<OneD, NekDouble> sortedinarray(ncols);
        Array<OneD, NekDouble> sortedoutarray(nrows);

        //Shuffle the data
        for(int j= 0; j < m_slices; ++j)
        {
            Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(sortedinarray[j]),m_slices);
        }

        // Create NekVectors from the given data arrays
        NekVector<NekDouble> in (ncols,sortedinarray,eWrapper);
        NekVector<NekDouble> out(nrows,sortedoutarray,eWrapper);

        // Perform matrix-vector multiply.
        out = (*blkmat)*in;

        //Reshuffle data
        for(int s = 0; s < m_slices; ++s)
        {
            Vmath::Vcopy(npoints,&sortedoutarray[s],m_slices,&m_interp[k][s*npoints],1);
        }

        Vmath::Zero(sortedinarray.num_elements(),&sortedinarray[0],1);
        Vmath::Zero(sortedoutarray.num_elements(),&sortedoutarray[0],1);

#endif

        //scaling of the Fourier coefficients
        NekDouble j=-1;
        for (int i = 2; i < m_slices; i += 2)
        {
            Vmath::Smul(2*npoints,j,&m_interp[k][i*npoints],1,&m_interp[k][i*npoints],1);
            j=-j;

        }

    }

    if(m_session->DefinesParameter("period"))
    {
        m_period=m_session->GetParameter("period");
    }
    else
    {
        m_period=(m_session->GetParameter("TimeStep")*m_slices)/(m_slices-1.);
    }
}

DNekBlkMatSharedPtr Nektar::AdjointAdvection::GenFloquetBlockMatrix ( FloquetMatType  mattype, bool  UseContCoeffs = false  ) const

protected

DNekBlkMatSharedPtr Nektar::AdjointAdvection::GetFloquetBlockMatrix ( FloquetMatType  mattype, bool  UseContCoeffs = false  ) const

protected

Definition at line 627 of file AdjointAdvection.cpp.

References Nektar::StdRegions::StdExpansion::DetShapeType(), Nektar::eDIAGONAL, Nektar::LibUtilities::eFourier, Nektar::LibUtilities::eFourierEvenlySpaced, Nektar::StdRegions::eFwdTrans, Nektar::StdRegions::StdExpansion::GetStdMatrix(), m_baseflow, and m_slices.

Referenced by DFT().

{
    DNekMatSharedPtr    loc_mat;
    DNekBlkMatSharedPtr BlkMatrix;
    int n_exp = 0;

    n_exp = m_baseflow[0].num_elements(); // will operatore on m_phys

    Array<OneD,unsigned int> nrows(n_exp);
    Array<OneD,unsigned int> ncols(n_exp);

    nrows = Array<OneD, unsigned int>(n_exp,m_slices);
    ncols = Array<OneD, unsigned int>(n_exp,m_slices);

    MatrixStorage blkmatStorage = eDIAGONAL;
    BlkMatrix = MemoryManager<DNekBlkMat>
    ::AllocateSharedPtr(nrows,ncols,blkmatStorage);


    const LibUtilities::PointsKey Pkey(m_slices,LibUtilities::eFourierEvenlySpaced);
    const LibUtilities::BasisKey  BK(LibUtilities::eFourier,m_slices,Pkey);
    StdRegions::StdSegExp StdSeg(BK);

    StdRegions::StdMatrixKey matkey(StdRegions::eFwdTrans,
                                                StdSeg.DetShapeType(),
                                                StdSeg);

    loc_mat = StdSeg.GetStdMatrix(matkey);

    // set up array of block matrices.
    for(int i = 0; i < n_exp; ++i)
    {
        BlkMatrix->SetBlock(i,i,loc_mat);
    }

    return BlkMatrix;
}

void Nektar::AdjointAdvection::ImportFldBase ( std::string  pInfile, Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields, int  slice  )

protected

Import Base flow.

Import field from infile and load into m_fields. This routine will also perform a BwdTrans to ensure data is in both the physical and coefficient storage.

Parameters

infile

Filename to read.

Definition at line 497 of file AdjointAdvection.cpp.

References ASSERTL0, m_baseflow, m_HalfMode, m_interp, m_npointsZ, m_session, m_SingleMode, Vmath::Vcopy(), and Vmath::Zero().

Referenced by DFT(), and v_InitObject().

{
    std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
    std::vector<std::vector<NekDouble> > FieldData;

    int nqtot = pFields[0]->GetTotPoints();
    Array<OneD, NekDouble> tmp_coeff(pFields[0]->GetNcoeffs(), 0.0);

    //Get Homogeneous
    LibUtilities::FieldIOSharedPtr fld =
            MemoryManager<LibUtilities::FieldIO>::AllocateSharedPtr(
                                                    m_session->GetComm());
    fld->Import(pInfile, FieldDef, FieldData);

    int nvar = m_session->GetVariables().size();
    int s;

    if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
    {
        std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
    }
    
    // copy FieldData into m_fields
    for(int j = 0; j < nvar; ++j)
    {
        for(int i = 0; i < FieldDef.size(); ++i)
        {
            if ((m_session->DefinesSolverInfo("HOMOGENEOUS") &&
                (m_session->GetSolverInfo("HOMOGENEOUS")=="HOMOGENEOUS1D" ||
                 m_session->GetSolverInfo("HOMOGENEOUS")=="1D" ||
                 m_session->GetSolverInfo("HOMOGENEOUS")=="Homo1D")) &&
                nvar==4)
            {
                // w-component must be ignored and set to zero.
                if (j != nvar - 2)
                {
                    // p component (it is 4th variable of the 3D and corresponds 3nd variable of 2D)
                    s = (j == nvar - 1) ? 2 : j;

                    //extraction of the 2D
                    pFields[j]->ExtractDataToCoeffs(
                                        FieldDef[i],
                                        FieldData[i],
                                        FieldDef[i]->m_fields[s],
                                        tmp_coeff);
                }

                // Put zero on higher modes
                int ncplane = (pFields[0]->GetNcoeffs()) / m_npointsZ;

                if (m_npointsZ > 2)
                {
                    Vmath::Zero(ncplane*(m_npointsZ-2),
                                &tmp_coeff[2*ncplane],1);
                }
            }
            //2D cases and Homogeneous1D Base Flows
            else
            {
                bool flag = FieldDef[i]->m_fields[j] ==
                                                m_session->GetVariable(j);

                ASSERTL0(flag, (std::string("Order of ") + pInfile
                                + std::string(" data and that defined in "
                                              "m_boundaryconditions differs")).c_str());

                pFields[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
                                               FieldDef[i]->m_fields[j],
                                               tmp_coeff);
            }
        }

        if(m_SingleMode || m_HalfMode)
        {
            pFields[j]->GetPlane(0)->BwdTrans(tmp_coeff, m_baseflow[j]);

            if(m_SingleMode)
            {
                //copy the bwd into the second plane for single Mode Analysis
                int ncplane=(pFields[0]->GetNpoints())/m_npointsZ;
                Vmath::Vcopy(ncplane,&m_baseflow[j][0],1,&m_baseflow[j][ncplane],1);
            }
        }
        else
        {
            pFields[j]->BwdTrans(tmp_coeff, m_baseflow[j]);
        }
    }

    if(m_session->DefinesParameter("N_slices"))
    {
        int n = pFields.num_elements()-1;

        for(int i=0; i<n;++i)
        {
            Vmath::Vcopy(nqtot, &m_baseflow[i][0], 1, &m_interp[i][slice*nqtot], 1);
        }

    }
}

void Nektar::AdjointAdvection::UpdateBase ( const NekDouble  m_slices, const Array< OneD, const NekDouble > &  inarray, Array< OneD, NekDouble > &  outarray, const NekDouble  m_time, const NekDouble  m_period  )

protected

Definition at line 600 of file AdjointAdvection.cpp.

References m_baseflow, m_period, m_slices, Vmath::Svtvp(), and Vmath::Vcopy().

Referenced by v_Advect().

{

    int npoints=m_baseflow[0].num_elements();

    NekDouble BetaT=2*M_PI*fmod (m_time, m_period) / m_period;
    NekDouble phase;
    Array<OneD, NekDouble> auxiliary(npoints);

    Vmath::Vcopy(npoints,&inarray[0],1,&outarray[0],1);
    Vmath::Svtvp(npoints, cos(0.5*m_slices*BetaT),&inarray[npoints],1,&outarray[0],1,&outarray[0],1);

    for (int i = 2; i < m_slices; i += 2)
    {
        phase = (i>>1) * BetaT;

        Vmath::Svtvp(npoints, cos(phase),&inarray[i*npoints],1,&outarray[0],1,&outarray[0],1);
        Vmath::Svtvp(npoints, sin(phase), &inarray[(i+1)*npoints], 1, &outarray[0], 1,&outarray[0],1);
    }

}

void Nektar::AdjointAdvection::v_Advect ( const int  nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields, const Array< OneD, Array< OneD, NekDouble > > &  advVel, const Array< OneD, Array< OneD, NekDouble > > &  inarray, Array< OneD, Array< OneD, NekDouble > > &  outarray, const NekDouble &  time  )

protectedvirtual

Advects a vector field.

Implements Nektar::SolverUtils::Advection.

Definition at line 263 of file AdjointAdvection.cpp.

References ASSERTL0, ASSERTL1, Nektar::LibUtilities::eFunctionTypeFile, m_baseflow, m_HalfMode, m_interp, m_period, m_session, m_SingleMode, m_slices, Vmath::Neg(), UpdateBase(), Vmath::Vmul(), and Vmath::Vvtvp().

{
    int nqtot            = fields[0]->GetTotPoints();
    ASSERTL1(nConvectiveFields == inarray.num_elements(),"Number of convective fields and Inarray are not compatible");

    Array<OneD, NekDouble > Deriv = Array<OneD, NekDouble> (nqtot*nConvectiveFields);

    for(int n = 0; n < nConvectiveFields; ++n)
    {
        //v_ComputeAdvectionTerm(fields,advVel,inarray[i],outarray[i],i,time,Deriv);
        int ndim       = advVel.num_elements();
        int nPointsTot = fields[0]->GetNpoints();
        Array<OneD, NekDouble> grad0,grad1,grad2;

        //Evaluation of the gradiend of each component of the base flow
        //\nabla U
        Array<OneD, NekDouble> grad_base_u0,grad_base_u1,grad_base_u2;
        // \nabla V
        Array<OneD, NekDouble> grad_base_v0,grad_base_v1,grad_base_v2;
        // \nabla W
        Array<OneD, NekDouble> grad_base_w0,grad_base_w1,grad_base_w2;


        grad0 = Array<OneD, NekDouble> (nPointsTot);
        grad_base_u0 = Array<OneD, NekDouble> (nPointsTot);
        grad_base_v0 = Array<OneD, NekDouble> (nPointsTot);
        grad_base_w0 = Array<OneD, NekDouble> (nPointsTot);


        //Evaluation of the base flow for periodic cases
        //(it requires fld files)
        if(m_slices>1)
        {
            if (m_session->GetFunctionType("BaseFlow", 0)
                == LibUtilities::eFunctionTypeFile)
            {
                for(int i=0; i<ndim;++i)
                {
                    UpdateBase(m_slices,m_interp[i],m_baseflow[i],time,m_period);
                }
            }
            else
            {
                ASSERTL0(false, "Periodic Base flow requires .fld files");
            }
        }

        //Evaluate the Adjoint advection term
        switch(ndim)
        {
                    // 1D
                case 1:
                    fields[0]->PhysDeriv(inarray[n],grad0);
                    fields[0]->PhysDeriv(m_baseflow[0],grad_base_u0);
                    //Evaluate  U du'/dx
                    Vmath::Vmul(nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
                    //Evaluate U du'/dx+ u' dU/dx
                    Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
                    break;

                    //2D
                case 2:

                    grad1 = Array<OneD, NekDouble> (nPointsTot);
                    grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
                    grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);

                    fields[0]->PhysDeriv(inarray[n],grad0,grad1);

                    //Derivates of the base flow
                    fields[0]-> PhysDeriv(m_baseflow[0], grad_base_u0, grad_base_u1);
                    fields[0]-> PhysDeriv(m_baseflow[1], grad_base_v0, grad_base_v1);

                    //Since the components of the velocity are passed one by one, it is necessary to distinguish which
                    //term is consider
                    switch (n)
                    {
                        //x-equation
                    case 0:
                        // Evaluate U du'/dx
                        Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
                        //Evaluate U du'/dx+ V du'/dy
                        Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
                        //Evaluate - (U du'/dx+ V du'/dy)
                        Vmath::Neg(nPointsTot,outarray[n],1);
                        //Evaluate -(U du'/dx+ V du'/dy)+u' dU/dx
                        Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
                        //Evaluate -(U du'/dx+ V du'/dy) +u' dU/dx +v' dV/dx
                        Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[1],1,outarray[n],1,outarray[n],1);
                        break;

                        //y-equation
                    case 1:
                        // Evaluate U dv'/dx
                        Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
                        //Evaluate U dv'/dx+ V dv'/dy
                        Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
                        //Evaluate -(U dv'/dx+ V dv'/dy)
                        Vmath::Neg(nPointsTot,outarray[n],1);
                        //Evaluate (U dv'/dx+ V dv'/dy)+u' dU/dy
                        Vmath::Vvtvp(nPointsTot,grad_base_u1,1,advVel[0],1,outarray[n],1,outarray[n],1);
                        //Evaluate (U dv'/dx+ V dv'/dy +u' dv/dx)+v' dV/dy
                        Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,outarray[n],1,outarray[n],1);
                        break;
                }
                    break;


                    //3D
                case 3:

                    grad1 = Array<OneD, NekDouble> (nPointsTot);
                    grad2 = Array<OneD, NekDouble> (nPointsTot);
                    grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
                    grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
                    grad_base_w1 = Array<OneD, NekDouble> (nPointsTot);

                    grad_base_u2 = Array<OneD, NekDouble> (nPointsTot);
                    grad_base_v2 = Array<OneD, NekDouble> (nPointsTot);
                    grad_base_w2 = Array<OneD, NekDouble> (nPointsTot);

                    fields[0]->PhysDeriv(m_baseflow[0], grad_base_u0, grad_base_u1,grad_base_u2);
                    fields[0]->PhysDeriv(m_baseflow[1], grad_base_v0, grad_base_v1,grad_base_v2);
                    fields[0]->PhysDeriv(m_baseflow[2], grad_base_w0, grad_base_w1, grad_base_w2);

                    //HalfMode has W(x,y,t)=0
                    if(m_HalfMode || m_SingleMode)
                    {
                        for(int i=0; i<grad_base_u2.num_elements();++i)
                        {
                            grad_base_u2[i]=0;
                            grad_base_v2[i]=0;
                            grad_base_w2[i]=0;

                        }
                    }

                fields[0]->PhysDeriv(inarray[n], grad0, grad1, grad2);

                switch (n)
                {
                        //x-equation
                    case 0:
                        //Evaluate U du'/dx
                        Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
                        //Evaluate U du'/dx+ V du'/dy
                        Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
                        //Evaluate U du'/dx+ V du'/dy+W du'/dz
                        Vmath::Vvtvp(nPointsTot,grad2,1,m_baseflow[2],1,outarray[n],1,outarray[n],1);
                        //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)
                        Vmath::Neg(nPointsTot,outarray[n],1);
                        //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)+u' dU/dx
                        Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
                        //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)+u'dU/dx+ v' dV/dx
                        Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[1],1,outarray[n],1,outarray[n],1);
                        //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)+u'dU/dx+ v' dV/dx+ w' dW/dz
                        Vmath::Vvtvp(nPointsTot,grad_base_w0,1,advVel[2],1,outarray[n],1,outarray[n],1);
                        break;
                        //y-equation
                    case 1:
                        //Evaluate U dv'/dx
                        Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
                        //Evaluate U dv'/dx+ V dv'/dy
                        Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
                        //Evaluate U dv'/dx+ V dv'/dy+W dv'/dz
                        Vmath::Vvtvp(nPointsTot,grad2,1,m_baseflow[2],1,outarray[n],1,outarray[n],1);
                        //Evaluate -(U dv'/dx+ V dv'/dy+W dv'/dz)
                        Vmath::Neg(nPointsTot,outarray[n],1);
                        //Evaluate  -(U dv'/dx+ V dv'/dy+W dv'/dz)+u' dU/dy
                        Vmath::Vvtvp(nPointsTot,grad_base_u1,1,advVel[0],1,outarray[n],1,outarray[n],1);
                        //Evaluate  -(U dv'/dx+ V dv'/dy+W dv'/dz)+u' dU/dy +v' dV/dy
                        Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,outarray[n],1,outarray[n],1);
                        //Evaluate  -(U dv'/dx+ V dv'/dy+W dv'/dz)+u' dU/dy +v' dV/dy+ w' dW/dy
                        Vmath::Vvtvp(nPointsTot,grad_base_w1,1,advVel[2],1,outarray[n],1,outarray[n],1);
                        break;

                        //z-equation
                    case 2:
                        //Evaluate U dw'/dx
                        Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
                        //Evaluate U dw'/dx+ V dw'/dx
                        Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
                        //Evaluate U dw'/dx+ V dw'/dx+ W dw'/dz
                        Vmath::Vvtvp(nPointsTot,grad2,1,m_baseflow[2],1,outarray[n],1,outarray[n],1);
                        //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)
                        Vmath::Neg(nPointsTot,outarray[n],1);
                        //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)+u' dU/dz
                        Vmath::Vvtvp(nPointsTot,grad_base_u2,1,advVel[0],1,outarray[n],1,outarray[n],1);
                        //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)+u' dU/dz+v'dV/dz
                        Vmath::Vvtvp(nPointsTot,grad_base_v2,1,advVel[1],1,outarray[n],1,outarray[n],1);
                        //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)+u' dU/dz+v'dV/dz + w' dW/dz
                        Vmath::Vvtvp(nPointsTot,grad_base_w2,1,advVel[2],1,outarray[n],1,outarray[n],1);
                        break;
                }
                    break;
        default:
            ASSERTL0(false,"dimension unknown");
        }

        Vmath::Neg(nqtot,outarray[n],1);
    }

}

void Nektar::AdjointAdvection::v_InitObject ( LibUtilities::SessionReaderSharedPtr  pSession, Array< OneD, MultiRegions::ExpListSharedPtr >  pFields  )

protectedvirtual

Initialises the advection object.

This function should be overridden in derived classes to initialise the specific advection data members. However, this base class function should be called as the first statement of the overridden function to ensure the base class is correctly initialised in order.

Parameters

pSession	Session information.
pFields	Expansion lists for scalar fields.

Reimplemented from Nektar::SolverUtils::Advection.

Definition at line 66 of file AdjointAdvection.cpp.

References ASSERTL0, DFT(), Nektar::MultiRegions::eDiscontinuous, Nektar::LibUtilities::eFunctionTypeFile, Nektar::MultiRegions::eGalerkin, eHomogeneous1D, eHomogeneous2D, eHomogeneous3D, Nektar::MultiRegions::eLocal, eNotHomogeneous, ImportFldBase(), m_baseflow, m_boundaryConditions, m_CoeffState, m_expdim, m_HalfMode, m_HomoDirec, m_homogen_dealiasing, m_HomogeneousType, m_LhomX, m_LhomY, m_LhomZ, m_MultipleModes, m_npointsX, m_npointsY, m_npointsZ, m_projectionType, m_session, m_SingleMode, m_slices, m_spacedim, and m_useFFT.

{

    Advection::v_InitObject(pSession, pFields);

    m_session            = pSession;
    m_boundaryConditions = MemoryManager<SpatialDomains::BoundaryConditions>
                    ::AllocateSharedPtr(m_session, pFields[0]->GetGraph());
    m_spacedim           = pFields[0]->GetGraph()->GetSpaceDimension();
    m_expdim             = pFields[0]->GetGraph()->GetMeshDimension();
    m_CoeffState         = MultiRegions::eLocal;

    //Setting parameters for homogeneous problems
    m_HomoDirec          = 0;
    m_useFFT             = false;
    m_HomogeneousType    = eNotHomogeneous;
    m_SingleMode         = false;
    m_HalfMode           = false;
    m_MultipleModes      = false;
    m_homogen_dealiasing = false;

    if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
    {
        std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
        m_spacedim          = 3;

        if((HomoStr == "HOMOGENEOUS1D")||(HomoStr == "Homogeneous1D")||
           (HomoStr == "1D")||(HomoStr == "Homo1D"))
        {
            m_HomogeneousType = eHomogeneous1D;
            m_LhomZ           = m_session->GetParameter("LZ");
            m_HomoDirec       = 1;

            m_homogen_dealiasing = pSession->DefinesSolverInfo("DEALIASING");

            if(m_session->DefinesSolverInfo("ModeType"))
            {
                m_session->MatchSolverInfo("ModeType","SingleMode",m_SingleMode,false);
                m_session->MatchSolverInfo("ModeType","HalfMode",m_HalfMode,false);
                m_session->MatchSolverInfo("ModeType","MultipleModes",m_MultipleModes,false);
            }
            if(m_session->DefinesSolverInfo("ModeType"))
            {
                if(m_SingleMode)
                {
                    m_npointsZ=2;
                }
                else if(m_HalfMode)
                {
                    m_npointsZ=1;
                }
                else if(m_MultipleModes)
                {
                    m_npointsZ        = m_session->GetParameter("HomModesZ");
                }
                else
                {
                    ASSERTL0(false, "SolverInfo ModeType not valid");
                }
            }
            else
            {
                m_session->LoadParameter("HomModesZ",m_npointsZ);

            }

        }

        if((HomoStr == "HOMOGENEOUS2D")||(HomoStr == "Homogeneous2D")||
           (HomoStr == "2D")||(HomoStr == "Homo2D"))
        {
            m_HomogeneousType = eHomogeneous2D;
            m_session->LoadParameter("HomModesY", m_npointsY);
            m_session->LoadParameter("LY",        m_LhomY);
            m_session->LoadParameter("HomModesZ", m_npointsZ);
            m_session->LoadParameter("LZ",        m_LhomZ);
            m_HomoDirec       = 2;
        }

        if((HomoStr == "HOMOGENEOUS3D")||(HomoStr == "Homogeneous3D")||
           (HomoStr == "3D")||(HomoStr == "Homo3D"))
        {
            m_HomogeneousType = eHomogeneous3D;
            m_session->LoadParameter("HomModesX",m_npointsX);
            m_session->LoadParameter("LX",       m_LhomX   );
            m_session->LoadParameter("HomModesY",m_npointsY);
            m_session->LoadParameter("LY",       m_LhomY   );
            m_session->LoadParameter("HomModesZ",m_npointsZ);
            m_session->LoadParameter("LZ",       m_LhomZ   );
            m_HomoDirec       = 3;
        }

        if(m_session->DefinesSolverInfo("USEFFT"))
        {
            m_useFFT = true;
        }
    }
    else
    {
        m_npointsZ = 1; // set to default value so can use to identify 2d or 3D (homogeneous) expansions
    }

    if(m_session->DefinesSolverInfo("PROJECTION"))
    {
        std::string ProjectStr
        = m_session->GetSolverInfo("PROJECTION");

        if((ProjectStr == "Continuous")||(ProjectStr == "Galerkin")||
           (ProjectStr == "CONTINUOUS")||(ProjectStr == "GALERKIN"))
        {
            m_projectionType = MultiRegions::eGalerkin;
        }
        else if(ProjectStr == "DisContinuous")
        {
            m_projectionType = MultiRegions::eDiscontinuous;
        }
        else
        {
            ASSERTL0(false,"PROJECTION value not recognised");
        }
    }
    else
    {
        cerr << "Projection type not specified in SOLVERINFO,"
        "defaulting to continuous Galerkin" << endl;
        m_projectionType = MultiRegions::eGalerkin;
    }

    int nvar = m_session->GetVariables().size();
    m_baseflow = Array<OneD, Array<OneD, NekDouble> >(nvar);
    for (int i = 0; i < nvar; ++i)
    {
        m_baseflow[i] = Array<OneD, NekDouble>(pFields[i]->GetTotPoints(), 0.0);
    }
    //SetUpBaseFields(pFields[0]->GetGraph());
    ASSERTL0(m_session->DefinesFunction("BaseFlow"),
             "Base flow must be defined for linearised forms.");
    string file = m_session->GetFunctionFilename("BaseFlow", 0);


    //Periodic base flows
    if(m_session->DefinesParameter("N_slices"))
    {
        m_session->LoadParameter("N_slices",m_slices);
        if(m_slices>1)
        {
            DFT(file,pFields,m_slices);
        }
        else
        {
            ASSERTL0(false,"Number of slices must be a positive number greater than 1");
        }
    }
    //Steady base-flow
    else
    {
        m_slices=1;

        //BaseFlow from file
        if (m_session->GetFunctionType("BaseFlow", m_session->GetVariable(0))
            == LibUtilities::eFunctionTypeFile)
        {
            ImportFldBase(file,pFields,1);

        }
        //analytic base flow
        else
        {
            int nq = pFields[0]->GetNpoints();
            Array<OneD,NekDouble> x0(nq);
            Array<OneD,NekDouble> x1(nq);
            Array<OneD,NekDouble> x2(nq);

            // get the coordinates (assuming all fields have the same
            // discretisation)
            pFields[0]->GetCoords(x0,x1,x2);
            for(unsigned int i = 0 ; i < m_baseflow.num_elements(); i++)
            {
                LibUtilities::EquationSharedPtr ifunc
                = m_session->GetFunction("BaseFlow", i);

                ifunc->Evaluate(x0,x1,x2,m_baseflow[i]);
            }

        }

    }
}

void Nektar::AdjointAdvection::v_SetBaseFlow ( const Array< OneD, Array< OneD, NekDouble > > &  inarray )

protectedvirtual

Overrides the base flow used during linearised advection.

Reimplemented from Nektar::SolverUtils::Advection.

Definition at line 474 of file AdjointAdvection.cpp.

References ASSERTL1, m_baseflow, npts, and Vmath::Vcopy().

{
    ASSERTL1(inarray.num_elements() == m_baseflow.num_elements(),
             "Number of base flow variables does not match what is"
             "expected.");

    int npts = inarray[0].num_elements();
    for (int i = 0; i < inarray.num_elements(); ++i)
    {
        ASSERTL1(npts == m_baseflow.num_elements(),
                "Size of base flow array does not match expected.");
        Vmath::Vcopy(npts, inarray[i], 1, m_baseflow[i], 1);
    }
}

Friends And Related Function Documentation

friend class MemoryManager< AdjointAdvection >

friend

Definition at line 61 of file AdjointAdvection.h.

Member Data Documentation

string Nektar::AdjointAdvection::className

static

Initial value:

 SolverUtils
        ::GetAdvectionFactory().RegisterCreatorFunction("Adjoint",
                                                    AdjointAdvection::create)

Name of class.

Definition at line 69 of file AdjointAdvection.h.

Array<OneD, Array<OneD, NekDouble> > Nektar::AdjointAdvection::m_baseflow

protected

Storage for base flow.

Definition at line 79 of file AdjointAdvection.h.

Referenced by DFT(), GetFloquetBlockMatrix(), ImportFldBase(), UpdateBase(), v_Advect(), v_InitObject(), and v_SetBaseFlow().

SpatialDomains::BoundaryConditionsSharedPtr Nektar::AdjointAdvection::m_boundaryConditions

private

Definition at line 178 of file AdjointAdvection.h.

Referenced by v_InitObject().

MultiRegions::CoeffState Nektar::AdjointAdvection::m_CoeffState

protected

Definition at line 100 of file AdjointAdvection.h.

Referenced by v_InitObject().

int Nektar::AdjointAdvection::m_expdim

protected

Definition at line 76 of file AdjointAdvection.h.

Referenced by v_InitObject().

LibUtilities::NektarFFTSharedPtr Nektar::AdjointAdvection::m_FFT

protected

Definition at line 88 of file AdjointAdvection.h.

Referenced by DFT().

FloquetBlockMatrixMapShPtr Nektar::AdjointAdvection::m_FloquetBlockMat

protected

Definition at line 110 of file AdjointAdvection.h.

bool Nektar::AdjointAdvection::m_HalfMode

protected

flag to determine if use half mode or not

Definition at line 96 of file AdjointAdvection.h.

Referenced by ImportFldBase(), v_Advect(), and v_InitObject().

int Nektar::AdjointAdvection::m_HomoDirec

private

number of homogenous directions

Definition at line 174 of file AdjointAdvection.h.

Referenced by v_InitObject().

bool Nektar::AdjointAdvection::m_homogen_dealiasing

protected

Definition at line 99 of file AdjointAdvection.h.

Referenced by v_InitObject().

enum HomogeneousType Nektar::AdjointAdvection::m_HomogeneousType

private

Definition at line 164 of file AdjointAdvection.h.

Referenced by v_InitObject().

Array<OneD, Array<OneD, NekDouble> > Nektar::AdjointAdvection::m_interp

protected

Definition at line 86 of file AdjointAdvection.h.

Referenced by DFT(), ImportFldBase(), and v_Advect().

NekDouble Nektar::AdjointAdvection::m_LhomX

private

physical length in X direction (if homogeneous)

Definition at line 166 of file AdjointAdvection.h.

Referenced by v_InitObject().

NekDouble Nektar::AdjointAdvection::m_LhomY

private

physical length in Y direction (if homogeneous)

Definition at line 167 of file AdjointAdvection.h.

Referenced by v_InitObject().

NekDouble Nektar::AdjointAdvection::m_LhomZ

private

physical length in Z direction (if homogeneous)

Definition at line 168 of file AdjointAdvection.h.

Referenced by v_InitObject().

bool Nektar::AdjointAdvection::m_MultipleModes

protected

flag to determine if use multiple mode or not

Definition at line 98 of file AdjointAdvection.h.

Referenced by v_InitObject().

int Nektar::AdjointAdvection::m_npointsX

private

number of points in X direction (if homogeneous)

Definition at line 170 of file AdjointAdvection.h.

Referenced by v_InitObject().

int Nektar::AdjointAdvection::m_npointsY

private

number of points in Y direction (if homogeneous)

Definition at line 171 of file AdjointAdvection.h.

Referenced by v_InitObject().

int Nektar::AdjointAdvection::m_npointsZ

private

number of points in Z direction (if homogeneous)

Definition at line 172 of file AdjointAdvection.h.

Referenced by ImportFldBase(), and v_InitObject().

int Nektar::AdjointAdvection::m_NumMode

private

Mode to use in case of single mode analysis.

Definition at line 176 of file AdjointAdvection.h.

NekDouble Nektar::AdjointAdvection::m_period

protected

Definition at line 84 of file AdjointAdvection.h.

Referenced by DFT(), UpdateBase(), and v_Advect().

MultiRegions::ProjectionType Nektar::AdjointAdvection::m_projectionType

protected

Definition at line 74 of file AdjointAdvection.h.

Referenced by v_InitObject().

LibUtilities::SessionReaderSharedPtr Nektar::AdjointAdvection::m_session

protected

Definition at line 72 of file AdjointAdvection.h.

Referenced by DFT(), ImportFldBase(), v_Advect(), and v_InitObject().

bool Nektar::AdjointAdvection::m_SingleMode

protected

flag to determine if use single mode or not

Definition at line 94 of file AdjointAdvection.h.

Referenced by ImportFldBase(), v_Advect(), and v_InitObject().

int Nektar::AdjointAdvection::m_slices

protected

Definition at line 82 of file AdjointAdvection.h.

Referenced by DFT(), GetFloquetBlockMatrix(), UpdateBase(), v_Advect(), and v_InitObject().

int Nektar::AdjointAdvection::m_spacedim

protected

Definition at line 75 of file AdjointAdvection.h.

Referenced by v_InitObject().

Array<OneD,NekDouble> Nektar::AdjointAdvection::m_tmpIN

protected

Definition at line 89 of file AdjointAdvection.h.

Referenced by DFT().

Array<OneD,NekDouble> Nektar::AdjointAdvection::m_tmpOUT

protected

Definition at line 90 of file AdjointAdvection.h.

Referenced by DFT().

bool Nektar::AdjointAdvection::m_useFFT

private

flag to determine if use or not the FFT for transformations

Definition at line 162 of file AdjointAdvection.h.

Referenced by v_InitObject().

bool Nektar::AdjointAdvection::m_useFFTW

protected

Definition at line 91 of file AdjointAdvection.h.