Advection.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: Advection.cpp
//
// 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 advection.
//
///////////////////////////////////////////////////////////////////////////////
 
#include <boost/core/ignore_unused.hpp>
 
#include <SolverUtils/Advection/Advection.h>
 
namespace Nektar
{
namespace SolverUtils
{
 
/**
 * @returns The advection factory.
 */
AdvectionFactory& GetAdvectionFactory()
{
    static AdvectionFactory instance;
    return instance;
}
 
 
/**
 * @param   pSession            Session configuration data.
 * @param   pFields             Array of ExpList objects.
 */
void Advection::InitObject(
    const LibUtilities::SessionReaderSharedPtr        pSession,
    Array<OneD, MultiRegions::ExpListSharedPtr>       pFields)
{
    v_InitObject(pSession, pFields);
}
 
 
/**
 * @param   nConvectiveFields   Number of velocity components.
 * @param   pFields             Expansion lists for scalar fields.
 * @param   pAdvVel             Advection velocity.
 * @param   pInarray            Scalar data to advect.
 * @param   pOutarray           Advected scalar data.
 * @param   pTime               Simulation time.
 */
void Advection::Advect(
    const int                                          nConvectiveFields,
    const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
    const Array<OneD, Array<OneD, NekDouble> >        &pAdvVel,
    const Array<OneD, Array<OneD, NekDouble> >        &pInarray,
    Array<OneD, Array<OneD, NekDouble> >              &pOutarray,
    const NekDouble                                   &pTime,
    const Array<OneD, Array<OneD, NekDouble> >        &pFwd,
    const Array<OneD, Array<OneD, NekDouble> >        &pBwd)
{
    v_Advect(nConvectiveFields, pFields, pAdvVel, pInarray,
            pOutarray, pTime, pFwd, pBwd);
}
 
template<typename DataType, typename TypeNekBlkMatSharedPtr>
void Advection::AddTraceJacToMat(
        const int                                         nConvectiveFields,
        const int                                         nSpaceDim,
        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
        const Array<OneD, TypeNekBlkMatSharedPtr>         &TracePntJacCons,
        Array<OneD, Array<OneD, TypeNekBlkMatSharedPtr> > &gmtxarray,
        const Array<OneD, TypeNekBlkMatSharedPtr>         &TracePntJacGrad,
        const Array<OneD, Array<OneD, DataType> >         &TracePntJacGradSign)
{
    MultiRegions::ExpListSharedPtr tracelist = pFields[0]->GetTrace();
    std::shared_ptr<LocalRegions::ExpansionVector> traceExp 
        = tracelist->GetExp();
    int ntotTrac = (*traceExp).size();
    int nTracePnts = tracelist->GetTotPoints();
    int nTracPnt,nTracCoef;
 
    DNekMatSharedPtr tmp2Add;
    Array<OneD, DataType>    MatData0;
    Array<OneD, NekDouble>    MatData1;
    Array<OneD, NekDouble>    MatData2;
 
    Array<OneD, DNekMatSharedPtr>  TraceJacFwd(ntotTrac);
    Array<OneD, DNekMatSharedPtr>  TraceJacBwd(ntotTrac);
 
    for(int nelmt = 0; nelmt < ntotTrac; nelmt++)
    {
        nTracCoef           = (*traceExp)[nelmt]->GetNcoeffs();
        nTracPnt            = (*traceExp)[nelmt]->GetTotPoints();
        TraceJacFwd[nelmt]    =MemoryManager<DNekMat>
                            ::AllocateSharedPtr(nTracCoef, nTracPnt,0.0);
        TraceJacBwd[nelmt]    =MemoryManager<DNekMat>
                            ::AllocateSharedPtr(nTracCoef, nTracPnt,0.0);
    }
 
    std::shared_ptr<LocalRegions::ExpansionVector> fieldExp 
        = pFields[0]->GetExp();
    int nTotElmt = (*fieldExp).size();
    int nElmtPnt,nElmtCoef;
 
    Array<OneD, DNekMatSharedPtr>  ElmtJacQuad(nTotElmt);
    Array<OneD, DNekMatSharedPtr>  ElmtJacCoef(nTotElmt);
 
    Array<OneD, NekDouble> SymmMatData;
    
    for(int nelmt = 0; nelmt < nTotElmt; nelmt++)
    {
        nElmtCoef          = (*fieldExp)[nelmt]->GetNcoeffs();
        nElmtPnt           = (*fieldExp)[nelmt]->GetTotPoints();
        
        ElmtJacQuad[nelmt]    =MemoryManager<DNekMat>
                            ::AllocateSharedPtr(nElmtCoef, nElmtPnt,0.0);
        ElmtJacCoef[nelmt]    =MemoryManager<DNekMat>
                            ::AllocateSharedPtr(nElmtCoef, nElmtCoef,0.0);
    }
 
    bool TracePntJacGradflag = true;
 
    Array<OneD, Array<OneD, DataType> > TraceJacConsSign(2);
    for(int i = 0; i < 2; i++)
    {
        TraceJacConsSign[i] =   Array<OneD, DataType>(nTracePnts,1.0);
    }
 
    if(0==TracePntJacGrad.size())
    {
        TracePntJacGradflag = false;
    }
 
    for(int m = 0; m < nConvectiveFields; m++)
    {
        for(int n = 0; n < nConvectiveFields; n++)
        {
            // ElmtJacCons to set 0
            for(int nelmt = 0; nelmt < nTotElmt; nelmt++)
            {
                (*ElmtJacCoef[nelmt]) =  0.0;
                (*ElmtJacQuad[nelmt]) =  0.0;
            }
 
            if(TracePntJacGradflag)
            {
                //only BJac is stored here because BJac = - FJac
                for(int ndir = 0; ndir < nSpaceDim; ndir++)
                {
                    // ElmtJacGrad to set 0
                    for(int nelmt = 0; nelmt < nTotElmt; nelmt++)
                    {
                        (*ElmtJacQuad[nelmt]) =  0.0;
                    }
                    int ngrad = n*nSpaceDim+ndir;
                    
                    CalcJacobTraceInteg(pFields, m,ngrad, TracePntJacGrad, 
                        TracePntJacGradSign, TraceJacFwd, TraceJacBwd);
                    pFields[0]->AddTraceJacToElmtJac(TraceJacFwd, TraceJacBwd, 
                        ElmtJacQuad);
                    // cost can be lower down by doing 3 directions together
                    pFields[0]->AddRightIPTPhysDerivBase(ndir, ElmtJacQuad, 
                        ElmtJacCoef);
                }
                for(int nelmt = 0; nelmt < nTotElmt; nelmt++)
                {
                    nElmtCoef          = (*fieldExp)[nelmt]->GetNcoeffs();
                    if(SymmMatData.size()<nElmtCoef)
                    {
                        SymmMatData = 
                            Array<OneD, NekDouble> (nElmtCoef*nElmtCoef);
                    }
                    MatData1 =  ElmtJacCoef[nelmt]->GetPtr();
                    for(int i=0;i<nElmtCoef;i++)
                    {
                        Vmath::Vcopy(nElmtCoef, &MatData1[i*nElmtCoef], 1, 
                            &SymmMatData[i], nElmtCoef);
                    }
                    Vmath::Vadd(nElmtCoef*nElmtCoef, SymmMatData, 1, 
                        MatData1, 1, MatData1, 1);
                }
            }
 
            CalcJacobTraceInteg(pFields, m, n, TracePntJacCons, 
                TraceJacConsSign, TraceJacFwd, TraceJacBwd);
                
            pFields[0]->AddTraceJacToElmtJac(TraceJacFwd, TraceJacBwd,
                ElmtJacQuad);
            //Memory can be save to explore reusing ElmtJacQuad as output
            pFields[0]->AddRightIPTBaseMatrix(ElmtJacQuad, ElmtJacCoef);
 
            // add ElmtJacCons to gmtxarray[m][n]
            for(int nelmt = 0; nelmt < nTotElmt; nelmt++)
            {
                tmp2Add = ElmtJacCoef[nelmt];
                MatData0 = gmtxarray[m][n]->GetBlock(nelmt,nelmt)->GetPtr();
                MatData1 = tmp2Add->GetPtr();
                for(int i=0;i<MatData0.size();i++)
                {
                    MatData0[i] += DataType(MatData1[i]);
                }
            }
        }
    }
}
 
template SOLVER_UTILS_EXPORT void Advection::AddTraceJacToMat( 
        const int                                         nConvectiveFields,
        const int                                         nSpaceDim,
        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
        const Array<OneD, DNekBlkMatSharedPtr>            &TracePntJacCons,
        Array<OneD, Array<OneD, DNekBlkMatSharedPtr> >    &gmtxarray,
        const Array<OneD, DNekBlkMatSharedPtr>            &TracePntJacGrad,
        const Array<OneD, Array<OneD, NekDouble> >        &TracePntJacGradSign);
template SOLVER_UTILS_EXPORT void Advection::AddTraceJacToMat( 
        const int                                         nConvectiveFields,
        const int                                         nSpaceDim,
        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
        const Array<OneD, SNekBlkMatSharedPtr>            &TracePntJacCons,
        Array<OneD, Array<OneD, SNekBlkMatSharedPtr> >    &gmtxarray,
        const Array<OneD, SNekBlkMatSharedPtr>            &TracePntJacGrad,
        const Array<OneD, Array<OneD, NekSingle> >        &TracePntJacGradSign);
 
 
// TODO: To change it into one by one
template<typename DataType, typename TypeNekBlkMatSharedPtr>
void Advection::CalcJacobTraceInteg(
    const Array<OneD, MultiRegions::ExpListSharedPtr>   &pFields,
    const int                                             m,
    const int                                             n,
    const Array<OneD, const TypeNekBlkMatSharedPtr>       & PntJac,
    const Array<OneD, const Array<OneD, DataType> >       & PntJacSign,
    Array<OneD, DNekMatSharedPtr>                         & TraceJacFwd,
    Array<OneD, DNekMatSharedPtr>                         & TraceJacBwd)
{
    MultiRegions::ExpListSharedPtr tracelist = pFields[0]->GetTrace();
    std::shared_ptr<LocalRegions::ExpansionVector> traceExp = 
        tracelist->GetExp();
    int ntotTrac            = (*traceExp).size();
    int nTracPnt,noffset,pntoffset;
 
    Array<OneD, int > tracepnts(ntotTrac);
    Array<OneD, Array<OneD, NekDouble> > JacFwd(ntotTrac);
    Array<OneD, Array<OneD, NekDouble> > JacBwd(ntotTrac);
 
    for(int  nelmt = 0; nelmt < ntotTrac; nelmt++)
    {
        nTracPnt            = (*traceExp)[nelmt]->GetTotPoints();
        tracepnts[nelmt]     =   nTracPnt;
 
        JacFwd[nelmt]     =Array<OneD, NekDouble>(nTracPnt,0.0);
        JacBwd[nelmt]     =Array<OneD, NekDouble>(nTracPnt,0.0);
    }
 
    // DNekMatSharedPtr FtmpMat,BtmpMat;
    DataType ftmp,btmp;
    for(int  nelmt = 0; nelmt < ntotTrac; nelmt++)
    {
        nTracPnt            =   tracepnts[nelmt];
        noffset             =   tracelist->GetPhys_Offset(nelmt);
        for(int npnt = 0; npnt < nTracPnt; npnt++)
        {
            pntoffset = noffset+npnt;
            ftmp = (*PntJac[0]->GetBlock(pntoffset,pntoffset))(m,n);
            JacFwd[nelmt][npnt] = NekDouble(PntJacSign[0][pntoffset]*ftmp);
 
            btmp = (*PntJac[1]->GetBlock(pntoffset,pntoffset))(m,n);
            JacBwd[nelmt][npnt] = NekDouble(PntJacSign[1][pntoffset]*btmp);
        }
    }
    tracelist->GetDiagMatIpwrtBase(JacFwd,TraceJacFwd);
    tracelist->GetDiagMatIpwrtBase(JacBwd,TraceJacBwd);
}
 
 
void Advection::v_AdvectVolumeFlux(
    const int nConvectiveFields,
    const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
    const Array<OneD, Array<OneD, NekDouble>>         &pAdvVel,
    const Array<OneD, Array<OneD, NekDouble>>         &pInarray,
    TensorOfArray3D<NekDouble>                        &pVolumeFlux,
    const NekDouble                                   &pTime)
{
    boost::ignore_unused(nConvectiveFields, pFields, pAdvVel, pInarray,
                        pVolumeFlux, pTime);
    ASSERTL0(false, "Not defined for AdvectVolumeFlux.");
}
 
/**
 * @brief calculate the advection flux in the cell the trace  integration
 */
void Advection::v_AdvectTraceFlux(
    const int nConvectiveFields,
    const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
    const Array<OneD, Array<OneD, NekDouble>>         &pAdvVel,
    const Array<OneD, Array<OneD, NekDouble>>         &pInarray,
    Array<OneD, Array<OneD, NekDouble>>               &pTraceFlux,
    const NekDouble                                   &pTime,
    const Array<OneD, Array<OneD, NekDouble>>         &pFwd,
    const Array<OneD, Array<OneD, NekDouble>>         &pBwd)
{
    boost::ignore_unused(nConvectiveFields, pFields, pAdvVel, pInarray,
                        pTraceFlux, pTime, pFwd, pBwd);
    ASSERTL0(false, "Not defined for AdvectTraceFlux.");
}
 
/**
 * @brief Similar with Advection::Advect(): calculate the advection flux
 * The difference is in the outarray:
 *  it is the coefficients of basis for AdvectCoeffs()
 *  it is the physics on quadrature points for Advect()
 *
 * @param   nConvectiveFields   Number of velocity components.
 * @param   pFields             Expansion lists for scalar fields.
 * @param   pAdvVel             Advection velocity.
 * @param   pInarray            Scalar data to advect.
 * @param   pOutarray           Advected scalar data.
 * @param   pTime               Simulation time.
 */
void Advection::AdvectCoeffs(
    const int                                          nConvectiveFields,
    const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
    const Array<OneD, Array<OneD, NekDouble> >        &pAdvVel,
    const Array<OneD, Array<OneD, NekDouble> >        &pInarray,
    Array<OneD, Array<OneD, NekDouble> >              &pOutarray,
    const NekDouble                                   &pTime,
    const Array<OneD, Array<OneD, NekDouble> >        &pFwd,
    const Array<OneD, Array<OneD, NekDouble> >        &pBwd)
{
    v_AdvectCoeffs(nConvectiveFields, pFields, pAdvVel, pInarray,
                   pOutarray, pTime, pFwd, pBwd);
}
 
/**
 * 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.
 *
 * @param   pSession            Session information.
 * @param   pFields             Expansion lists for scalar fields.
 */
void Advection::v_InitObject(
    const LibUtilities::SessionReaderSharedPtr        pSession,
    Array<OneD, MultiRegions::ExpListSharedPtr>       pFields)
{
    m_spaceDim = pFields[0]->GetCoordim(0);
 
    if (pSession->DefinesSolverInfo("HOMOGENEOUS"))
    {
        std::string HomoStr = pSession->GetSolverInfo("HOMOGENEOUS");
        if (HomoStr == "HOMOGENEOUS1D" || HomoStr == "Homogeneous1D" ||
            HomoStr == "1D"            || HomoStr == "Homo1D" ||
            HomoStr == "HOMOGENEOUS2D" || HomoStr == "Homogeneous2D" ||
            HomoStr == "2D"            || HomoStr == "Homo2D")
        {
            m_spaceDim++;
        }
        else
        {
            NEKERROR(ErrorUtil::efatal,
                     "Only 1D homogeneous dimension supported.");
        }
    }
}
 
 
/**
 *
 */
void Advection::v_SetBaseFlow(
        const Array<OneD, Array<OneD, NekDouble> >        &inarray,
        const Array<OneD, MultiRegions::ExpListSharedPtr> &fields)
{
    boost::ignore_unused(inarray, fields);
    NEKERROR(ErrorUtil::efatal,
            "A baseflow is not appropriate for this advection type.");
}
 
void Advection::v_AdvectCoeffs(
    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,
    const Array<OneD, Array<OneD, NekDouble> >        &pFwd,
    const Array<OneD, Array<OneD, NekDouble> >        &pBwd)
{
    boost::ignore_unused(nConvectiveFields, fields, advVel, inarray, outarray,
                        time, pFwd, pBwd);
    ASSERTL0(false, "v_AdvectCoeffs not defined");
}
 
void Advection::v_AddVolumJacToMat( 
        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
        const int &nConvectiveFields,
        const TensorOfArray5D<NekDouble> &ElmtJacArray,
        Array<OneD, Array<OneD, SNekBlkMatSharedPtr>> &gmtxarray)
    {
        boost::ignore_unused(pFields,nConvectiveFields,ElmtJacArray,gmtxarray);
        ASSERTL0(false,"v_AddVolumJacToMat NOT SPECIFIED");
        return;
    }
 
}
}