PulseWaveSystem.h

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File PulseWaveSystem.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).
//
// 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: Generic timestepping for PulseWaveSolver
//
///////////////////////////////////////////////////////////////////////////////
 
#ifndef NEKTAR_SOLVERS_PULSEWAVESOLVER_EQUATIONSYSTEMS_PULSEWAVESYSTEM_H
#define NEKTAR_SOLVERS_PULSEWAVESOLVER_EQUATIONSYSTEMS_PULSEWAVESYSTEM_H
 
#include <SolverUtils/UnsteadySystem.h>
#include <PulseWaveSolver/EquationSystems/PulseWavePressureArea.h>
 
using namespace Nektar::SolverUtils;
 
namespace Nektar
{
 
    enum UpwindTypePulse
    {
        eNotSetPulse,            ///< flux not defined
        eUpwindPulse,            ///< simple upwinding scheme
        SIZE_UpwindTypePulse     ///< Length of enum list
    };
 
    const char* const UpwindTypeMapPulse[] =
        {
            "NoSetPulse",
            "UpwindPulse",
        };
 
    struct InterfacePoint
    {
        InterfacePoint(const int vid,
                       const int domain,
                       const int elmt,
                       const int elmtVert,
                       const int traceId,
                       const int bcpos):
        m_vid(vid),
            m_domain(domain),
            m_elmt(elmt),
            m_elmtVert(elmtVert),
            m_traceId(traceId),
            m_bcPosition(bcpos)
        {
        };
        int m_vid;        // Global Vid of interface point
        int m_domain;     // domain interface point belongs to
        int m_elmt;       // element id of vertex
        int m_elmtVert;   // vertex id within local element
        int m_traceId;    // Element id  within the trace
        int m_bcPosition; // Position of boundary condition in region
    };
 
    typedef std::shared_ptr<InterfacePoint> InterfacePointShPtr;
 
    /// Base class for unsteady solvers.
    class PulseWaveSystem : public UnsteadySystem
    {
    public:
        /// Destructor
        virtual ~PulseWaveSystem();
 
        int  GetNdomains()
        {
            return m_nDomains;
        }
 
        Array<OneD, MultiRegions::ExpListSharedPtr> UpdateVessels(void)
        {
            return m_vessels;
        }
 
    protected:
        Array<OneD, MultiRegions::ExpListSharedPtr>     m_vessels;
        int                                             m_nDomains;
        int                                             m_currentDomain;
        int                                             m_nVariables;
        UpwindTypePulse                                 m_upwindTypePulse;
 
        Array<OneD, int>                                m_fieldPhysOffset;
        NekDouble                                       m_rho;
        NekDouble                                       m_pext;
 
        NekDouble m_C;
        NekDouble m_RT;
        NekDouble m_pout;
 
        Array<OneD, Array<OneD, NekDouble> >        m_A_0;
        Array<OneD, Array<OneD, NekDouble> >        m_A_0_trace;
        Array<OneD, Array<OneD, NekDouble> >        m_beta;
        Array<OneD, Array<OneD, NekDouble> >        m_beta_trace;
        Array<OneD, Array<OneD, NekDouble> >        m_gamma;
        Array<OneD, Array<OneD, NekDouble> >        m_alpha;
        Array<OneD, Array<OneD, NekDouble> >        m_alpha_trace;
        Array<OneD, Array<OneD, NekDouble> >        m_trace_fwd_normal;
 
        Array<OneD, Array<OneD, NekDouble> >        m_pressure;
        PulseWavePressureAreaSharedPtr              m_pressureArea;
 
        // Extra Fields
        bool extraFields = false;
        Array<OneD, Array<OneD, NekDouble> >        m_PWV;
        Array<OneD, Array<OneD, NekDouble> >        m_W1;
        Array<OneD, Array<OneD, NekDouble> >        m_W2;
 
        std::vector<std::vector<InterfacePointShPtr> >  m_vesselJcts;
        std::vector<std::vector<InterfacePointShPtr> >  m_bifurcations;
        std::vector<std::vector<InterfacePointShPtr> >  m_mergingJcts;
 
        /// Initialises PulseWaveSystem class members.
        PulseWaveSystem(const LibUtilities::SessionReaderSharedPtr& pSession,
                        const SpatialDomains::MeshGraphSharedPtr& pGraph);
 
        virtual void v_InitObject();
 
        /// Sets up initial conditions.
        virtual void v_DoInitialise();
 
        /// Solves an unsteady problem.
        virtual void v_DoSolve();
 
        /// Links the subdomains
        void LinkSubdomains(Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
                                                                       &fields);
 
        /// Riemann Problem for Bifurcation
        void BifurcationRiemann(Array<OneD, NekDouble> &Au, Array<OneD, NekDouble> &uu,
                      Array<OneD, NekDouble> &beta, Array<OneD, NekDouble> &A_0,
                                                 Array<OneD, NekDouble> &alpha);
 
        /// Riemann Problem for Merging Flow
        void MergingRiemann(Array<OneD, NekDouble> &Au, Array<OneD, NekDouble> &uu,
                      Array<OneD, NekDouble> &beta, Array<OneD, NekDouble> &A_0,
                                                 Array<OneD, NekDouble> &alpha);
 
        /// Riemann Problem for Junction
        void JunctionRiemann(Array<OneD, NekDouble> &Au, Array<OneD, NekDouble> &uu,
                      Array<OneD, NekDouble> &beta, Array<OneD, NekDouble> &A_0,
                                                 Array<OneD, NekDouble> &alpha);
 
        // Ouptut field information
        virtual void v_Output(void);
 
        // Checkpoint field output
        void CheckPoint_Output(const int n);
 
        /// Compute the L2 error between fields and a given exact solution.
        NekDouble v_L2Error(unsigned int field,
                  const Array<OneD,NekDouble> &exactsoln = NullNekDouble1DArray,
                                                       bool Normalised = false);
 
        /// Compute the L_inf error between fields and a given exact solution.
        NekDouble v_LinfError(unsigned int field,
                 const Array<OneD,NekDouble> &exactsoln = NullNekDouble1DArray);
 
        /// Write input fields to the given filename.
        void WriteVessels(const std::string &outname);
 
        void EnforceInterfaceConditions(const Array<OneD, const Array<OneD, NekDouble> > &fields);
 
    private:
        void SetUpDomainInterfaces(void);
        void FillDataFromInterfacePoint(InterfacePointShPtr &I,
           const Array<OneD, const Array<OneD, NekDouble> >&field, NekDouble &A,
               NekDouble &u, NekDouble &beta, NekDouble &A_0, NekDouble &alpha);
    };
 
    typedef std::shared_ptr<PulseWaveSystem> PulseWaveSystemSharedPtr;
}
 
#endif