bear_scheidegger.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2025 Gesellschaft fuer Anlagen- und Reaktorsicherheit gGmbH
// SPDX-License-Identifier: EUPL-1.2
// SPDX-FileContributor: Andreas Vogel
// SPDX-FileContributor: Goethe Universität Frankfurt
// SPDX-FileType: SOURCE
 
#ifndef __H__UG__LIB_DISC__DENSITY_DRIVEN_FLOW__BEAR_SCHEIDEGGER__
#define __H__UG__LIB_DISC__DENSITY_DRIVEN_FLOW__BEAR_SCHEIDEGGER__
 
#include "lib_disc/spatial_disc/user_data/linker/linker.h"
#include "d3f_limex.h"
 
#define LIMEX_PATCH
 
namespace ug{
 
template <int dim>
class BearScheidegger
  : public StdDataLinker< BearScheidegger<dim>, MathMatrix<dim,dim>, dim>
{
    typedef StdDataLinker< BearScheidegger<dim>, MathMatrix<dim,dim>, dim> base_type;
 
  public:
    BearScheidegger(SmartPtr<CplUserData<MathVector<dim>, dim> > vel,
                    number alphaTrans, number alphaLong) :
        m_alphaT(alphaTrans), m_alphaL(alphaLong),
      m_spVelocity(NULL), m_spDVelocity(NULL),
      m_damp_limex(1.0)
    {
      this->set_num_input(1);
      m_spVelocity = vel;
      m_spDVelocity = vel.template cast_dynamic<DependentUserData<MathVector<dim>, dim> >();
      base_type::set_input(_Q_, vel, vel);
    }
 
    BearScheidegger(SmartPtr<CplUserData<MathVector<dim>, dim> > vel,
                        number alphaTrans, number alphaLong, number damp) :
      m_alphaT(alphaTrans), m_alphaL(alphaLong),
      m_spVelocity(NULL), m_spDVelocity(NULL),
      m_damp_limex(damp)
    {
      this->set_num_input(1);
      m_spVelocity = vel;
      m_spDVelocity = vel.template cast_dynamic<DependentUserData<MathVector<dim>, dim> >();
      base_type::set_input(_Q_, vel, vel);
    }
    
    double alphaT () const {return m_alphaT;}
    double alphaL () const {return m_alphaL;}
    SmartPtr<CplUserData<MathVector<dim>, dim> > velocity () const {return m_spVelocity;}
 
    inline void evaluate (MathMatrix<dim,dim>& D,
                          const MathVector<dim>& globIP,
                          number time, int si) const
    {
      MathVector<dim> q;
      (*m_spVelocity)(q, globIP, time, si);
 
      const number qabs = VecTwoNorm(q);
      if(qabs == 0.0){
        D = 0.0;
      }
      else{
        const number alpha = (m_alphaL - m_alphaT)/qabs;
 
        for(size_t d1 = 0; d1 < dim; ++d1)
          for(size_t d2 = 0; d2 < dim; ++d2)
            D(d1,d2) = alpha*q[d1]*q[d2];
 
        for(size_t d = 0; d < dim; ++d)
          D(d,d) += qabs*m_alphaT;
      }
    }
 
    template <int refDim>
    inline void evaluate(MathMatrix<dim,dim> vD[],
                         const MathVector<dim> vGlobIP[],
                         number time, int si,
                         GridObject* elem,
                         const MathVector<dim> vCornerCoords[],
                         const MathVector<refDim> vLocIP[],
                         const size_t nip,
                         LocalVector* u,
                         const MathMatrix<refDim, dim>* vJT = NULL) const
    {
      std::vector<MathVector<dim> > vQ(nip);
      (*m_spVelocity)(&vQ[0], vGlobIP, time, si,
              elem, vCornerCoords, vLocIP, nip, u, vJT);
 
      for(size_t ip = 0; ip < nip; ++ip)
      {
        MathMatrix<dim,dim>& D = vD[ip];
        MathVector<dim>& q = vQ[ip];
 
        const number qabs = VecTwoNorm(q);
 
        if(qabs == 0.0){
          D = 0.0;
        }
        else{
          const number alpha = (m_alphaL - m_alphaT)/qabs;
 
          for(size_t d1 = 0; d1 < dim; ++d1)
            for(size_t d2 = 0; d2 < dim; ++d2)
              D(d1,d2) = alpha*q[d1]*q[d2];
 
          for(size_t d = 0; d < dim; ++d)
            D(d,d) += qabs*m_alphaT;
        }
      }
    }
 
    template <int refDim>
    void eval_and_deriv(MathMatrix<dim,dim> vD[],
                        const MathVector<dim> vGlobIP[],
                        number time, int si,
                        GridObject* elem,
                        const MathVector<dim> vCornerCoords[],
                        const MathVector<refDim> vLocIP[],
                        const size_t nip,
                        LocalVector* u,
                        bool bDeriv,
                        int subset,
                        std::vector<std::vector<MathMatrix<dim,dim> > > vvvDeriv[],
                        const MathMatrix<refDim, dim>* vJT = NULL) const
    {
      const bool noDeriv =  !bDeriv ||
                  !m_spDVelocity.valid() ||
                  m_spDVelocity->zero_derivative();
 
      // a check
      if(!noDeriv)  {
        for(size_t fct = 0; fct < m_spDVelocity->num_fct(); ++fct)
          if(this->input_common_fct(_Q_, fct) != fct)
            UG_THROW("BearScheidegger: commonFct <-> fct should be equal");
      }
 
      const MathVector<dim>* vQ = m_spVelocity->values(subset);
      for(size_t ip = 0; ip < nip; ++ip)
      {
        MathMatrix<dim,dim>& D = vD[ip];
        const MathVector<dim>& q = vQ[ip];
 
        const number qabs = VecTwoNorm(q);
        if(qabs == 0.0){
          D = 0.0;
          if(noDeriv) continue;
          for(size_t fct = 0; fct < m_spDVelocity->num_fct(); ++fct)
            for(size_t sh = 0; sh < this->num_sh(fct); ++sh)
              vvvDeriv[ip][fct][sh] = 0.0;
          continue;
        }
 
        const number alpha = (m_alphaL - m_alphaT)/qabs;
 
        for(size_t d1 = 0; d1 < dim; ++d1)
          for(size_t d2 = 0; d2 < dim; ++d2)
            D(d1,d2) = alpha*q[d1]*q[d2];
 
        for(size_t d = 0; d < dim; ++d)
          D(d,d) += qabs*m_alphaT;
 
      // Derivative   
#ifndef DF_LIMEX_PATCH
        if(noDeriv) continue;
#else
        if(noDeriv || m_damp_limex==0.0) continue;              
#endif        
 
        for(size_t fct = 0; fct < m_spDVelocity->num_fct(); ++fct)
        {
          const MathVector<dim>* vDVelocity = m_spDVelocity->deriv(subset, ip, fct);
 
          for(size_t sh = 0; sh < this->num_sh(fct); ++sh)
          {
            const MathVector<dim>& qLin = vDVelocity[sh];
            const number qabslin = VecDot(q, qLin) / qabs;
            const number qabslinOverQabs = qabslin / qabs;
 
            for(size_t d1 = 0; d1 < dim; ++d1)
              for(size_t d2 = 0; d2 < dim; ++d2)
                vvvDeriv[ip][fct][sh](d1,d2) =
                    alpha*(qLin[d1]*q[d2] + q[d1]*qLin[d2]
#ifdef DF_LIMEX_PATCH
                          - q[d1]*q[d2] * qabslinOverQabs)*m_damp_limex;
#else
                        - q[d1]*q[d2] * qabslinOverQabs);
#endif
 
            for(size_t d = 0; d < dim; ++d)
#ifdef DF_LIMEX_PATCH
              vvvDeriv[ip][fct][sh](d,d) += qabslin*m_alphaT*m_damp_limex;
#else
              vvvDeriv[ip][fct][sh](d,d) += qabslin*m_alphaT;
#endif
          }
        }
      }
    }
#ifdef DF_LIMEX_PATCH
    void set_damp(number x) {m_damp_limex = x;}
#endif
 
  protected:
    const number m_alphaT;
    const number m_alphaL;
 
    static const size_t _Q_ = 0;
 
    SmartPtr<CplUserData<MathVector<dim>, dim> > m_spVelocity;
    SmartPtr<DependentUserData<MathVector<dim>, dim> > m_spDVelocity;
#ifdef DF_LIMEX_PATCH
    number m_damp_limex;
#endif
};
 
} // end namespace ug
 
#endif /* __H__UG__LIB_DISC__DENSITY_DRIVEN_FLOW__BEAR_SCHEIDEGGER__ */