poisson.h

Go to the documentation of this file.

#ifndef _DG_POISSON_CUH
#define _DG_POISSON_CUH
 
#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
{
 
template< class Geometry, class Matrix, class Container >
struct Poisson
{
    using geometry_type = Geometry;
    using matrix_type = Matrix;
    using container_type = Container;
    using value_type = get_value_type<Container>;
    Poisson() = default;
    Poisson( const Geometry& g);
    Poisson( const Geometry& g, bc bcx, bc bcy);
    Poisson( const Geometry& g, bc bcxlhs, bc bcylhs, bc bcxrhs, bc bcyrhs );
    template<class ContainerType0, class ContainerType1, class ContainerType2>
    void operator()( const ContainerType0& lhs, const ContainerType1& rhs, ContainerType2& result);
    template<class ContainerType0>
    void set_chi( const ContainerType0& new_chi) {
        dg::blas1::pointwiseDivide( new_chi, m_perp_vol, m_chi);
    }
 
    const Matrix& dxlhs() const {
        return m_dxlhs;
    }
    const Matrix& dylhs() const {
        return m_dylhs;
    }
    const Matrix& dxrhs() const {
        return m_dxrhs;}
    const Matrix& dyrhs() {
        return m_dyrhs;
    }
 
  private:
    Container m_dxlhslhs, m_dxrhsrhs, m_dylhslhs, m_dyrhsrhs, m_helper;
    Matrix m_dxlhs, m_dylhs, m_dxrhs, m_dyrhs;
    Container m_chi, m_perp_vol;
};
 
//idea: backward transform lhs and rhs and then use bdxf and bdyf , then forward transform
//needs less memory!! and is faster
template< class Geometry, class Matrix, class Container>
Poisson<Geometry, Matrix, Container>::Poisson( const Geometry& g ):
    Poisson( g, g.bcx(), g.bcy(), g.bcx(), g.bcy()){}
 
template< class Geometry, class Matrix, class Container>
Poisson<Geometry, Matrix, Container>::Poisson( const Geometry& g, bc bcx, bc bcy):
    Poisson( g, bcx, bcy, bcx, bcy){}
 
template< class Geometry, class Matrix, class Container>
Poisson<Geometry, Matrix, Container>::Poisson(  const Geometry& g, bc bcxlhs, bc bcylhs, bc bcxrhs, bc bcyrhs):
    m_dxlhslhs( dg::evaluate( one, g) ), m_dxrhsrhs(m_dxlhslhs), m_dylhslhs(m_dxlhslhs), m_dyrhsrhs( m_dxlhslhs), m_helper( m_dxlhslhs),
    m_dxlhs(dg::create::dx( g, bcxlhs, dg::centered)),
    m_dylhs(dg::create::dy( g, bcylhs, dg::centered)),
    m_dxrhs(dg::create::dx( g, bcxrhs, dg::centered)),
    m_dyrhs(dg::create::dy( g, bcyrhs, dg::centered))
{
    m_chi = m_perp_vol = dg::tensor::volume2d(g.metric());
    dg::blas1::pointwiseDivide( 1., m_perp_vol, m_chi);
}
 
template< class Geometry, class Matrix, class Container>
template<class ContainerType0, class ContainerType1, class ContainerType2>
void Poisson< Geometry, Matrix, Container>::operator()( const ContainerType0& lhs, const ContainerType1& rhs, ContainerType2& result)
{
    blas2::symv(  m_dxlhs, lhs,  m_dxlhslhs); //dx_lhs lhs
    blas2::symv(  m_dylhs, lhs,  m_dylhslhs); //dy_lhs lhs
    blas2::symv(  m_dxrhs, rhs,  m_dxrhsrhs); //dx_rhs rhs
    blas2::symv(  m_dyrhs, rhs,  m_dyrhsrhs); //dy_rhs rhs
 
    blas1::pointwiseDot( 1., m_dxlhslhs, m_dyrhsrhs, -1., m_dylhslhs, m_dxrhsrhs, 0., result);
    blas1::pointwiseDot( m_chi, result, result);
}
 
}//namespace dg
 
#endif //_DG_POISSON_CUH