EulerADCFE.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File EulerArtificialDiffusionCFE.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).
//
// License for the specific language governing rights and limitations under
// 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: Euler equations in conservative variables with artificial
// diffusion
//
///////////////////////////////////////////////////////////////////////////////

#include <CompressibleFlowSolver/EquationSystems/EulerADCFE.h>
#include <boost/algorithm/string.hpp>

using namespace std;

namespace Nektar
{
    string EulerADCFE::className =
    SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
        "EulerADCFE", EulerADCFE::create,
        "Euler equations in conservative variables with "
        "artificial diffusion.");

    EulerADCFE::EulerADCFE(
        const LibUtilities::SessionReaderSharedPtr& pSession)
    : CompressibleFlowSystem(pSession)
    {
    }

    void EulerADCFE::v_InitObject()
    {
        CompressibleFlowSystem::v_InitObject();

        if (m_shockCaptureType == "Smooth")
        {
            ASSERTL0(m_fields.num_elements() == m_spacedim + 3,
                     "Not enough variables for smooth shock capturing; "
                     "make sure you have added eps to variable list.");
            m_smoothDiffusion = true;
        }

        m_diffusion->SetArtificialDiffusionVector(
            &EulerADCFE::GetArtificialDynamicViscosity, this);

        if(m_session->DefinesSolverInfo("PROBLEMTYPE"))
        {
            std::string ProblemTypeStr = m_session->GetSolverInfo("PROBLEMTYPE");
            int i;
            for(i = 0; i < (int) SIZE_ProblemType; ++i)
            {
                if(boost::iequals(ProblemTypeMap[i], ProblemTypeStr))
                {
                    m_problemType = (ProblemType)i;
                    break;
                }
            }
        }
        else
        {
            m_problemType = (ProblemType)0;
        }

        if (m_explicitAdvection)
        {
            m_ode.DefineOdeRhs    (&EulerADCFE::
            DoOdeRhs, this);
            m_ode.DefineProjection(&EulerADCFE::
            DoOdeProjection, this);
        }
        else
        {
            ASSERTL0(false, "Implicit CFE not set up.");
        }
    }

    EulerADCFE::~EulerADCFE()
    {

    }

    void EulerADCFE::v_GenerateSummary(SolverUtils::SummaryList& s)
    {
        CompressibleFlowSystem::v_GenerateSummary(s);
        SolverUtils::AddSummaryItem(
            s, "Problem Type", ProblemTypeMap[m_problemType]);
    }

    void EulerADCFE::v_SetInitialConditions(
        NekDouble initialtime,
        bool      dumpInitialConditions,
        const int domain)
    {
        EquationSystem::v_SetInitialConditions(initialtime, false);
        CompressibleFlowSystem::v_SetInitialConditions();

        if(dumpInitialConditions)
        {
            // Dump initial conditions to file
            Checkpoint_Output(0);
        }
    }

    void EulerADCFE::DoOdeRhs(
        const Array<OneD, const Array<OneD, NekDouble> > &inarray,
        Array<OneD,       Array<OneD, NekDouble> > &outarray,
        const NekDouble                                   time)
    {
        int i;
        int nvariables = inarray.num_elements();
        int npoints    = GetNpoints();

        Array<OneD, Array<OneD, NekDouble> > advVel;
        Array<OneD, Array<OneD, NekDouble> > outarrayAdv(nvariables);
        Array<OneD, Array<OneD, NekDouble> > outarrayDiff(nvariables);

        for (i = 0; i < nvariables; ++i)
        {
            outarrayAdv[i] = Array<OneD, NekDouble>(npoints, 0.0);
            outarrayDiff[i] = Array<OneD, NekDouble>(npoints, 0.0);
        }

        m_advection->Advect(nvariables, m_fields, advVel, inarray,
                                        outarrayAdv, m_time);

        for (i = 0; i < nvariables; ++i)
        {
            Vmath::Neg(npoints, outarrayAdv[i], 1);
        }

        m_diffusion->Diffuse(nvariables, m_fields, inarray, outarrayDiff);

        if (m_shockCaptureType == "NonSmooth")
        {
            for (i = 0; i < nvariables; ++i)
            {
                Vmath::Vadd(npoints,
                            outarrayAdv[i], 1,
                            outarrayDiff[i], 1,
                            outarray[i], 1);
            }
        }
        if(m_shockCaptureType == "Smooth")
        {
            const Array<OneD, int> ExpOrder = GetNumExpModesPerExp();

            NekDouble pOrder = Vmath::Vmax(ExpOrder.num_elements(), ExpOrder, 1);

            Array <OneD, NekDouble > a_vel  (npoints, 0.0);
            Array <OneD, NekDouble > u_abs  (npoints, 0.0);
            Array <OneD, NekDouble > pres   (npoints, 0.0);
            Array <OneD, NekDouble > wave_sp(npoints, 0.0);

            GetPressure(inarray, pres);
            GetSoundSpeed(inarray, pres, a_vel);
            GetAbsoluteVelocity(inarray, u_abs);

            Vmath::Vadd(npoints, a_vel, 1, u_abs, 1, wave_sp, 1);

            NekDouble max_wave_sp = Vmath::Vmax(npoints, wave_sp, 1);

            Vmath::Smul(npoints,
                        m_C2,
                        outarrayDiff[nvariables-1], 1,
                        outarrayDiff[nvariables-1], 1);

            Vmath::Smul(npoints,
                        max_wave_sp,
                        outarrayDiff[nvariables-1], 1,
                        outarrayDiff[nvariables-1], 1);

            Vmath::Smul(npoints,
                        pOrder,
                        outarrayDiff[nvariables-1], 1,
                        outarrayDiff[nvariables-1], 1);

            for (i = 0; i < nvariables; ++i)
            {
                Vmath::Vadd(npoints,
                            outarrayAdv[i], 1,
                            outarrayDiff[i], 1,
                            outarray[i], 1);
            }

            Array<OneD, Array<OneD, NekDouble> > outarrayForcing(nvariables);

            for (i = 0; i < nvariables; ++i)
            {
                outarrayForcing[i] = Array<OneD, NekDouble>(npoints, 0.0);
            }

            GetForcingTerm(inarray, outarrayForcing);

            for (i = 0; i < nvariables; ++i)
            {
                // Add Forcing Term
                Vmath::Vadd(npoints,
                            outarray[i], 1,
                            outarrayForcing[i], 1,
                            outarray[i], 1);
            }
        }

        // Add sponge layer if defined in the session file
        std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
        for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
        {
            (*x)->Apply(m_fields, inarray, outarray, time);
        }
    }

    void EulerADCFE::DoOdeProjection(
        const Array<OneD, const Array<OneD, NekDouble> > &inarray,
        Array<OneD,       Array<OneD, NekDouble> > &outarray,
        const NekDouble                                   time)
    {
        int i;
        int nvariables = inarray.num_elements();

        switch(m_projectionType)
        {
            case MultiRegions::eDiscontinuous:
            {
                // Just copy over array
                int npoints = GetNpoints();

                for(i = 0; i < nvariables; ++i)
                {
                    Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
                }
                SetBoundaryConditions(outarray, time);
                break;
            }
            case MultiRegions::eGalerkin:
            case MultiRegions::eMixed_CG_Discontinuous:
            {
                ASSERTL0(false, "No Continuous Galerkin for Euler equations");
                break;
            }
            default:
                ASSERTL0(false, "Unknown projection scheme");
                break;
        }
    }

    void EulerADCFE::SetBoundaryConditions(
        Array<OneD, Array<OneD, NekDouble> > &inarray,
        NekDouble                             time)
    {
        std::string varName;
        int cnt        = 0;
        int nTracePts  = GetTraceTotPoints();
        int nvariables = inarray.num_elements();

        Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
        for (int i = 0; i < nvariables; ++i)
        {
            Fwd[i] = Array<OneD, NekDouble>(nTracePts);
            m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
        }

        // loop over Boundary Regions
        for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
        {
            std::string type = m_fields[0]->GetBndConditions()[n]->GetUserDefined();

            // Wall Boundary Condition
            if (boost::iequals(type,"WallViscous"))
            {
                // Wall Boundary Condition
                ASSERTL0(false, "WallViscous is a wrong bc for the "
                         "Euler equations");
            }
            else
            {
                SetCommonBC(type, n, time, cnt, Fwd, inarray);
            }

            // no User Defined conditions provided so skip cnt
            // this line is left in case solver specific condition is added.
            cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
        }
    }
}