stabilization.h

Go to the documentation of this file.

/*
 * Copyright (c) 2010-2015:  G-CSC, Goethe University Frankfurt
 * Author: Andreas Vogel
 * 
 * This file is part of UG4.
 * 
 * UG4 is free software: you can redistribute it and/or modify it under the
 * terms of the GNU Lesser General Public License version 3 (as published by the
 * Free Software Foundation) with the following additional attribution
 * requirements (according to LGPL/GPL v3 §7):
 * 
 * (1) The following notice must be displayed in the Appropriate Legal Notices
 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (2) The following notice must be displayed at a prominent place in the
 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (3) The following bibliography is recommended for citation and must be
 * preserved in all covered files:
 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
 *   parallel geometric multigrid solver on hierarchically distributed grids.
 *   Computing and visualization in science 16, 4 (2013), 151-164"
 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
 *   flexible software system for simulating pde based models on high performance
 *   computers. Computing and visualization in science 16, 4 (2013), 165-179"
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU Lesser General Public License for more details.
 */
 
#ifndef __H__UG__PLUGINS__NAVIER_STOKES__INCOMPRESSIBLE__FV1__STABILIZATION__
#define __H__UG__PLUGINS__NAVIER_STOKES__INCOMPRESSIBLE__FV1__STABILIZATION__
 
//#define UG_NSSTAB_ASSERT(cond, exp)
// include define below to assert arrays used in stabilization
#define UG_NSSTAB_ASSERT(cond, exp) UG_ASSERT((cond), (exp))
 
#include "../../upwind_interface.h"
#include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
#include "lib_disc/spatial_disc/user_data/data_import.h"
 
 
namespace ug{
namespace NavierStokes{
 
template <int dim>
class INavierStokesFV1Stabilization
{
  protected:
    typedef fv1_dim_traits<dim, dim> traits;
 
  public:
    static const size_t maxNumSCV = traits::maxNumSCV;
 
    static const size_t maxNumSCVF = traits::maxNumSCVF;
 
    static const size_t maxNumSH = traits::maxNSH;
 
  private:
    typedef INavierStokesFV1Stabilization<dim> this_type;
 
  public:
    INavierStokesFV1Stabilization()
    : m_numScvf(0), m_numSh(0),
      m_spUpwind(NULL)
    {
      m_vUpdateFunc.clear();
    }
 
    void set_upwind(SmartPtr<INavierStokesUpwind<dim> > spUpwind)
    {
      m_spUpwind = spUpwind;
      m_spConstUpwind = spUpwind;
    }
 
    const ConstSmartPtr<INavierStokesUpwind<dim> >& upwind() const {return m_spConstUpwind;}
    
    bool upwind_valid() const {return m_spConstUpwind.valid();}
 
    template <typename TFVGeom>
    void set_geometry_type()
    {
    //  get unique geometry id
      size_t id = GetUniqueFVGeomID<TFVGeom>();
    //  check that function exists
      if(id >= m_vUpdateFunc.size() || m_vUpdateFunc[id] == NULL)
        UG_THROW("No update function registered for Geometry "<<id);
    //  set current geometry
      m_id = id;
 
    //  set sizes
      TFVGeom& geo = GeomProvider<TFVGeom>::get();
      m_numScvf = geo.num_scvf();
      m_numSh = geo.num_sh();
    //  set sizes in upwind
      if(m_spUpwind.valid()) m_spUpwind->template set_geometry_type<TFVGeom>();
    }
  
  // interface to the upwind
  
  protected:
  
    template <typename TFVGeom>
    void compute_upwind(const TFVGeom& geo,
                        const MathVector<dim> vStdVel[])
    {
      UG_NSSTAB_ASSERT(m_spUpwind.valid(), "No upwind object");
      m_spUpwind->update_upwind(geo, vStdVel);
    };
 
    template <typename TFVGeom>
    void compute_downwind(const TFVGeom& geo,
                          const MathVector<dim> vStdVel[])
    {
      UG_NSSTAB_ASSERT(m_spUpwind.valid(), "No upwind object");
      m_spUpwind->update_downwind(geo, vStdVel);
    };
    
  // the data interface (for the NS discretization)
 
  public:
  
    size_t num_ip () const {return m_numScvf;}
    
    size_t num_sh () const {return m_numSh;}
  
    const MathVector<dim>& stab_vel(size_t scvf) const
    {
      UG_NSSTAB_ASSERT(scvf < m_numScvf, "Invalid index.");
      return m_vStabVel[scvf];
    }
 
    bool vel_comp_connected() const {return m_bVelCompConnected;}
 
    number stab_shape_vel(size_t scvf, size_t compOut, size_t compIn, size_t sh) const
    {
      UG_NSSTAB_ASSERT(scvf < m_numScvf, "Invalid index.");
      UG_NSSTAB_ASSERT(compOut < dim, "Invalid index.");
      UG_NSSTAB_ASSERT(compIn < dim, "Invalid index.");
      UG_NSSTAB_ASSERT(sh < m_numSh, "Invalid index.");
      return m_vvvvStabShapeVel[scvf][compOut][compIn][sh];
    }
 
    number stab_shape_p(size_t scvf, size_t compOut, size_t sh) const
    {
      UG_NSSTAB_ASSERT(scvf < m_numScvf, "Invalid index.");
      UG_NSSTAB_ASSERT(compOut < dim, "Invalid index.");
      UG_NSSTAB_ASSERT(sh < m_numSh, "Invalid index.");
      return m_vvvStabShapePressure[scvf][compOut][sh];
    }
 
 
    void update(const FVGeometryBase* geo,
                const LocalVector& vCornerValue,
                const MathVector<dim> vStdVel[],
                const bool bStokes,
                const DataImport<number, dim>& kinVisco,
                const DataImport<number, dim>& density,
                const DataImport<MathVector<dim>, dim>* pSource,
                const LocalVector* pvCornerValueOldTime, number dt)
      {(this->*(m_vUpdateFunc[m_id]))(geo, vCornerValue, vStdVel,
                      bStokes, kinVisco, density, pSource,
                      pvCornerValueOldTime, dt);}
 
  // the data interface (for the implementation)
 
  protected:
  
    MathVector<dim>& stab_vel(size_t scvf)
    {
      UG_NSSTAB_ASSERT(scvf < m_numScvf, "Invalid index.");
      return m_vStabVel[scvf];
    }
 
    void set_vel_comp_connected(bool bVelCompConnected) {m_bVelCompConnected = bVelCompConnected;}
 
    number& stab_shape_vel(size_t scvf, size_t compOut, size_t compIn, size_t sh)
    {
      UG_NSSTAB_ASSERT(scvf < m_numScvf, "Invalid index.");
      UG_NSSTAB_ASSERT(compOut < dim, "Invalid index.");
      UG_NSSTAB_ASSERT(compIn < dim, "Invalid index.");
      UG_NSSTAB_ASSERT(sh < m_numSh, "Invalid index.");
      return m_vvvvStabShapeVel[scvf][compOut][compIn][sh];
    }
 
    number& stab_shape_p(size_t scvf, size_t compOut, size_t sh)
    {
      UG_NSSTAB_ASSERT(scvf < m_numScvf, "Invalid index.");
      UG_NSSTAB_ASSERT(compOut < dim, "Invalid index.");
      UG_NSSTAB_ASSERT(sh < m_numSh, "Invalid index.");
      return m_vvvStabShapePressure[scvf][compOut][sh];
    }
 
  // the data
  
  private:
 
    size_t m_numScvf;
 
    size_t m_numSh;
 
    MathVector<dim> m_vStabVel[maxNumSCVF];
 
    bool m_bVelCompConnected;
 
    number m_vvvvStabShapeVel[maxNumSCVF][dim][dim][maxNumSH];
 
    number m_vvvStabShapePressure[maxNumSCVF][dim][maxNumSH];
 
    int m_id;
 
    SmartPtr<INavierStokesUpwind<dim> > m_spUpwind;
    ConstSmartPtr<INavierStokesUpwind<dim> > m_spConstUpwind;
  
  // registering mechanism
 
  protected:
    typedef void (this_type::*UpdateFunc)(  const FVGeometryBase* geo,
                        const LocalVector& vCornerValue,
                        const MathVector<dim> vStdVel[],
                        const bool bStokes,
                        const DataImport<number, dim>& kinVisco,
                        const DataImport<number, dim>& density,
                        const DataImport<MathVector<dim>, dim>* pSource,
                        const LocalVector* pvCornerValueOldTime, number dt);
 
  public:
    template <typename TFVGeom, typename TAssFunc>
    void register_update_func(TAssFunc func);
 
  protected:
    std::vector<UpdateFunc> m_vUpdateFunc;
};
 
/*-------- Schneider-Raw-type stabilizations --------*/
 
template <int dim>
class INavierStokesSRFV1Stabilization
  : public INavierStokesFV1Stabilization<dim>
{
  private:
    typedef INavierStokesSRFV1Stabilization<dim> this_type;
 
  protected:
    enum DiffusionLength
    {
        RAW = 0,
        FIVEPOINT,
        COR
    };
 
  public:
    INavierStokesSRFV1Stabilization()
    {
    //  default setup
      set_diffusion_length("RAW");
    }
 
    void set_diffusion_length(std::string diffLength);
 
  protected:
  
    number diff_length_sq_inv(size_t scvf) const
    {
      UG_NSSTAB_ASSERT(scvf < this_type::num_ip(), "Invalid index.");
      return m_vDiffLengthSqInv[scvf];
    }
 
  // forward methods of Upwind Velocity
 
    number upwind_conv_length(size_t scvf) const
    {
      UG_NSSTAB_ASSERT(this_type::upwind_valid(), "No upwind object");
      return this_type::upwind()->upwind_conv_length(scvf);
    }
 
    number downwind_conv_length(size_t scvf) const
    {
      UG_NSSTAB_ASSERT(this_type::upwind_valid(), "No upwind object");
      return this_type::upwind()->downwind_conv_length(scvf);
    }
 
    number upwind_shape_sh(size_t scvf, size_t sh) const
    {
      UG_NSSTAB_ASSERT(this_type::upwind_valid(), "No upwind object");
      return this_type::upwind()->upwind_shape_sh(scvf, sh);
    }
 
    number downwind_shape_sh(size_t scvf, size_t sh) const
    {
      UG_NSSTAB_ASSERT(this_type::upwind_valid(), "No upwind object");
      return this_type::upwind()->downwind_shape_sh(scvf, sh);
    }
 
    bool non_zero_shape_ip() const
    {
      UG_NSSTAB_ASSERT(this_type::upwind_valid(), "No upwind object");
      return this_type::upwind()->non_zero_shape_ip();
    }
 
    number upwind_shape_ip(size_t scvf, size_t scvf2) const
    {
      UG_NSSTAB_ASSERT(this_type::upwind_valid(), "No upwind object");
      return this_type::upwind()->upwind_shape_ip(scvf, scvf2);
    }
 
    number downwind_shape_ip(size_t scvf, size_t scvf2) const
    {
      UG_NSSTAB_ASSERT(this_type::upwind_valid(), "No upwind object");
      return this_type::upwind()->downwind_shape_ip(scvf, scvf2);
    }
 
  // internal handling
 
  protected:
    template <typename TFVGeom>
    void compute_diff_length(const TFVGeom& geo);
 
  // data
 
  protected:
    DiffusionLength m_diffLengthType;
 
    number m_vDiffLengthSqInv[this_type::maxNumSCVF];
 
};
 
template <int dim>
SmartPtr<INavierStokesSRFV1Stabilization<dim> > CreateNavierStokesStabilization(const std::string& name);
 
// FIELDS
 
template <int TDim>
class NavierStokesFIELDSStabilization
  : public INavierStokesSRFV1Stabilization<TDim>
{
  public:
    typedef INavierStokesSRFV1Stabilization<TDim> base_type;
 
    typedef NavierStokesFIELDSStabilization<TDim> this_type;
 
    static const int dim = TDim;
 
  protected:
  //  explicitly forward some function
    using base_type::register_update_func;
    using base_type::diff_length_sq_inv;
    using base_type::stab_shape_vel;
    using base_type::stab_shape_p;
    using base_type::stab_vel;
    using base_type::set_vel_comp_connected;
 
  //  functions from upwind
    using base_type::upwind_conv_length;
    using base_type::downwind_conv_length;
    using base_type::upwind_shape_sh;
    using base_type::downwind_shape_sh;
    using base_type::non_zero_shape_ip;
    using base_type::upwind_shape_ip;
    using base_type::downwind_shape_ip;
 
  public:
    NavierStokesFIELDSStabilization()
    {
    //  vel comp not coupled
      set_vel_comp_connected(false);
 
    //  register evaluation function
      register_func();
    }
 
    template <typename TElem>
    void update(const FV1Geometry<TElem, dim>* geo,
                const LocalVector& vCornerValue,
                const MathVector<dim> vStdVel[],
                const bool bStokes,
                const DataImport<number, dim>& kinVisco,
          const DataImport<number, dim>& density,
                const DataImport<MathVector<dim>, dim>* pSource,
                const LocalVector* pvCornerValueOldTime, number dt);
 
  private:
    void register_func();
 
    template <typename TElem>
    void register_func()
    {
      typedef FV1Geometry<TElem, dim> TGeom;
      typedef void (this_type::*TFunc)(const TGeom* geo,
                       const LocalVector& vCornerValue,
                       const MathVector<dim> vStdVel[],
                       const bool bStokes,
                       const DataImport<number, dim>& kinVisco,
                       const DataImport<number, dim>& density,
                       const DataImport<MathVector<dim>, dim>* pSource,
                       const LocalVector* pvCornerValueOldTime, number dt);
 
      this->template register_update_func<TGeom, TFunc>(&this_type::template update<TElem>);
    }
};
 
 
// FLOW
 
template <int TDim>
class NavierStokesFLOWStabilization
  : public INavierStokesSRFV1Stabilization<TDim>
{
  public:
    typedef INavierStokesSRFV1Stabilization<TDim> base_type;
 
    typedef NavierStokesFLOWStabilization<TDim> this_type;
 
    static const int dim = TDim;
 
  protected:
  //  explicitly forward some function
    using base_type::register_update_func;
    using base_type::set_vel_comp_connected;
    using base_type::diff_length_sq_inv;
    using base_type::stab_shape_vel;
    using base_type::stab_shape_p;
    using base_type::stab_vel;
 
  //  functions from upwind
    using base_type::upwind_conv_length;
    using base_type::downwind_conv_length;
    using base_type::upwind_shape_sh;
    using base_type::downwind_shape_sh;
    using base_type::non_zero_shape_ip;
    using base_type::upwind_shape_ip;
    using base_type::downwind_shape_ip;
 
  public:
    NavierStokesFLOWStabilization()
    {
    //  vel comp coupled
      set_vel_comp_connected(true);
 
    //  register evaluation function
      register_func();
    }
 
    template <typename TElem>
    void update(const FV1Geometry<TElem, dim>* geo,
          const LocalVector& vCornerValue,
          const MathVector<dim> vStdVel[],
          const bool bStokes,
          const DataImport<number, dim>& kinVisco,
          const DataImport<number, dim>& density,
          const DataImport<MathVector<dim>, dim>* pSource,
          const LocalVector* pvCornerValueOldTime, number dt);
 
  private:
    void register_func();
 
    template <typename TElem>
    void register_func()
    {
      typedef FV1Geometry<TElem, dim> TGeom;
      typedef void (this_type::*TFunc)(const TGeom* geo,
                       const LocalVector& vCornerValue,
                       const MathVector<dim> vStdVel[],
                       const bool bStokes,
                       const DataImport<number, dim>& kinVisco,
                       const DataImport<number, dim>& density,
                       const DataImport<MathVector<dim>, dim>* pSource,
                       const LocalVector* pvCornerValueOldTime, number dt);
 
      this->template register_update_func<TGeom, TFunc>(&this_type::template update<TElem>);
    }
};
 
// NO STABILIZATION (Note: The discretization is then unstable!)
 
template <int TDim>
class NavierStokesFV1WithoutStabilization
  : public INavierStokesFV1Stabilization<TDim>
{
  public:
    typedef INavierStokesFV1Stabilization<TDim> base_type;
 
    typedef NavierStokesFV1WithoutStabilization<TDim> this_type;
 
    static const int dim = TDim;
 
  protected:
  //  explicitly forward some function
    using base_type::register_update_func;
    using base_type::set_vel_comp_connected;
    using base_type::stab_shape_vel;
    using base_type::stab_shape_p;
    using base_type::stab_vel;
    
  public:
    NavierStokesFV1WithoutStabilization ()
    {
    //  vel comp not interconnected
      set_vel_comp_connected(false);
 
    //  register evaluation function
      register_func();
    }
  
    template <typename TElem>
    void update(const FV1Geometry<TElem, dim>* geo,
          const LocalVector& vCornerValue,
          const MathVector<dim> vStdVel[],
          const bool bStokes,
          const DataImport<number, dim>& kinVisco,
          const DataImport<number, dim>& density,
          const DataImport<MathVector<dim>, dim>* pSource,
          const LocalVector* pvCornerValueOldTime, number dt);
 
  private:
    void register_func();
 
    template <typename TElem>
    void register_func()
    {
      typedef FV1Geometry<TElem, dim> TGeom;
      typedef void (this_type::*TFunc)(const TGeom* geo,
                       const LocalVector& vCornerValue,
                       const MathVector<dim> vStdVel[],
                       const bool bStokes,
                       const DataImport<number, dim>& kinVisco,
                       const DataImport<number, dim>& density,
                       const DataImport<MathVector<dim>, dim>* pSource,
                       const LocalVector* pvCornerValueOldTime, number dt);
 
      this->template register_update_func<TGeom, TFunc>(&this_type::template update<TElem>);
    }
};
 
} // namespace NavierStokes
} // end namespace ug
 
#endif /* __H__UG__PLUGINS__NAVIER_STOKES__FV__STABILIZATION__ */