Nektar::ForcingMovingBody Class Reference

#include <ForcingMovingBody.h>

Inheritance diagram for Nektar::ForcingMovingBody:

Static Public Member Functions
static SolverUtils::ForcingSharedPtr	create (const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::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 std::string	className
	Name of the class. More...

Protected Member Functions
virtual void	v_InitObject (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce) override
virtual 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) override
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
GlobalMapping::MappingSharedPtr	m_mapping
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...

Private Member Functions
	ForcingMovingBody (const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::EquationSystem > &pEquation)
void	CheckIsFromFile (const TiXmlElement *pForce)
void	InitialiseCableModel (const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
void	InitialiseFilter (const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pForce)
void	Newmark_betaSolver (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &FcePhysinArray, Array< OneD, NekDouble > &MotPhysinArray)
void	EvaluateStructDynModel (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &Hydroforces, NekDouble time)
void	SetDynEqCoeffMatrix (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
void	RollOver (Array< OneD, Array< OneD, NekDouble > > &input)

Private Attributes
size_t	m_movingBodyCalls
	number of times the movbody have been called More...
size_t	m_np
	number of planes per processors More...
size_t	m_vdim
	vibration dimension More...
NekDouble	m_structrho
	mass of the cable per unit length More...
NekDouble	m_structdamp
	damping ratio of the cable More...
NekDouble	m_lhom
	length ratio of the cable More...
NekDouble	m_kinvis
	fluid viscous More...
NekDouble	m_timestep
	time step More...
LibUtilities::NektarFFTSharedPtr	m_FFT
FilterMovingBodySharedPtr	m_MovBodyfilter
Array< OneD, NekDouble >	m_Aeroforces
	storage for the cable's force(x,y) variables More...
Array< OneD, Array< OneD, NekDouble > >	m_MotionVars
	storage for the cable's motion(x,y) variables More...
Array< OneD, Array< OneD, Array< OneD, NekDouble > > >	m_fV
	fictitious velocity storage More...
Array< OneD, Array< OneD, Array< OneD, NekDouble > > >	m_fA
	fictitious acceleration storage More...
Array< OneD, DNekMatSharedPtr >	m_CoeffMat_A
	matrices in Newmart-beta method More...
Array< OneD, DNekMatSharedPtr >	m_CoeffMat_B
	matrices in Newmart-beta method More...
Array< OneD, std::string >	m_funcName
	[0] is displacements, [1] is velocities, [2] is accelerations More...
Array< OneD, std::string >	m_motion
	motion direction: [0] is 'x' and [1] is 'y' More...
Array< OneD, bool >	m_IsFromFile
	do determine if the the body motion come from an extern file More...
Array< OneD, Array< OneD, NekDouble > >	m_zta
	Store the derivatives of motion variables in x-direction. More...
Array< OneD, Array< OneD, NekDouble > >	m_eta
	Store the derivatives of motion variables in y-direction. More...
unsigned int	m_outputFrequency
Array< OneD, std::ofstream >	m_outputStream
std::string	m_outputFile_fce
std::string	m_outputFile_mot

Friends
class	MemoryManager< ForcingMovingBody >

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 49 of file ForcingMovingBody.h.

Constructor & Destructor Documentation

ForcingMovingBody()

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

private

Definition at line 48 of file ForcingMovingBody.cpp.

    : Forcing(pSession, pEquation)
{
}

Member Function Documentation

CheckIsFromFile()

void Nektar::ForcingMovingBody::CheckIsFromFile ( const TiXmlElement *  pForce )

private

Definition at line 1076 of file ForcingMovingBody.cpp.

{
 
    m_funcName  = Array<OneD, std::string>(3);
    m_motion    = Array<OneD, std::string>(2);
    m_motion[0] = "x";
    m_motion[1] = "y";
 
    m_IsFromFile = Array<OneD, bool>(6);
    // Loading the x-dispalcement (m_zta) and the y-displacement (m_eta)
    // Those two variables are bith functions of z and t and the may come
    // from an equation (forced vibration) or from another solver which, given
    // the aerodynamic forces at the previous step, calculates the
    // displacements.
 
    // Get the body displacement: m_zta and m_eta
    const TiXmlElement *funcNameElmt_D =
        pForce->FirstChildElement("DISPLACEMENTS");
    ASSERTL0(funcNameElmt_D,
             "MOVINGBODYFORCE tag has to specify a function name which "
             "prescribes the body displacement as d(z,t).");
 
    m_funcName[0] = funcNameElmt_D->GetText();
    ASSERTL0(m_session->DefinesFunction(m_funcName[0]),
             "Function '" + m_funcName[0] + "' not defined.");
 
    // Get the body velocity of movement: d(m_zta)/dt and d(m_eta)/dt
    const TiXmlElement *funcNameElmt_V =
        pForce->FirstChildElement("VELOCITIES");
    ASSERTL0(funcNameElmt_D,
             "MOVINGBODYFORCE tag has to specify a function name which "
             "prescribes the body velocity of movement as v(z,t).");
 
    m_funcName[1] = funcNameElmt_V->GetText();
    ASSERTL0(m_session->DefinesFunction(m_funcName[1]),
             "Function '" + m_funcName[1] + "' not defined.");
 
    // Get the body acceleration: dd(m_zta)/ddt and dd(m_eta)/ddt
    const TiXmlElement *funcNameElmt_A =
        pForce->FirstChildElement("ACCELERATIONS");
    ASSERTL0(funcNameElmt_A,
             "MOVINGBODYFORCE tag has to specify a function name which "
             "prescribes the body acceleration as a(z,t).");
 
    m_funcName[2] = funcNameElmt_A->GetText();
    ASSERTL0(m_session->DefinesFunction(m_funcName[2]),
             "Function '" + m_funcName[2] + "' not defined.");
 
    LibUtilities::FunctionType vType;
 
    // Check Displacement x
    vType = m_session->GetFunctionType(m_funcName[0], m_motion[0]);
    if (vType == LibUtilities::eFunctionTypeExpression)
    {
        m_IsFromFile[0] = false;
    }
    else if (vType == LibUtilities::eFunctionTypeFile)
    {
        m_IsFromFile[0] = true;
    }
    else
    {
        ASSERTL0(false, "The displacements in x must be specified via an "
                        "equation <E> or a file <F>");
    }
 
    // Check Displacement y
    vType = m_session->GetFunctionType(m_funcName[0], m_motion[1]);
    if (vType == LibUtilities::eFunctionTypeExpression)
    {
        m_IsFromFile[1] = false;
    }
    else if (vType == LibUtilities::eFunctionTypeFile)
    {
        m_IsFromFile[1] = true;
    }
    else
    {
        ASSERTL0(false, "The displacements in y must be specified via an "
                        "equation <E> or a file <F>");
    }
 
    // Check Velocity x
    vType = m_session->GetFunctionType(m_funcName[1], m_motion[0]);
    if (vType == LibUtilities::eFunctionTypeExpression)
    {
        m_IsFromFile[2] = false;
    }
    else if (vType == LibUtilities::eFunctionTypeFile)
    {
        m_IsFromFile[2] = true;
    }
    else
    {
        ASSERTL0(false, "The velocities in x must be specified via an equation "
                        "<E> or a file <F>");
    }
 
    // Check Velocity y
    vType = m_session->GetFunctionType(m_funcName[1], m_motion[1]);
    if (vType == LibUtilities::eFunctionTypeExpression)
    {
        m_IsFromFile[3] = false;
    }
    else if (vType == LibUtilities::eFunctionTypeFile)
    {
        m_IsFromFile[3] = true;
    }
    else
    {
        ASSERTL0(false, "The velocities in y must be specified via an equation "
                        "<E> or a file <F>");
    }
 
    // Check Acceleration x
    vType = m_session->GetFunctionType(m_funcName[2], m_motion[0]);
    if (vType == LibUtilities::eFunctionTypeExpression)
    {
        m_IsFromFile[4] = false;
    }
    else if (vType == LibUtilities::eFunctionTypeFile)
    {
        m_IsFromFile[4] = true;
    }
    else
    {
        ASSERTL0(false, "The accelerations in x must be specified via an "
                        "equation <E> or a file <F>");
    }
 
    // Check Acceleration y
    vType = m_session->GetFunctionType(m_funcName[2], m_motion[1]);
    if (vType == LibUtilities::eFunctionTypeExpression)
    {
        m_IsFromFile[5] = false;
    }
    else if (vType == LibUtilities::eFunctionTypeFile)
    {
        m_IsFromFile[5] = true;
    }
    else
    {
        ASSERTL0(false, "The accelerations in y must be specified via an "
                        "equation <E> or a file <F>");
    }
}

References ASSERTL0, Nektar::LibUtilities::eFunctionTypeExpression, Nektar::LibUtilities::eFunctionTypeFile, m_funcName, m_IsFromFile, m_motion, and Nektar::SolverUtils::Forcing::m_session.

Referenced by v_InitObject().

create()

static SolverUtils::ForcingSharedPtr Nektar::ForcingMovingBody::create ( const LibUtilities::SessionReaderSharedPtr &  pSession, const std::weak_ptr< SolverUtils::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 55 of file ForcingMovingBody.h.

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

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

EvaluateStructDynModel()

void Nektar::ForcingMovingBody::EvaluateStructDynModel ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields, Array< OneD, NekDouble > &  Hydroforces, NekDouble  time  )

private

Definition at line 229 of file ForcingMovingBody.cpp.

{
    boost::ignore_unused(time);
 
    LibUtilities::CommSharedPtr vcomm = pFields[0]->GetComm();
    size_t colrank                    = vcomm->GetColumnComm()->GetRank();
    size_t nproc                      = vcomm->GetColumnComm()->GetSize();
 
    bool homostrip;
    m_session->MatchSolverInfo("HomoStrip", "True", homostrip, false);
 
    // number of structral modes and number of strips
    size_t npts, nstrips;
    if (!homostrip) // full resolutions
    {
        npts = m_session->GetParameter("HomModesZ");
    }
    else
    {
        m_session->LoadParameter("HomStructModesZ", npts);
        m_session->LoadParameter("Strip_Z", nstrips);
 
        // ASSERTL0(nstrips < = npts,
        //     "the number of struct. modes must be larger than that of the
        //     strips.");
    }
 
    // the hydrodynamic forces
    Array<OneD, Array<OneD, NekDouble>> fces(2);
    // forces in x-direction
    fces[0] = Array<OneD, NekDouble>(npts, 0.0);
    // forces in y-direction
    fces[1] = Array<OneD, NekDouble>(npts, 0.0);
 
    // fill the force vectors
    if (colrank == 0)
    {
        if (!homostrip) // full resolutions
        {
            Vmath::Vcopy(m_np, Hydroforces, 1, fces[0], 1);
            Vmath::Vcopy(m_np, Hydroforces + m_np, 1, fces[1], 1);
        }
        else // strip modelling
        {
            fces[0][0] = Hydroforces[0];
            fces[1][0] = Hydroforces[m_np];
        }
 
        if (!homostrip) // full resolutions
        {
            Array<OneD, NekDouble> tmp(2 * m_np);
            for (size_t i = 1; i < nproc; ++i)
            {
                vcomm->GetColumnComm()->Recv(i, tmp);
                for (size_t n = 0; n < m_np; n++)
                {
                    for (size_t j = 0; j < 2; ++j)
                    {
                        fces[j][i * m_np + n] = tmp[j * m_np + n];
                    }
                }
            }
        }
        else // strip modelling
        // if the body is submerged partly, the fces are filled partly
        // by the flow induced forces
        {
            Array<OneD, NekDouble> tmp(2);
            for (size_t i = 1; i < nstrips; ++i)
            {
                vcomm->GetColumnComm()->Recv(i, tmp);
 
                for (size_t j = 0; j < 2; ++j)
                {
                    fces[j][i] = tmp[j];
                }
            }
        }
    }
    else
    {
        if (!homostrip) // full resolutions
        {
            vcomm->GetColumnComm()->Send(0, Hydroforces);
        }
        else // strip modelling
        {
            for (size_t i = 1; i < nstrips; ++i)
            {
                if (colrank == i)
                {
                    Array<OneD, NekDouble> tmp(2);
                    tmp[0] = Hydroforces[0];
                    tmp[1] = Hydroforces[m_np];
                    vcomm->GetColumnComm()->Send(0, tmp);
                }
            }
        }
    }
 
    if (colrank == 0)
    {
        // Fictitious mass method used to stablize the explicit coupling at
        // relatively lower mass ratio
        bool fictmass;
        m_session->MatchSolverInfo("FictitiousMassMethod", "True", fictmass,
                                   false);
        if (fictmass)
        {
            NekDouble fictrho, fictdamp;
            m_session->LoadParameter("FictMass", fictrho);
            m_session->LoadParameter("FictDamp", fictdamp);
 
            static NekDouble Betaq_Coeffs[2][2] = {{1.0, 0.0}, {2.0, -1.0}};
 
            // only consider second order approximation for fictitious variables
            size_t intOrder = 2;
            size_t nint     = min(m_movingBodyCalls + 1, intOrder);
            size_t nlevels  = m_fV[0].size();
 
            for (size_t i = 0; i < m_motion.size(); ++i)
            {
                RollOver(m_fV[i]);
                RollOver(m_fA[i]);
 
                size_t Voffset = npts;
                size_t Aoffset = 2 * npts;
 
                Vmath::Vcopy(npts, m_MotionVars[i] + Voffset, 1, m_fV[i][0], 1);
                Vmath::Vcopy(npts, m_MotionVars[i] + Aoffset, 1, m_fA[i][0], 1);
 
                // Extrapolate to n+1
                Vmath::Smul(npts, Betaq_Coeffs[nint - 1][nint - 1],
                            m_fV[i][nint - 1], 1, m_fV[i][nlevels - 1], 1);
                Vmath::Smul(npts, Betaq_Coeffs[nint - 1][nint - 1],
                            m_fA[i][nint - 1], 1, m_fA[i][nlevels - 1], 1);
 
                for (size_t n = 0; n < nint - 1; ++n)
                {
                    Vmath::Svtvp(npts, Betaq_Coeffs[nint - 1][n], m_fV[i][n], 1,
                                 m_fV[i][nlevels - 1], 1, m_fV[i][nlevels - 1],
                                 1);
                    Vmath::Svtvp(npts, Betaq_Coeffs[nint - 1][n], m_fA[i][n], 1,
                                 m_fA[i][nlevels - 1], 1, m_fA[i][nlevels - 1],
                                 1);
                }
 
                // Add the fictitious forces on the RHS of the equation
                Vmath::Svtvp(npts, fictdamp, m_fV[i][nlevels - 1], 1, fces[i],
                             1, fces[i], 1);
                Vmath::Svtvp(npts, fictrho, m_fA[i][nlevels - 1], 1, fces[i], 1,
                             fces[i], 1);
            }
        }
    }
    // structural solver is implemented on the root process
    if (colrank == 0)
    {
        // Tensioned cable model is evaluated in wave space
        for (size_t n = 0, cn = 1; n < m_vdim; n++, cn--)
        {
            Newmark_betaSolver(pFields, fces[cn], m_MotionVars[cn]);
        }
    }
 
    Array<OneD, NekDouble> Motvars(2 * 2 * m_np);
    // send physical coeffients to all planes of each processor
    if (!homostrip) // full resolutions
    {
        Array<OneD, NekDouble> tmp(2 * 2 * m_np);
 
        if (colrank != 0)
        {
            vcomm->GetColumnComm()->Recv(0, tmp);
            Vmath::Vcopy(2 * 2 * m_np, tmp, 1, Motvars, 1);
        }
        else
        {
            for (size_t i = 1; i < nproc; ++i)
            {
                for (size_t j = 0; j < 2; j++) // moving dimensions
                {
                    for (size_t k = 0; k < 2; k++) // disp. and vel.
                    {
                        for (size_t n = 0; n < m_np; n++)
                        {
                            tmp[j * 2 * m_np + k * m_np + n] =
                                m_MotionVars[j][k * npts + i * m_np + n];
                        }
                    }
                }
                vcomm->GetColumnComm()->Send(i, tmp);
            }
 
            for (size_t j = 0; j < 2; j++)
            {
                for (size_t k = 0; k < 2; k++)
                {
                    for (size_t n = 0; n < m_np; n++)
                    {
                        tmp[j * 2 * m_np + k * m_np + n] =
                            m_MotionVars[j][k * npts + n];
                    }
                    Vmath::Vcopy(2 * 2 * m_np, tmp, 1, Motvars, 1);
                }
            }
        }
    }
    else // strip modelling
    {
        Array<OneD, NekDouble> tmp(2 * 2);
 
        if (colrank != 0)
        {
            for (size_t j = 1; j < nproc / nstrips; j++)
            {
                if (colrank == j * nstrips)
                {
                    vcomm->GetColumnComm()->Recv(0, tmp);
 
                    for (size_t plane = 0; plane < m_np; plane++)
                    {
                        for (size_t var = 0; var < 2; var++)
                        {
                            for (size_t k = 0; k < 2; k++)
                            {
                                Motvars[var * 2 * m_np + k * m_np + plane] =
                                    tmp[var * 2 + k];
                            }
                        }
                    }
                }
            }
 
            for (size_t i = 1; i < nstrips; i++)
            {
                for (size_t j = 0; j < nproc / nstrips; j++)
                {
                    if (colrank == i + j * nstrips)
                    {
                        vcomm->GetColumnComm()->Recv(0, tmp);
 
                        for (size_t plane = 0; plane < m_np; plane++)
                        {
                            for (size_t var = 0; var < 2; var++)
                            {
                                for (size_t k = 0; k < 2; k++)
                                {
                                    Motvars[var * 2 * m_np + k * m_np + plane] =
                                        tmp[var * 2 + k];
                                }
                            }
                        }
                    }
                }
            }
        }
        else
        {
            for (size_t j = 0; j < 2; ++j)
            {
                for (size_t k = 0; k < 2; ++k)
                {
                    tmp[j * 2 + k] = m_MotionVars[j][k * npts];
                }
            }
 
            for (size_t j = 1; j < nproc / nstrips; j++)
            {
                vcomm->GetColumnComm()->Send(j * nstrips, tmp);
            }
 
            for (size_t plane = 0; plane < m_np; plane++)
            {
                for (size_t j = 0; j < 2; j++)
                {
                    for (size_t k = 0; k < 2; ++k)
                    {
                        Motvars[j * 2 * m_np + k * m_np + plane] =
                            m_MotionVars[j][k * npts];
                    }
                }
            }
 
            for (size_t i = 1; i < nstrips; ++i)
            {
                for (size_t j = 0; j < 2; ++j)
                {
                    for (size_t k = 0; k < 2; ++k)
                    {
                        tmp[j * 2 + k] = m_MotionVars[j][i + k * npts];
                    }
                }
 
                for (size_t j = 0; j < nproc / nstrips; j++)
                {
                    vcomm->GetColumnComm()->Send(i + j * nstrips, tmp);
                }
            }
        }
    }
 
    // Set the m_forcing term based on the motion of the cable
    for (size_t var = 0; var < 2; var++)
    {
        for (size_t plane = 0; plane < m_np; plane++)
        {
            size_t n = pFields[0]->GetPlane(plane)->GetTotPoints();
 
            Array<OneD, NekDouble> tmp;
 
            size_t offset  = plane * n;
            size_t xoffset = var * m_np + plane;
            size_t yoffset = 2 * m_np + xoffset;
 
            Vmath::Fill(n, Motvars[xoffset], tmp = m_zta[var] + offset, 1);
            Vmath::Fill(n, Motvars[yoffset], tmp = m_eta[var] + offset, 1);
        }
    }
}

References Vmath::Fill(), m_eta, m_fA, m_fV, m_motion, m_MotionVars, m_movingBodyCalls, m_np, Nektar::SolverUtils::Forcing::m_session, m_vdim, m_zta, Newmark_betaSolver(), RollOver(), Vmath::Smul(), Vmath::Svtvp(), and Vmath::Vcopy().

Referenced by v_Apply().

InitialiseCableModel()

void Nektar::ForcingMovingBody::InitialiseCableModel ( const LibUtilities::SessionReaderSharedPtr &  pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields  )

private

Definition at line 685 of file ForcingMovingBody.cpp.

{
    boost::ignore_unused(pSession, pFields);
 
    m_movingBodyCalls = 0;
    m_session->LoadParameter("Kinvis", m_kinvis);
    m_session->LoadParameter("TimeStep", m_timestep, 0.01);
 
    LibUtilities::CommSharedPtr vcomm = pFields[0]->GetComm();
    size_t colrank                    = vcomm->GetColumnComm()->GetRank();
    size_t nproc                      = vcomm->GetColumnComm()->GetSize();
 
    // number of structral modes
    size_t npts;
    bool homostrip;
    m_session->MatchSolverInfo("HomoStrip", "True", homostrip, false);
 
    if (!homostrip) // full resolutions
    {
        npts = m_session->GetParameter("HomModesZ");
    }
    else
    {
        m_session->LoadParameter("HomStructModesZ", npts);
    }
 
    m_MotionVars    = Array<OneD, Array<OneD, NekDouble>>(2);
    m_MotionVars[0] = Array<OneD, NekDouble>(3 * npts, 0.0);
    m_MotionVars[1] = Array<OneD, NekDouble>(3 * npts, 0.0);
 
    Array<OneD, unsigned int> ZIDs;
    ZIDs = pFields[0]->GetZIDs();
    m_np = ZIDs.size();
 
    std::string vibtype = m_session->GetSolverInfo("VibrationType");
 
    if (boost::iequals(vibtype, "Constrained"))
    {
        m_vdim = 1;
    }
    else if (boost::iequals(vibtype, "Free"))
    {
        m_vdim = 2;
    }
    else if (boost::iequals(vibtype, "Forced"))
    {
        m_vdim = 0;
        return;
    }
 
    if (!homostrip)
    {
        m_session->LoadParameter("LZ", m_lhom);
        int nplanes = m_session->GetParameter("HomModesZ");
        m_FFT = LibUtilities::GetNektarFFTFactory().CreateInstance("NekFFTW",
                                                                   nplanes);
    }
    else
    {
        int nstrips;
        NekDouble DistStrip;
 
        m_session->LoadParameter("DistStrip", DistStrip);
        m_session->LoadParameter("Strip_Z", nstrips);
        m_lhom = nstrips * DistStrip;
        m_FFT  = LibUtilities::GetNektarFFTFactory().CreateInstance("NekFFTW",
                                                                   nstrips);
    }
 
    // load the structural dynamic parameters from xml file
    m_session->LoadParameter("StructRho", m_structrho);
    m_session->LoadParameter("StructDamp", m_structdamp, 0.0);
 
    // Identify whether the fictitious mass method is active for explicit
    // coupling of fluid solver and structural dynamics solver
    bool fictmass;
    m_session->MatchSolverInfo("FictitiousMassMethod", "True", fictmass, false);
    if (fictmass)
    {
        NekDouble fictrho, fictdamp;
        m_session->LoadParameter("FictMass", fictrho);
        m_session->LoadParameter("FictDamp", fictdamp);
        m_structrho += fictrho;
        m_structdamp += fictdamp;
 
        // Storage array of Struct Velocity and Acceleration used for
        // extrapolation of fictitious force
        m_fV = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(2);
        m_fA = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(2);
        for (size_t i = 0; i < m_motion.size(); ++i)
        {
            m_fV[i] = Array<OneD, Array<OneD, NekDouble>>(2);
            m_fA[i] = Array<OneD, Array<OneD, NekDouble>>(2);
 
            for (size_t n = 0; n < 2; ++n)
            {
                m_fV[i][n] = Array<OneD, NekDouble>(npts, 0.0);
                m_fA[i][n] = Array<OneD, NekDouble>(npts, 0.0);
            }
        }
    }
 
    // Setting the coefficient matrices for solving structural dynamic ODEs
    SetDynEqCoeffMatrix(pFields);
 
    // Set initial condition for cable's motion
    size_t cnt = 0;
 
    for (size_t j = 0; j < m_funcName.size(); j++)
    {
        // loading from the specified files through inputstream
        if (m_IsFromFile[cnt] && m_IsFromFile[cnt + 1])
        {
            std::ifstream inputStream;
            size_t nzpoints;
 
            if (homostrip)
            {
                m_session->LoadParameter("HomStructModesZ", nzpoints);
            }
            else
            {
                nzpoints = pFields[0]->GetHomogeneousBasis()->GetNumModes();
            }
 
            if (vcomm->GetRank() == 0)
            {
                std::string filename =
                    m_session->GetFunctionFilename(m_funcName[j], m_motion[0]);
 
                // Open intputstream for cable motions
                inputStream.open(filename.c_str());
 
                // Import the head string from the file
                Array<OneD, std::string> tmp(9);
                for (size_t n = 0; n < tmp.size(); n++)
                {
                    inputStream >> tmp[n];
                }
 
                NekDouble time, z_cds;
                // Import the motion variables from the file
                for (size_t n = 0; n < nzpoints; n++)
                {
                    inputStream >> setprecision(6) >> time;
                    inputStream >> setprecision(6) >> z_cds;
                    inputStream >> setprecision(8) >> m_MotionVars[0][n];
                    inputStream >> setprecision(8) >>
                        m_MotionVars[0][n + nzpoints];
                    inputStream >> setprecision(8) >>
                        m_MotionVars[0][n + 2 * nzpoints];
                    inputStream >> setprecision(8) >> m_MotionVars[1][n];
                    inputStream >> setprecision(8) >>
                        m_MotionVars[1][n + nzpoints];
                    inputStream >> setprecision(8) >>
                        m_MotionVars[1][n + 2 * nzpoints];
                }
                // Close inputstream for cable motions
                inputStream.close();
            }
            cnt = cnt + 2;
        }
        else // Evaluate from the functions specified in xml file
        {
            if (!homostrip)
            {
                size_t nzpoints =
                    pFields[0]->GetHomogeneousBasis()->GetNumModes();
                Array<OneD, NekDouble> z_coords(nzpoints, 0.0);
                Array<OneD, const NekDouble> pts =
                    pFields[0]->GetHomogeneousBasis()->GetZ();
 
                Vmath::Smul(nzpoints, m_lhom / 2.0, pts, 1, z_coords, 1);
                Vmath::Sadd(nzpoints, m_lhom / 2.0, z_coords, 1, z_coords, 1);
 
                Array<OneD, NekDouble> x0(m_np, 0.0);
                Array<OneD, NekDouble> x1(m_np, 0.0);
                Array<OneD, NekDouble> x2(m_np, 0.0);
                for (size_t plane = 0; plane < m_np; plane++)
                {
                    x2[plane] = z_coords[ZIDs[plane]];
                }
 
                Array<OneD, NekDouble> tmp(m_np, 0.0);
                Array<OneD, NekDouble> tmp0(m_np, 0.0);
                Array<OneD, NekDouble> tmp1(m_np, 0.0);
                LibUtilities::EquationSharedPtr ffunc0, ffunc1;
 
                ffunc0 = m_session->GetFunction(m_funcName[j], m_motion[0]);
                ffunc1 = m_session->GetFunction(m_funcName[j], m_motion[1]);
                ffunc0->Evaluate(x0, x1, x2, 0.0, tmp0);
                ffunc1->Evaluate(x0, x1, x2, 0.0, tmp1);
 
                size_t offset = j * npts;
                Vmath::Vcopy(m_np, tmp0, 1, tmp = m_MotionVars[0] + offset, 1);
                Vmath::Vcopy(m_np, tmp1, 1, tmp = m_MotionVars[1] + offset, 1);
 
                if (colrank == 0)
                {
                    for (size_t i = 1; i < nproc; ++i)
                    {
                        vcomm->GetColumnComm()->Recv(i, tmp0);
                        vcomm->GetColumnComm()->Recv(i, tmp1);
                        Vmath::Vcopy(m_np, tmp0, 1,
                                     tmp = m_MotionVars[0] + offset + i * m_np,
                                     1);
                        Vmath::Vcopy(m_np, tmp1, 1,
                                     tmp = m_MotionVars[1] + offset + i * m_np,
                                     1);
                    }
                }
                else
                {
                    vcomm->GetColumnComm()->Send(0, tmp0);
                    vcomm->GetColumnComm()->Send(0, tmp1);
                }
            }
            else
            {
                if (colrank == 0)
                {
                    size_t nstrips;
                    m_session->LoadParameter("Strip_Z", nstrips);
 
                    ASSERTL0(
                        m_session->DefinesSolverInfo("USEFFT"),
                        "Fourier transformation of cable motion is currently "
                        "implemented only for FFTW module.");
 
                    NekDouble DistStrip;
                    m_session->LoadParameter("DistStrip", DistStrip);
 
                    Array<OneD, NekDouble> x0(npts, 0.0);
                    Array<OneD, NekDouble> x1(npts, 0.0);
                    Array<OneD, NekDouble> x2(npts, 0.0);
                    Array<OneD, NekDouble> tmp(npts, 0.0);
                    Array<OneD, NekDouble> tmp0(npts, 0.0);
                    Array<OneD, NekDouble> tmp1(npts, 0.0);
                    for (size_t plane = 0; plane < npts; plane++)
                    {
                        x2[plane] = plane * DistStrip;
                    }
                    LibUtilities::EquationSharedPtr ffunc0, ffunc1;
                    ffunc0 = m_session->GetFunction(m_funcName[j], m_motion[0]);
                    ffunc1 = m_session->GetFunction(m_funcName[j], m_motion[1]);
                    ffunc0->Evaluate(x0, x1, x2, 0.0, tmp0);
                    ffunc1->Evaluate(x0, x1, x2, 0.0, tmp1);
 
                    size_t offset = j * npts;
                    Vmath::Vcopy(npts, tmp0, 1, tmp = m_MotionVars[0] + offset,
                                 1);
                    Vmath::Vcopy(npts, tmp1, 1, tmp = m_MotionVars[1] + offset,
                                 1);
                }
            }
 
            cnt = cnt + 2;
        }
    }
}

References ASSERTL0, Nektar::LibUtilities::NekFactory< tKey, tBase, tParam >::CreateInstance(), Nektar::LibUtilities::GetNektarFFTFactory(), m_fA, m_FFT, m_funcName, m_fV, m_IsFromFile, m_kinvis, m_lhom, m_motion, m_MotionVars, m_movingBodyCalls, m_np, Nektar::SolverUtils::Forcing::m_session, m_structdamp, m_structrho, m_timestep, m_vdim, Vmath::Sadd(), SetDynEqCoeffMatrix(), Vmath::Smul(), and Vmath::Vcopy().

Referenced by v_InitObject().

InitialiseFilter()

void Nektar::ForcingMovingBody::InitialiseFilter ( const LibUtilities::SessionReaderSharedPtr &  pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields, const TiXmlElement *  pForce  )

private

Definition at line 1226 of file ForcingMovingBody.cpp.

{
    // Get the outputfile name, output frequency and
    // the boundary's ID for the cable's wall
    std::string typeStr = pForce->Attribute("TYPE");
    std::map<std::string, std::string> vParams;
 
    const TiXmlElement *param = pForce->FirstChildElement("PARAM");
    while (param)
    {
        ASSERTL0(param->Attribute("NAME"),
                 "Missing attribute 'NAME' for parameter in filter " + typeStr +
                     "'.");
        std::string nameStr = param->Attribute("NAME");
 
        ASSERTL0(param->GetText(), "Empty value string for param.");
        std::string valueStr = param->GetText();
 
        vParams[nameStr] = valueStr;
 
        param = param->NextSiblingElement("PARAM");
    }
 
    // Creat the filter for MovingBody, where we performed the calculation of
    // fluid forces and write both the aerodynamic forces and motion variables
    // into the output files
    m_MovBodyfilter = MemoryManager<FilterMovingBody>::AllocateSharedPtr(
        pSession, m_equ, vParams);
 
    // Initialise the object of MovingBody filter
    m_MovBodyfilter->Initialise(pFields, 0.0);
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::SolverUtils::Forcing::m_equ, and m_MovBodyfilter.

Referenced by v_InitObject().

Newmark_betaSolver()

void Nektar::ForcingMovingBody::Newmark_betaSolver ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields, Array< OneD, NekDouble > &  FcePhysinArray, Array< OneD, NekDouble > &  MotPhysinArray  )

private

Definition at line 555 of file ForcingMovingBody.cpp.

{
    boost::ignore_unused(pFields);
 
    std::string supptype = m_session->GetSolverInfo("SupportType");
 
    size_t npts = HydroForces.size();
 
    Array<OneD, Array<OneD, NekDouble>> fft_i(4);
    Array<OneD, Array<OneD, NekDouble>> fft_o(4);
 
    for (size_t i = 0; i < 4; i++)
    {
        fft_i[i] = Array<OneD, NekDouble>(npts, 0.0);
        fft_o[i] = Array<OneD, NekDouble>(npts, 0.0);
    }
 
    Vmath::Vcopy(npts, HydroForces, 1, fft_i[0], 1);
    for (size_t i = 0; i < 3; i++)
    {
        Vmath::Vcopy(npts, BodyMotions + i * npts, 1, fft_i[i + 1], 1);
    }
 
    // Implement Fourier transformation of the motion variables
    if (boost::iequals(supptype, "Free-Free"))
    {
        for (size_t j = 0; j < 4; ++j)
        {
            m_FFT->FFTFwdTrans(fft_i[j], fft_o[j]);
        }
    }
    else if (boost::iequals(supptype, "Pinned-Pinned"))
    {
        // TODO:
        size_t N = fft_i[0].size();
 
        for (size_t var = 0; var < 4; var++)
        {
            for (size_t i = 0; i < N; i++)
            {
                fft_o[var][i] = 0;
 
                for (size_t k = 0; k < N; k++)
                {
                    fft_o[var][i] +=
                        fft_i[var][k] * sin(M_PI / (N) * (k + 1 / 2) * (i + 1));
                }
            }
        }
    }
    else
    {
        ASSERTL0(false, "Unrecognized support type for cable's motion");
    }
 
    // solve the ODE in the wave space
    for (size_t i = 0; i < npts; ++i)
    {
        size_t nrows = 3;
 
        Array<OneD, NekDouble> tmp0, tmp1, tmp2;
        tmp0 = Array<OneD, NekDouble>(3, 0.0);
        tmp1 = Array<OneD, NekDouble>(3, 0.0);
        tmp2 = Array<OneD, NekDouble>(3, 0.0);
 
        for (size_t var = 0; var < 3; var++)
        {
            tmp0[var] = fft_o[var + 1][i];
        }
 
        tmp2[0] = fft_o[0][i];
 
        Blas::Dgemv('N', nrows, nrows, 1.0, &(m_CoeffMat_B[i]->GetPtr())[0],
                    nrows, &tmp0[0], 1, 0.0, &tmp1[0], 1);
        Blas::Dgemv('N', nrows, nrows, 1.0 / m_structrho,
                    &(m_CoeffMat_A[i]->GetPtr())[0], nrows, &tmp2[0], 1, 1.0,
                    &tmp1[0], 1);
 
        for (size_t var = 0; var < 3; var++)
        {
            fft_i[var][i] = tmp1[var];
        }
    }
 
    // get physical coeffients via Backward fourier transformation of wave
    // coefficients
    if (boost::iequals(supptype, "Free-Free"))
    {
        for (size_t var = 0; var < 3; var++)
        {
            m_FFT->FFTBwdTrans(fft_i[var], fft_o[var]);
        }
    }
    else if (boost::iequals(supptype, "Pinned-Pinned"))
    {
        // TODO:
        size_t N = fft_i[0].size();
 
        for (size_t var = 0; var < 3; var++)
        {
            for (size_t i = 0; i < N; i++)
            {
                fft_o[var][i] = 0;
 
                for (size_t k = 0; k < N; k++)
                {
                    fft_o[var][i] += fft_i[var][k] *
                                     sin(M_PI / (N) * (k + 1) * (i + 1 / 2)) *
                                     2 / N;
                }
            }
        }
    }
    else
    {
        ASSERTL0(false, "Unrecognized support type for cable's motion");
    }
 
    for (size_t i = 0; i < 3; i++)
    {
        Array<OneD, NekDouble> tmp(npts, 0.0);
        Vmath::Vcopy(npts, fft_o[i], 1, tmp = BodyMotions + i * npts, 1);
    }
}

References ASSERTL0, Blas::Dgemv(), m_CoeffMat_A, m_CoeffMat_B, m_FFT, Nektar::SolverUtils::Forcing::m_session, m_structrho, and Vmath::Vcopy().

Referenced by EvaluateStructDynModel().

RollOver()

void Nektar::ForcingMovingBody::RollOver ( Array< OneD, Array< OneD, NekDouble > > &  input )

private

Function to roll time-level storages to the next step layout. The stored data associated with the oldest time-level (not required anymore) are moved to the top, where they will be overwritten as the solution process progresses.

Definition at line 1058 of file ForcingMovingBody.cpp.

{
    size_t nlevels = input.size();
 
    Array<OneD, NekDouble> tmp;
    tmp = input[nlevels - 1];
 
    for (size_t n = nlevels - 1; n > 0; --n)
    {
        input[n] = input[n - 1];
    }
 
    input[0] = tmp;
}

Referenced by EvaluateStructDynModel().

SetDynEqCoeffMatrix()

void Nektar::ForcingMovingBody::SetDynEqCoeffMatrix ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields )

private

Definition at line 951 of file ForcingMovingBody.cpp.

{
    boost::ignore_unused(pFields);
 
    size_t nplanes;
 
    bool homostrip;
    m_session->MatchSolverInfo("HomoStrip", "True", homostrip, false);
 
    if (!homostrip)
    {
        nplanes = m_session->GetParameter("HomModesZ");
    }
    else
    {
        m_session->LoadParameter("Strip_Z", nplanes);
    }
 
    m_CoeffMat_A = Array<OneD, DNekMatSharedPtr>(nplanes);
    m_CoeffMat_B = Array<OneD, DNekMatSharedPtr>(nplanes);
 
    NekDouble tmp1, tmp2, tmp3;
    NekDouble tmp4, tmp5, tmp6, tmp7;
 
    // load the structural dynamic parameters from xml file
    NekDouble cabletension;
    NekDouble bendingstiff;
    NekDouble structstiff;
    m_session->LoadParameter("StructStiff", structstiff, 0.0);
    m_session->LoadParameter("CableTension", cabletension, 0.0);
    m_session->LoadParameter("BendingStiff", bendingstiff, 0.0);
 
    tmp1 = m_timestep * m_timestep;
    tmp2 = structstiff / m_structrho;
    tmp3 = m_structdamp / m_structrho;
    tmp4 = cabletension / m_structrho;
    tmp5 = bendingstiff / m_structrho;
 
    // solve the ODE in the wave space for cable motion to obtain disp, vel and
    // accel
 
    std::string supptype = m_session->GetSolverInfo("SupportType");
 
    for (size_t plane = 0; plane < nplanes; plane++)
    {
        int nel = 3;
        m_CoeffMat_A[plane] =
            MemoryManager<DNekMat>::AllocateSharedPtr(nel, nel);
        m_CoeffMat_B[plane] =
            MemoryManager<DNekMat>::AllocateSharedPtr(nel, nel);
 
        // Initialised to avoid compiler warnings.
        unsigned int K = 0;
        NekDouble beta = 0.0;
 
        if (boost::iequals(supptype, "Free-Free"))
        {
            K    = plane / 2;
            beta = 2.0 * M_PI / m_lhom;
        }
        else if (boost::iequals(supptype, "Pinned-Pinned"))
        {
            K    = plane + 1;
            beta = M_PI / m_lhom;
        }
        else
        {
            ASSERTL0(false, "Unrecognized support type for cable's motion");
        }
 
        tmp6 = beta * K;
        tmp6 = tmp6 * tmp6;
        tmp7 = tmp6 * tmp6;
        tmp7 = tmp2 + tmp4 * tmp6 + tmp5 * tmp7;
 
        (*m_CoeffMat_A[plane])(0, 0) = tmp7;
        (*m_CoeffMat_A[plane])(0, 1) = tmp3;
        (*m_CoeffMat_A[plane])(0, 2) = 1.0;
        (*m_CoeffMat_A[plane])(1, 0) = 0.0;
        (*m_CoeffMat_A[plane])(1, 1) = 1.0;
        (*m_CoeffMat_A[plane])(1, 2) = -m_timestep / 2.0;
        (*m_CoeffMat_A[plane])(2, 0) = 1.0;
        (*m_CoeffMat_A[plane])(2, 1) = 0.0;
        (*m_CoeffMat_A[plane])(2, 2) = -tmp1 / 4.0;
        (*m_CoeffMat_B[plane])(0, 0) = 0.0;
        (*m_CoeffMat_B[plane])(0, 1) = 0.0;
        (*m_CoeffMat_B[plane])(0, 2) = 0.0;
        (*m_CoeffMat_B[plane])(1, 0) = 0.0;
        (*m_CoeffMat_B[plane])(1, 1) = 1.0;
        (*m_CoeffMat_B[plane])(1, 2) = m_timestep / 2.0;
        (*m_CoeffMat_B[plane])(2, 0) = 1.0;
        (*m_CoeffMat_B[plane])(2, 1) = m_timestep;
        (*m_CoeffMat_B[plane])(2, 2) = tmp1 / 4.0;
 
        m_CoeffMat_A[plane]->Invert();
        (*m_CoeffMat_B[plane]) =
            (*m_CoeffMat_A[plane]) * (*m_CoeffMat_B[plane]);
    }
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::LibUtilities::beta, m_CoeffMat_A, m_CoeffMat_B, m_lhom, Nektar::SolverUtils::Forcing::m_session, m_structdamp, m_structrho, and m_timestep.

Referenced by InitialiseCableModel().

v_Apply()

void Nektar::ForcingMovingBody::v_Apply ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields, const Array< OneD, Array< OneD, NekDouble > > &  inarray, Array< OneD, Array< OneD, NekDouble > > &  outarray, const NekDouble &  time  )

overrideprotectedvirtual

Implements Nektar::SolverUtils::Forcing.

Definition at line 104 of file ForcingMovingBody.cpp.

{
    boost::ignore_unused(inarray, outarray);
 
    // Update the forces from the calculation of fluid field, which is
    // implemented in the movingbody filter
    Array<OneD, NekDouble> Hydroforces(2 * m_np, 0.0);
    m_MovBodyfilter->UpdateForce(m_session, pFields, Hydroforces, time);
 
    // for "free" type (m_vdim = 2), the cable vibrates both in streamwise and
    // crossflow dimections, for "constrained" type (m_vdim = 1), the cable only
    // vibrates in crossflow dimection, and for "forced" type (m_vdim = 0), the
    // calbe vibrates specifically along a given function or file
    if (m_vdim == 1 || m_vdim == 2)
    {
        // For free vibration case, displacements, velocities and acceleartions
        // are obtained through solving structure dynamic model
        EvaluateStructDynModel(pFields, Hydroforces, time);
 
        // Convert result to format required by mapping
        size_t physTot = pFields[0]->GetTotPoints();
        Array<OneD, Array<OneD, NekDouble>> coords(3);
        Array<OneD, Array<OneD, NekDouble>> coordsVel(3);
        for (size_t i = 0; i < 3; i++)
        {
            coords[i]    = Array<OneD, NekDouble>(physTot, 0.0);
            coordsVel[i] = Array<OneD, NekDouble>(physTot, 0.0);
        }
        // Get original coordinates
        pFields[0]->GetCoords(coords[0], coords[1], coords[2]);
 
        // Add displacement to coordinates
        Vmath::Vadd(physTot, coords[0], 1, m_zta[0], 1, coords[0], 1);
        Vmath::Vadd(physTot, coords[1], 1, m_eta[0], 1, coords[1], 1);
        // Copy velocities
        Vmath::Vcopy(physTot, m_zta[1], 1, coordsVel[0], 1);
        Vmath::Vcopy(physTot, m_eta[1], 1, coordsVel[1], 1);
 
        // Update mapping
        m_mapping->UpdateMapping(time, coords, coordsVel);
    }
    else if (m_vdim == 0)
    {
        // For forced vibration case, load from specified file or function
        size_t cnt = 0;
        for (size_t j = 0; j < m_funcName.size(); j++)
        {
            if (m_IsFromFile[cnt] && m_IsFromFile[cnt + 1])
            {
                ASSERTL0(false, "Motion loading from file needs specific "
                                "implementation: Work in Progress!");
            }
            else
            {
                GetFunction(pFields, m_session, m_funcName[j], true)
                    ->Evaluate(m_motion[0], m_zta[j], time);
                GetFunction(pFields, m_session, m_funcName[j], true)
                    ->Evaluate(m_motion[1], m_eta[j], time);
                cnt = cnt + 2;
            }
        }
 
        // Update mapping
        m_mapping->UpdateMapping(time);
 
        // Convert result from mapping
        size_t physTot = pFields[0]->GetTotPoints();
        Array<OneD, Array<OneD, NekDouble>> coords(3);
        Array<OneD, Array<OneD, NekDouble>> coordsVel(3);
        for (size_t i = 0; i < 3; i++)
        {
            coords[i]    = Array<OneD, NekDouble>(physTot, 0.0);
            coordsVel[i] = Array<OneD, NekDouble>(physTot, 0.0);
        }
        // Get original coordinates
        pFields[0]->GetCoords(coords[0], coords[1], coords[2]);
 
        // Get Coordinates and coord velocity from mapping
        m_mapping->GetCartesianCoordinates(m_zta[0], m_eta[0], coords[2]);
        m_mapping->GetCoordVelocity(coordsVel);
 
        // Calculate displacement
        Vmath::Vsub(physTot, m_zta[0], 1, coords[0], 1, m_zta[0], 1);
        Vmath::Vsub(physTot, m_eta[0], 1, coords[1], 1, m_eta[0], 1);
 
        Vmath::Vcopy(physTot, coordsVel[0], 1, m_zta[1], 1);
        Vmath::Vcopy(physTot, coordsVel[1], 1, m_eta[1], 1);
 
        for (size_t var = 0; var < 3; var++)
        {
            for (size_t plane = 0; plane < m_np; plane++)
            {
                size_t n      = pFields[0]->GetPlane(plane)->GetTotPoints();
                size_t offset = plane * n;
                size_t Offset = var * m_np + plane;
 
                m_MotionVars[0][Offset] = m_zta[var][offset];
                m_MotionVars[1][Offset] = m_eta[var][offset];
            }
        }
    }
    else
    {
        ASSERTL0(false, "Unrecogized vibration type for cable's dynamic model");
    }
 
    LibUtilities::CommSharedPtr vcomm = pFields[0]->GetComm();
    size_t colrank                    = vcomm->GetColumnComm()->GetRank();
    // Pass the variables of the cable's motion to the movingbody filter
    if (colrank == 0)
    {
        size_t n = m_MotionVars[0].size();
        Array<OneD, NekDouble> tmpArray(2 * n), tmp(n);
        Vmath::Vcopy(n, m_MotionVars[0], 1, tmpArray, 1);
        Vmath::Vcopy(n, m_MotionVars[1], 1, tmp = tmpArray + n, 1);
        m_MovBodyfilter->UpdateMotion(m_session, pFields, tmpArray, time);
    }
}

References ASSERTL0, EvaluateStructDynModel(), Nektar::SolverUtils::Forcing::GetFunction(), m_eta, m_funcName, m_IsFromFile, m_mapping, m_motion, m_MotionVars, m_MovBodyfilter, m_np, Nektar::SolverUtils::Forcing::m_session, m_vdim, m_zta, Vmath::Vadd(), Vmath::Vcopy(), and Vmath::Vsub().

v_InitObject()

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

overrideprotectedvirtual

Implements Nektar::SolverUtils::Forcing.

Definition at line 55 of file ForcingMovingBody.cpp.

{
    boost::ignore_unused(pNumForcingFields);
 
    // Just 3D homogenous 1D problems can use this techinque
    ASSERTL0(pFields[0]->GetExpType() == MultiRegions::e3DH1D,
             "Moving body implemented just for 3D Homogenous 1D expansions.");
 
    // At this point we know in the xml file where those quantities
    // are declared (equation or file) - via a function name which is now
    // stored in funcNameD etc. We need now to fill in with this info the
    // m_zta and m_eta vectors (actuallythey are matrices) Array to control if
    // the motion is determined by an equation or is from a file.(not Nektar++)
    // check if we need to load a file or we have an equation
    CheckIsFromFile(pForce);
 
    // Initialise movingbody filter
    InitialiseFilter(m_session, pFields, pForce);
 
    // Initialise the cable model
    InitialiseCableModel(m_session, pFields);
 
    // Load mapping
    m_mapping = GlobalMapping::Mapping::Load(m_session, pFields);
    m_mapping->SetTimeDependent(true);
 
    if (m_vdim > 0)
    {
        m_mapping->SetFromFunction(false);
    }
 
    m_zta = Array<OneD, Array<OneD, NekDouble>>(3);
    m_eta = Array<OneD, Array<OneD, NekDouble>>(3);
    // What are this bi-dimensional vectors
    // ------------------------------------------ m_zta[0] = m_zta |  m_eta[0] =
    // m_eta                      | m_zta[1] = d(m_zta)/dt               |
    // m_eta[1] = d(m_eta)/dt                | m_zta[2] = dd(m_zta)/ddtt |
    // m_eta[2] = dd(m_eta)/ddtt             |
    //--------------------------------------------------------------------------------
    size_t phystot = pFields[0]->GetTotPoints();
    for (size_t i = 0; i < m_zta.size(); i++)
    {
        m_zta[i] = Array<OneD, NekDouble>(phystot, 0.0);
        m_eta[i] = Array<OneD, NekDouble>(phystot, 0.0);
    }
}

References ASSERTL0, CheckIsFromFile(), Nektar::MultiRegions::e3DH1D, InitialiseCableModel(), InitialiseFilter(), Nektar::GlobalMapping::Mapping::Load(), m_eta, m_mapping, Nektar::SolverUtils::Forcing::m_session, m_vdim, and m_zta.

Friends And Related Function Documentation

MemoryManager< ForcingMovingBody >

friend class MemoryManager< ForcingMovingBody >

friend

Definition at line 1 of file ForcingMovingBody.h.

Member Data Documentation

className

std::string Nektar::ForcingMovingBody::className

static

Initial value:

=
    SolverUtils::GetForcingFactory().RegisterCreatorFunction(
        "MovingBody", ForcingMovingBody::create, "Moving Body Forcing")

Name of the class.

Definition at line 69 of file ForcingMovingBody.h.

m_Aeroforces

Array<OneD, NekDouble> Nektar::ForcingMovingBody::m_Aeroforces

private

storage for the cable's force(x,y) variables

Definition at line 130 of file ForcingMovingBody.h.

m_CoeffMat_A

Array<OneD, DNekMatSharedPtr> Nektar::ForcingMovingBody::m_CoeffMat_A

private

matrices in Newmart-beta method

Definition at line 138 of file ForcingMovingBody.h.

Referenced by Newmark_betaSolver(), and SetDynEqCoeffMatrix().

m_CoeffMat_B

Array<OneD, DNekMatSharedPtr> Nektar::ForcingMovingBody::m_CoeffMat_B

private

matrices in Newmart-beta method

Definition at line 140 of file ForcingMovingBody.h.

Referenced by Newmark_betaSolver(), and SetDynEqCoeffMatrix().

m_eta

Array<OneD, Array<OneD, NekDouble> > Nektar::ForcingMovingBody::m_eta

private

Store the derivatives of motion variables in y-direction.

Definition at line 150 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), v_Apply(), and v_InitObject().

m_fA

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::ForcingMovingBody::m_fA

private

fictitious acceleration storage

Definition at line 136 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), and InitialiseCableModel().

m_FFT

LibUtilities::NektarFFTSharedPtr Nektar::ForcingMovingBody::m_FFT

private

Definition at line 126 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel(), and Newmark_betaSolver().

m_funcName

Array<OneD, std::string> Nektar::ForcingMovingBody::m_funcName

private

[0] is displacements, [1] is velocities, [2] is accelerations

Definition at line 142 of file ForcingMovingBody.h.

Referenced by CheckIsFromFile(), InitialiseCableModel(), and v_Apply().

m_fV

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::ForcingMovingBody::m_fV

private

fictitious velocity storage

Definition at line 134 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), and InitialiseCableModel().

m_IsFromFile

Array<OneD, bool> Nektar::ForcingMovingBody::m_IsFromFile

private

do determine if the the body motion come from an extern file

Definition at line 146 of file ForcingMovingBody.h.

Referenced by CheckIsFromFile(), InitialiseCableModel(), and v_Apply().

m_kinvis

NekDouble Nektar::ForcingMovingBody::m_kinvis

private

fluid viscous

Definition at line 123 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel().

m_lhom

NekDouble Nektar::ForcingMovingBody::m_lhom

private

length ratio of the cable

Definition at line 122 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel(), and SetDynEqCoeffMatrix().

m_mapping

GlobalMapping::MappingSharedPtr Nektar::ForcingMovingBody::m_mapping

protected

Definition at line 73 of file ForcingMovingBody.h.

Referenced by v_Apply(), and v_InitObject().

m_motion

Array<OneD, std::string> Nektar::ForcingMovingBody::m_motion

private

motion direction: [0] is 'x' and [1] is 'y'

Definition at line 144 of file ForcingMovingBody.h.

Referenced by CheckIsFromFile(), EvaluateStructDynModel(), InitialiseCableModel(), and v_Apply().

m_MotionVars

Array<OneD, Array<OneD, NekDouble> > Nektar::ForcingMovingBody::m_MotionVars

private

storage for the cable's motion(x,y) variables

Definition at line 132 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), InitialiseCableModel(), and v_Apply().

m_MovBodyfilter

FilterMovingBodySharedPtr Nektar::ForcingMovingBody::m_MovBodyfilter

private

Definition at line 128 of file ForcingMovingBody.h.

Referenced by InitialiseFilter(), and v_Apply().

m_movingBodyCalls

size_t Nektar::ForcingMovingBody::m_movingBodyCalls

private

number of times the movbody have been called

Definition at line 116 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), and InitialiseCableModel().

m_np

size_t Nektar::ForcingMovingBody::m_np

private

number of planes per processors

Definition at line 117 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), InitialiseCableModel(), and v_Apply().

m_outputFile_fce

std::string Nektar::ForcingMovingBody::m_outputFile_fce

private

Definition at line 154 of file ForcingMovingBody.h.

m_outputFile_mot

std::string Nektar::ForcingMovingBody::m_outputFile_mot

private

Definition at line 155 of file ForcingMovingBody.h.

m_outputFrequency

unsigned int Nektar::ForcingMovingBody::m_outputFrequency

private

Definition at line 152 of file ForcingMovingBody.h.

m_outputStream

Array<OneD, std::ofstream> Nektar::ForcingMovingBody::m_outputStream

private

Definition at line 153 of file ForcingMovingBody.h.

m_structdamp

NekDouble Nektar::ForcingMovingBody::m_structdamp

private

damping ratio of the cable

Definition at line 121 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel(), and SetDynEqCoeffMatrix().

m_structrho

NekDouble Nektar::ForcingMovingBody::m_structrho

private

mass of the cable per unit length

Definition at line 120 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel(), Newmark_betaSolver(), and SetDynEqCoeffMatrix().

m_timestep

NekDouble Nektar::ForcingMovingBody::m_timestep

private

time step

Definition at line 124 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel(), and SetDynEqCoeffMatrix().

m_vdim

size_t Nektar::ForcingMovingBody::m_vdim

private

vibration dimension

Definition at line 118 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), InitialiseCableModel(), v_Apply(), and v_InitObject().

m_zta

Array<OneD, Array<OneD, NekDouble> > Nektar::ForcingMovingBody::m_zta

private

Store the derivatives of motion variables in x-direction.

Definition at line 148 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), v_Apply(), and v_InitObject().