app-fvdgfem/advection.hh

Go to the documentation of this file.

#ifndef advection_hh
#define advection_hh

#include "basics/typedefs.hh"
#include "toolbox/inputOutput.hh"
#include "formula/boundary.hh"
#include "geometry/boundaryConditions.hh"
#include "operator/compositions.hh"
#include "operator/diagonal.hh"
#include "timestepping/vectors.hh"
#include "linDG3D/element.hh"
#include "linDG3D/space.hh"
#include "linDG3D/elementPair.hh"

#include <map>
#include <vector>
#include <utility>
#include <numeric>

#define DEBUG_ADVECTION 0
#define DEBUG_FLUXES 0
#define DEBUG_ADAPT_LIM_ELM 0
#define DEBUG_ADAPT_LIM_INDICATOR 0
#define DEBUG_REF_ASSIGNMENT 0

using concepts::Real;


// ****************************************************************** Fluxes **
class BurgerFlux {
public:
  BurgerFlux() :
    normal_(concepts::Real3d(1., 1., 1.)) {}
  concepts::Real3d flux(const Real u) const {
    return normal_*(0.5*u*u); 
  }
  Real flux(const Real u, const concepts::Real3d& direction) const {
    return normal_*direction*(0.5*u*u);
  }
  concepts::Real3d deriv(const Real u) const {
    return normal_*u;
  }
  Real deriv(const Real u, const concepts::Real3d& direction) const {
    return normal_*direction*u;
  }
  bool isConvex(const concepts::Real3d& normal) const {
    return normal*normal_ > 0.;
  }
  Real uMin(const concepts::Real3d& normal) const { return 0; }

private:
  const concepts::Real3d normal_;
};

class LinearAdvectionFlux {
 public:
  LinearAdvectionFlux() :
    normal_(concepts::Real3d(1., 1., 1.)) {}
  concepts::Real3d flux(const Real u) const { return normal_*u;}
  Real flux(const Real u, const concepts::Real3d& direction) const {
    return normal_*direction*u;
  }
  concepts::Real3d deriv(const Real u) const { return normal_; }
  Real deriv(const Real u, const concepts::Real3d& direction) const {
    return normal_*direction;
  }
private:
  const concepts::Real3d normal_;
};


// *************************************************************** Advection **

template<typename NumFlux>
class Advection : public concepts::Operator<Real> {
public:
  Advection(const NumFlux& fvFlux,
            const linDG3D::FvdgSpace& spc,
            const bool pcb,
            concepts::BoundaryConditions* bc = 0)
    : concepts::Operator<Real>(spc.dim(), spc.dim())
    , numFlux_(fvFlux),
      spc_(spc),
      innerList_(spc.innerElmPairList()),
      boundaryList_(spc.boundaryElmPairList()),
      time_(0),
      pcb_(pcb),
      bc_(bc) {};
  void operator()
    (const concepts::Function<Real>& fncY, concepts::Function<Real>& fncX);
  const concepts::Space<Real>& spaceX() const { return spc_; }
  const concepts::Space<Real>& spaceY() const { return spc_; }
protected:
  virtual std::ostream& info (std::ostream &os) const;
private:
  const NumFlux& numFlux_;
  const linDG3D::FvdgSpace& spc_;
  const concepts::ElementPairList<Real>& innerList_;
  const concepts::ElementPairList<Real>& boundaryList_;
  Real time_;
  const bool pcb_;
  concepts::BoundaryConditions* bc_;
  Real uBoundary_(const linDG3D::FvdgElement& elm,
      const uint faceIdx,
      const concepts::Real3d &x) const;
};

// Numerical Fluxes

//************************************************************ LaxFriedrichs **
template <typename T>
class LaxFriedrichs {
public:
  friend class Advection<LaxFriedrichs>;
  LaxFriedrichs() : fluxFct_() {}
  Real operator()
    (Real u1, Real u2, const concepts::Real3d& normal) const;
private:
  const T fluxFct_;
};

//******************************************************************* Upwind **
template <typename T>
class Upwind {
  friend class Advection<Upwind>;
public:
  Upwind() : fluxFct_() {}
  Real operator()
    (Real u1, Real u2, const concepts::Real3d& normal) const;
private:
  const T fluxFct_;
};

//****************************************************************** Godunov **
template <typename T>
class Godunov {
public:
  friend class Advection<Godunov>;
  Godunov() : fluxFct_() {}
  Real operator() (Real u1, Real u2,
       const concepts::Real3d& normal) const;
private:
  const T fluxFct_;
};

//******************************************************************** Osher **
template <typename T>
class Osher {
public:
  friend class Advection<Osher>;
  Osher() : fluxFct_() {}
  Real operator()
    (Real u1, Real u2, const concepts::Real3d& normal) const;
private:
  const T fluxFct_;
};

// Limiter

//*********************************************************** FvdgPi0Limiter **
class FvdgPi0Limiter : public concepts::Operator<Real> {
public:
  FvdgPi0Limiter(const linDG3D::FvdgSpace& spc);
  virtual void operator() (const concepts::Function<Real> &fncY,
         concepts::Function<Real> &fncX) {
    diag_(fncY, fncX);
  }
protected:
  concepts::DiagonalMatrix<Real> diag_;
};

//*********************************************************** FvdgPi0Limiter **
class FvdgAdaptiveLimiter : public concepts::Operator<Real> {
public:
  FvdgAdaptiveLimiter(const linDG3D::FvdgSpace& spc,
          const Real alpha,
          const Real lambda,
          concepts::BoundaryConditions* bc);
  virtual const concepts::Space<Real>& spaceX() const {
    return spc_;
  }
  virtual const concepts::Space<Real>& spaceY() const {
    return spc_; 
  }
  virtual void operator() (const concepts::Function<Real> &fncY,
         concepts::Function<Real> &fncX);
  void getIndicator(concepts::Vector<Real> &indicator);
protected:
  const linDG3D::FvdgSpace& spc_;
  concepts::BoundaryConditions* bc_;
  const Real alpha_;
  Real lambda_;
  std::map<const linDG3D::FvdgElement*, Real> elmIndicator_;
  // Private methods
  void limit_(concepts::Vector<Real>& vecX,
        const linDG3D::FvdgElement& elm) const;
  void reset_();
  void innerContrib_(const concepts::Vector<Real>& vecY, 
         concepts::Vector<Real>& vecX);
  void bndryContrib_(const concepts::Vector<Real>& vecY, 
         concepts::Vector<Real>& vecX);
  Real elmVolume_(const linDG3D::FvdgElement& elm) const;
};

#endif // advection_hh