Extrapolate.h

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: Extrapolate.h
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Abstract base class for Extrapolate.
//
///////////////////////////////////////////////////////////////////////////////
 
#ifndef NEKTAR_SOLVERS_EXTRAPOLATE_H
#define NEKTAR_SOLVERS_EXTRAPOLATE_H
 
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Memory/NekMemoryManager.hpp>
#include <LibUtilities/TimeIntegration/TimeIntegrationScheme.h>
#include <MultiRegions/ExpList.h>
#include <SolverUtils/AdvectionSystem.h>
#include <SolverUtils/Forcing/Forcing.h>
 
namespace Nektar
{
 
enum HBCType
{
    eNOHBC,
    eHBCNeumann, // Standard High Order BC
    eOBC, // High Order outflow BC (Neumann-Dirichlet) from Dong et al JCP 2014
    eConvectiveOBC // Convective High Order (Robin type) BC from Dong JCP 2015
};
 
// Forward declaration
class Extrapolate;
typedef std::shared_ptr<Extrapolate> ExtrapolateSharedPtr;
typedef LibUtilities::NekFactory<
    std::string, Extrapolate, const LibUtilities::SessionReaderSharedPtr &,
    Array<OneD, MultiRegions::ExpListSharedPtr> &,
    MultiRegions::ExpListSharedPtr &, const Array<OneD, int> &,
    const SolverUtils::AdvectionSharedPtr &>
    ExtrapolateFactory;
 
struct HighOrderOutflow;
typedef std::shared_ptr<HighOrderOutflow> HighOrderOutflowSharedPtr;
 
ExtrapolateFactory &GetExtrapolateFactory();
 
class Extrapolate
{
public:
    Extrapolate(const LibUtilities::SessionReaderSharedPtr pSession,
                Array<OneD, MultiRegions::ExpListSharedPtr> pFields,
                MultiRegions::ExpListSharedPtr pPressure,
                const Array<OneD, int> pVel,
                const SolverUtils::AdvectionSharedPtr advObject);
 
    virtual ~Extrapolate();
 
    void GenerateHOPBCMap(
        const LibUtilities::SessionReaderSharedPtr &pSsession);
 
    void UpdateRobinPrimCoeff(void);
 
    inline void SubSteppingTimeIntegration(
        const LibUtilities::TimeIntegrationSchemeSharedPtr &IntegrationScheme);
 
    inline void SubStepSaveFields(const int nstep);
 
    inline void SubStepSetPressureBCs(
        const Array<OneD, const Array<OneD, NekDouble>> &inarray,
        const NekDouble Aii_DT, NekDouble kinvis);
 
    inline void SubStepAdvance(const int nstep, NekDouble time);
 
    inline void MountHOPBCs(int HBCdata, NekDouble kinvis,
                            Array<OneD, NekDouble> &Q,
                            Array<OneD, const NekDouble> &Advection);
 
    inline void EvaluatePressureBCs(
        const Array<OneD, const Array<OneD, NekDouble>> &fields,
        const Array<OneD, const Array<OneD, NekDouble>> &N, NekDouble kinvis);
 
    inline void SetForcing(
        const std::vector<SolverUtils::ForcingSharedPtr> &forcing);
 
    void AddDuDt(void);
 
    void AddVelBC(void);
 
    void ExtrapolatePressureHBCs(void);
    void CopyPressureHBCsToPbndExp(void);
 
    Array<OneD, NekDouble> GetMaxStdVelocity(
        const Array<OneD, Array<OneD, NekDouble>> inarray);
 
    void CorrectPressureBCs(const Array<OneD, NekDouble> &pressure);
 
    void IProductNormVelocityOnHBC(
        const Array<OneD, const Array<OneD, NekDouble>> &Vel,
        Array<OneD, NekDouble> &IprodVn);
 
    void IProductNormVelocityBCOnHBC(Array<OneD, NekDouble> &IprodVn);
 
    std::string GetSubStepName(void);
 
    void ExtrapolateArray(Array<OneD, Array<OneD, NekDouble>> &array);
 
    void EvaluateBDFArray(Array<OneD, Array<OneD, NekDouble>> &array);
 
    void ExtrapolateArray(Array<OneD, Array<OneD, NekDouble>> &oldarrays,
                          Array<OneD, NekDouble> &newarray,
                          Array<OneD, NekDouble> &outarray);
 
    void AddNormVelOnOBC(const int nbcoeffs, const int nreg,
                         Array<OneD, Array<OneD, NekDouble>> &u);
 
    void AddPressureToOutflowBCs(NekDouble kinvis);
 
    void GenerateBndElmtExpansion(void);
 
protected:
    virtual void v_EvaluatePressureBCs(
        const Array<OneD, const Array<OneD, NekDouble>> &inarray,
        const Array<OneD, const Array<OneD, NekDouble>> &N,
        NekDouble kinvis) = 0;
 
    virtual void v_SubSteppingTimeIntegration(
        const LibUtilities::TimeIntegrationSchemeSharedPtr
            &IntegrationScheme) = 0;
 
    virtual void v_SubStepSaveFields(int nstep) = 0;
 
    virtual void v_SubStepSetPressureBCs(
        const Array<OneD, const Array<OneD, NekDouble>> &inarray,
        NekDouble Aii_DT, NekDouble kinvis) = 0;
 
    virtual void v_SubStepAdvance(int nstep, NekDouble time) = 0;
 
    virtual void v_MountHOPBCs(int HBCdata, NekDouble kinvis,
                               Array<OneD, NekDouble> &Q,
                               Array<OneD, const NekDouble> &Advection) = 0;
 
    virtual std::string v_GetSubStepName(void);
 
    virtual void v_AccelerationBDF(Array<OneD, Array<OneD, NekDouble>> &array);
 
    void CalcNeumannPressureBCs(
        const Array<OneD, const Array<OneD, NekDouble>> &fields,
        const Array<OneD, const Array<OneD, NekDouble>> &N, NekDouble kinvis)
    {
        v_CalcNeumannPressureBCs(fields, N, kinvis);
    }
 
    virtual void v_CalcNeumannPressureBCs(
        const Array<OneD, const Array<OneD, NekDouble>> &fields,
        const Array<OneD, const Array<OneD, NekDouble>> &N, NekDouble kinvis);
 
    virtual void v_CorrectPressureBCs(const Array<OneD, NekDouble> &pressure);
 
    virtual void v_AddNormVelOnOBC(const int nbcoeffs, const int nreg,
                                   Array<OneD, Array<OneD, NekDouble>> &u);
 
    void CalcOutflowBCs(const Array<OneD, const Array<OneD, NekDouble>> &fields,
                        NekDouble kinvis);
 
    void RollOver(Array<OneD, Array<OneD, NekDouble>> &input);
 
    LibUtilities::SessionReaderSharedPtr m_session;
 
    LibUtilities::CommSharedPtr m_comm;
 
    /// Array of type of high order BCs for splitting shemes
    Array<OneD, HBCType> m_hbcType;
 
    /// Velocity fields
    Array<OneD, MultiRegions::ExpListSharedPtr> m_fields;
 
    /// Pointer to field holding pressure field
    MultiRegions::ExpListSharedPtr m_pressure;
 
    /// int which identifies which components of m_fields contains the
    /// velocity (u,v,w);
    Array<OneD, int> m_velocity;
 
    SolverUtils::AdvectionSharedPtr m_advObject;
 
    std::vector<SolverUtils::ForcingSharedPtr> m_forcing;
 
    Array<OneD, Array<OneD, NekDouble>> m_previousVelFields;
 
    /// Curl-curl dimensionality
    int m_curl_dim;
 
    /// bounday dimensionality
    int m_bnd_dim;
 
    /// pressure boundary conditions container
    Array<OneD, const SpatialDomains::BoundaryConditionShPtr> m_PBndConds;
 
    /// pressure boundary conditions expansion container
    Array<OneD, MultiRegions::ExpListSharedPtr> m_PBndExp;
 
    /// Boundary expansions on each domain boundary
    Array<OneD, MultiRegions::ExpListSharedPtr> m_bndElmtExps;
 
    /// number of times the high-order pressure BCs have been called
    int m_pressureCalls;
 
    // Number of degrees of freedom (coefficients) for HOPbc
    int m_numHBCDof;
 
    // Number of HOPbcs
    int m_HBCnumber;
 
    /// Maximum points used in pressure BC evaluation
    int m_intSteps;
 
    NekDouble m_timestep;
 
    /// Storage for current and previous levels of high order pressure boundary
    /// conditions.
    Array<OneD, Array<OneD, NekDouble>> m_pressureHBCs;
 
    /// Storage for current and previous levels of the inner product of normal
    /// velocity
    Array<OneD, Array<OneD, NekDouble>> m_iprodnormvel;
 
    Array<OneD, Array<OneD, NekDouble>> m_traceNormals;
 
    // Velocity correction scheme coefficient required for extrapolation.
    static NekDouble StifflyStable_Betaq_Coeffs[3][3];
    static NekDouble StifflyStable_Alpha_Coeffs[3][3];
    static NekDouble StifflyStable_Gamma0_Coeffs[3];
 
    // data related to high order outflow.
    HighOrderOutflowSharedPtr m_houtflow;
 
private:
    static std::string def;
};
 
struct HighOrderOutflow
{
    HighOrderOutflow(const int numHOpts, const int outHBCnumber,
                     const int curldim,
                     const LibUtilities::SessionReaderSharedPtr &pSession)
        : m_numOutHBCPts(numHOpts), m_outHBCnumber(outHBCnumber)
    {
        m_outflowVel =
            Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble>>>>(
                outHBCnumber);
 
        m_outflowVelBnd =
            Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble>>>>(
                outHBCnumber);
 
        m_UBndExp =
            Array<OneD, Array<OneD, MultiRegions::ExpListSharedPtr>>(curldim);
 
        pSession->LoadParameter("OutflowBC_theta", m_obcTheta, 1.0);
        pSession->LoadParameter("OutflowBC_alpha1", m_obcAlpha1, 0.0);
        pSession->LoadParameter("OutflowBC_alpha2", m_obcAlpha2, 0.0);
 
        pSession->LoadParameter("U0_HighOrderBC", m_U0, 1.0);
        pSession->LoadParameter("Delta_HighOrderBC", m_delta, 1 / 20.0);
    }
 
    virtual ~HighOrderOutflow(){};
 
    /// Number of quadrature points for Outflow HOBC
    int m_numOutHBCPts;
 
    /// Number of Outflow HOBCs
    int m_outHBCnumber;
 
    /// Parameters for outflow boundary condition
    NekDouble m_obcTheta;
    NekDouble m_obcAlpha1;
    NekDouble m_obcAlpha2;
    NekDouble m_U0;
    NekDouble m_delta;
    std::string m_defVelPrimCoeff[3];
 
    /// Storage for current and previous velocity fields along the outflow
    Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble>>>> m_outflowVel;
 
    /// Storage for current and previous velocities along the outflow boundary
    Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble>>>>
        m_outflowVelBnd;
 
    /// Velocity boundary condition expansions on high order boundaries.
    Array<OneD, Array<OneD, MultiRegions::ExpListSharedPtr>> m_UBndExp;
 
    /// primitive coefficient for pressure  when using convetive like OBCs
    Array<OneD, NekDouble> m_pressurePrimCoeff;
 
    /// primitive coefficient for velocities when using convetive like OBCs
    Array<OneD, Array<OneD, NekDouble>> m_velocityPrimCoeff;
};
 
/**
 * Evaluate Pressure Boundary Conditions for Standard Extrapolation
 */
inline void Extrapolate::EvaluatePressureBCs(
    const Array<OneD, const Array<OneD, NekDouble>> &inarray,
    const Array<OneD, const Array<OneD, NekDouble>> &N, NekDouble kinvis)
{
    v_EvaluatePressureBCs(inarray, N, kinvis);
}
 
/**
 *
 */
inline void Extrapolate::SubSteppingTimeIntegration(
    const LibUtilities::TimeIntegrationSchemeSharedPtr &IntegrationScheme)
{
    v_SubSteppingTimeIntegration(IntegrationScheme);
}
 
/**
 *
 */
inline void Extrapolate::SubStepSaveFields(int nstep)
{
    v_SubStepSaveFields(nstep);
}
 
/**
 *
 */
inline void Extrapolate::SetForcing(
    const std::vector<SolverUtils::ForcingSharedPtr> &forcing)
{
    m_forcing = forcing;
}
 
/**
 *
 */
inline void Extrapolate::SubStepSetPressureBCs(
    const Array<OneD, const Array<OneD, NekDouble>> &inarray, NekDouble Aii_DT,
    NekDouble kinvis)
{
    v_SubStepSetPressureBCs(inarray, Aii_DT, kinvis);
}
 
/**
 *
 */
inline void Extrapolate::SubStepAdvance(int nstep, NekDouble time)
{
    v_SubStepAdvance(nstep, time);
}
 
/**
 *
 */
inline void Extrapolate::MountHOPBCs(int HBCdata, NekDouble kinvis,
                                     Array<OneD, NekDouble> &Q,
                                     Array<OneD, const NekDouble> &Advection)
{
    v_MountHOPBCs(HBCdata, kinvis, Q, Advection);
}
 
/**
 *
 */
inline std::string Extrapolate::GetSubStepName(void)
{
    return v_GetSubStepName();
}
 
/**
 *
 */
inline void Extrapolate::CorrectPressureBCs(
    const Array<OneD, NekDouble> &pressure)
{
    v_CorrectPressureBCs(pressure);
}
 
/**
 *
 */
inline void Extrapolate::AddNormVelOnOBC(const int nbcoeffs, const int nreg,
                                         Array<OneD, Array<OneD, NekDouble>> &u)
{
    v_AddNormVelOnOBC(nbcoeffs, nreg, u);
}
} // namespace Nektar
 
#endif