advection.h

Go to the documentation of this file.

#ifndef _DG_ADVECTION_H
#define _DG_ADVECTION_H
#include "blas.h"
#include "topology/geometry.h"
#include "enums.h"
#include "topology/evaluation.h"
#include "topology/derivatives.h"
#ifdef MPI_VERSION
#include "topology/mpi_derivatives.h"
#include "topology/mpi_evaluation.h"
#endif
namespace dg
{
 
    //MW this scheme cannot be formulated as a weak form
 
template< class Geometry, class Matrix, class Container >
struct Advection
{
    using geometry_type = Geometry;
    using matrix_type = Matrix;
    using container_type = Container;
    using value_type = get_value_type<Container>;
    Advection() = default;
    Advection( const Geometry& g);
    Advection( const Geometry& g, bc bcx, bc bcy);
 
    template<class ...Params>
    void construct( Params&& ...ps)
    {
        //construct and swap
        *this = Advection( std::forward<Params>( ps)...);
    }
 
    template<class ContainerType0, class ContainerType1, class ContainerType2, class ContainerType3>
    void upwind( value_type alpha, const ContainerType0& vx, const ContainerType1& vy, const ContainerType2& f, value_type beta, ContainerType3& result);
 
  private:
    Container m_temp0, m_temp1;
    Matrix m_dxf, m_dyf, m_dxb, m_dyb;
};
 
template<class Geometry, class Matrix, class Container>
Advection<Geometry, Matrix, Container>::Advection( const Geometry& g ):
    Advection( g, g.bcx(), g.bcy()) { }
 
template<class Geometry, class Matrix, class Container>
Advection<Geometry, Matrix, Container>::Advection( const Geometry& g, bc bcx, bc bcy):
    m_temp0( dg::construct<Container>(dg::evaluate( one, g)) ), m_temp1(m_temp0),
    m_dxf(dg::create::dx( g, bcx, dg::forward)),
    m_dyf(dg::create::dy( g, bcy, dg::forward)),
    m_dxb(dg::create::dx( g, bcx, dg::backward)),
    m_dyb(dg::create::dy( g, bcy, dg::backward))
{
}
 
template< class Geometry, class Matrix, class Container>
template<class ContainerType0, class ContainerType1, class ContainerType2, class ContainerType3>
void Advection<Geometry, Matrix, Container>::upwind( value_type alpha, const ContainerType0& vx, const ContainerType1& vy, const ContainerType2& f, value_type beta, ContainerType3& result)
{
    blas2::symv( m_dxb, f, m_temp0);
    blas2::symv( m_dxf, f, m_temp1);
    blas1::evaluate( result, dg::Axpby<value_type>( alpha, beta), dg::UpwindProduct(), vx, m_temp0, m_temp1);
    blas2::symv( m_dyb, f, m_temp0);
    blas2::symv( m_dyf, f, m_temp1);
    blas1::evaluate( result, dg::Axpby<value_type>( alpha, 1.), dg::UpwindProduct(), vy, m_temp0, m_temp1);
}
 
}//namespace dg
 
 
#endif //_DG_ADVECTION_H