fract_ddf_fv1.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: Dmitry Logashenko, Sabine Stichel
 * SPDX-FileContributor: Goethe Universität Frankfurt
 * SPDX-FileType: SOURCE
 * 
 * This file is part of d3f++.
 * d3f++ is an extension for UG4. Licensing information and citation requirements of UG4 are provided in LICENSES/UG4-LGPL_2.1
 */
 
/*
 * fract_ddf_fv1.h - discretization of the density driven flow in porous media
 * with fractures.
 *
 * Based on the modules by A. Vogel.
 */
 
#ifndef __H__UG__LIB_DISC__DENSITY_DRIVEN_FLOW__FRACT_DDF_FV1_
#define __H__UG__LIB_DISC__DENSITY_DRIVEN_FLOW__FRACT_DDF_FV1_
 
// ug4 headers
 
#include "common/common.h"
#include "lib_grid/lg_base.h"
#include "lib_grid/algorithms/deg_layer_mngr.h"
 
#include "lib_disc/spatial_disc/elem_disc/elem_disc_interface.h"
#include "lib_disc/spatial_disc/user_data/data_export.h"
#include "lib_disc/spatial_disc/user_data/data_import.h"
#include "lib_disc/spatial_disc/user_data/const_user_data.h"
#include "lib_disc/spatial_disc/disc_util/conv_shape_interface.h"
#include "lib_disc/spatial_disc/disc_util/consistent_gravity.h"
 
#ifdef UG_FOR_LUA
#include "bindings/lua/lua_user_data.h"
#endif
 
// d3f headers
#include "../../extravol.h"
#include "../../d3f_sss.h"
 
namespace ug{
namespace d3f{
 
 
 
template <typename TDomain>
class FractDDF_FV1
  : public IElemDisc<TDomain>
{
private:
  typedef IElemDisc<TDomain> base_type;
 
  typedef FractDDF_FV1<TDomain> this_type;
 
  typedef typename base_type::domain_type domain_type;
 
  typedef typename base_type::position_type position_type;
  
  static const int dim = base_type::dim;
 
  static const int low_dim = dim - 1;
  
  typedef DimFV1Geometry<low_dim, dim> TFractFVGeom;
  
  typedef typename domain_traits<dim>::DimElemList AssembleElemList;
  
  typedef DegeneratedLayerManager<dim> fract_manager_type;
  
  typedef typename fract_manager_type::side_type side_type;
  
  static const size_t maxFractSideCorners = fract_manager_type::maxLayerSideCorners;
  
  static const size_t _C_ = 0;
 
  static const size_t _P_ = 1;
  
public:
 
  FractDDF_FV1
  (
    const char* functions, 
    const char* subsets 
  );
  FractDDF_FV1
  (
    const std::vector<std::string>& vFct,
      const std::vector<std::string>& vSubset
  );
  void init();
  
//---- Settings: ----
public:
  
  void set_fract_manager (SmartPtr<fract_manager_type> fract_manager)
  {
    m_spFractManager = fract_manager;
  }
 
  void set_boussinesq_transport (bool bUse) {m_BoussinesqTransport = bUse;}
 
  void set_boussinesq_flow (bool bUse) {m_BoussinesqFlow = bUse;}
  
  void set_boussinesq (bool bUse) {m_BoussinesqFlow = m_BoussinesqTransport = bUse;};
 
  void set_upwind (SmartPtr<IConvectionShapes<dim> > upwind)
  {
    m_spUpwind = upwind;
  }
 
  void set_density(SmartPtr<CplUserData<number,dim> > data)
  {
  //  remove old data
    SmartPtr<ICplUserData<dim> > oldData = m_imDensityCo.data();
    if (oldData.valid())
    {
      m_exFractDarcyVel->remove_needed_data(oldData);
      m_exOrthoFractDarcyVel->remove_needed_data(oldData);
    }
    oldData = m_imDensityIP.data();
    if (oldData.valid())
    {
      m_exFractDarcyVel->remove_needed_data(oldData);
      m_exOrthoFractDarcyVel->remove_needed_data(oldData);
    }
 
  //  connect to import
    m_imDensityCo.set_data(data);
    m_imDensityIP.set_data(data);
    if (m_spVolStabData.valid())
      m_imOldDensityCo.set_data(data);
    else
      m_imOldDensityCo.set_data(SPNULL);
  
  //  reset the data in the old density
    set_volume_stabilization(m_spVolStabData);
 
  //  darcy velocity depends on density
    m_exFractDarcyVel->add_needed_data(data);
    m_exOrthoFractDarcyVel->add_needed_data(data);
  }
  void set_density(number val)
  {
    set_density(make_sp(new ConstUserNumber<dim>(val)));
  }
  void set_viscosity(SmartPtr<CplUserData<number, dim> > user)
  {
  //  remove old data
    SmartPtr<ICplUserData<dim> > oldData = m_imViscosityCo.data();
    if (oldData.valid())
    {
      m_exFractDarcyVel->remove_needed_data(oldData);
      m_exOrthoFractDarcyVel->remove_needed_data(oldData);
    }
    oldData = m_imViscosityIP.data();
    if (oldData.valid())
    {
      m_exFractDarcyVel->remove_needed_data(oldData);
      m_exOrthoFractDarcyVel->remove_needed_data(oldData);
    }
 
  //  connect to import
    m_imViscosityCo.set_data(user);
    m_imViscosityIP.set_data(user);
 
  //  darcy velocity depends on viscosity
    m_exFractDarcyVel->add_needed_data(user);
    m_exOrthoFractDarcyVel->add_needed_data(user);
 
  }
  void set_viscosity(number val)
  {
    set_viscosity(make_sp(new ConstUserNumber<dim>(val)));
  }
  void set_gravity(SmartPtr<CplUserData<MathVector<dim>, dim> > user)
  {
    m_imConstGravity.set_data(user);
  }
  void set_gravity(const std::vector<number>& vGravity)
  {
    set_gravity(SmartPtr<ConstUserVector<dim> >(new ConstUserVector<dim>(vGravity)));
  }
  
  void set_porosity(SmartPtr<CplUserData<number, dim> > user)
  {
    m_imPorosity.set_data(user);
  }
  void set_porosity(number val)
  {
    set_porosity(make_sp(new ConstUserNumber<dim>(val)));
  }
#ifdef UG_FOR_LUA
  void set_porosity(const char* fctName)
  {
    set_porosity(LuaUserDataFactory<number, dim>::create(fctName));
  }
  void set_porosity(LuaFunctionHandle fct)
  {
    set_porosity(make_sp(new LuaUserData<number, dim>(fct)));
  }
#endif
  
 
  void set_permeability(SmartPtr<CplUserData<MathMatrix<dim, dim>, dim> > user)
  {
    m_imPermeability.set_data(user);
  }
  void set_permeability(number val)
  {
    set_permeability(make_sp(new ConstUserMatrix<dim>(val)));
  }
#ifdef UG_FOR_LUA
  void set_permeability(const char* fctName)
  {
    set_permeability(LuaUserDataFactory<MathMatrix<dim,dim>, dim>::create(fctName));
  }
  void set_permeability(LuaFunctionHandle fct)
  {
    set_permeability(make_sp(new LuaUserData<MathMatrix<dim,dim>, dim>(fct)));
  }
#endif
  
  void set_molecular_diffusion(SmartPtr<CplUserData<MathMatrix<dim, dim>, dim> > user)
  {
    m_imDiffusion.set_data(user);
  }
  void set_molecular_diffusion(number val)
  {
    set_molecular_diffusion(make_sp(new ConstUserMatrix<dim>(val)));
  }
#ifdef UG_FOR_LUA
  void set_molecular_diffusion(const char* fctName)
  {
    set_molecular_diffusion(LuaUserDataFactory<MathMatrix<dim,dim>, dim>::create(fctName));
  }
  void set_molecular_diffusion(LuaFunctionHandle fct)
  {
    set_molecular_diffusion(make_sp(new LuaUserData<MathMatrix<dim,dim>, dim>(fct)));
  }
#endif
  
 
  void set_aperture(SmartPtr<CplUserData<number, dim> > user)
  {
    m_imAperture.set_data(user);
  }
  void set_aperture(number val)
  {
    set_aperture(make_sp(new ConstUserNumber<dim>(val)));
  }
#ifdef UG_FOR_LUA
  void set_aperture(const char* fctName)
  {
    set_aperture(LuaUserDataFactory<number, dim>::create(fctName));
  }
  void set_aperture(LuaFunctionHandle fct)
  {
    set_aperture(make_sp(new LuaUserData<number, dim>(fct)));
  }
#endif
  
  void set_fract_permeability(SmartPtr<CplUserData<number, dim> > user)
  {
    m_imFractPermeability.set_data(user);
  }
  void set_fract_permeability(number val)
  {
    set_fract_permeability(make_sp(new ConstUserNumber<dim>(val)));
  }
#ifdef UG_FOR_LUA
  void set_fract_permeability(const char* fctName)
  {
    set_fract_permeability(LuaUserDataFactory<number, dim>::create(fctName));
  }
  void set_fract_permeability(LuaFunctionHandle fct)
  {
    set_fract_permeability(make_sp(new LuaUserData<number, dim>(fct)));
  }
#endif
  
  void set_orthofract_permeability(SmartPtr<CplUserData<number, dim> > user)
  {
    m_imOrthoPermeability.set_data(user);
  }
  void set_orthofract_permeability(number val)
  {
    set_orthofract_permeability(make_sp(new ConstUserNumber<dim>(val)));
  }
#ifdef UG_FOR_LUA
  void set_orthofract_permeability(const char* fctName)
  {
    set_orthofract_permeability(LuaUserDataFactory<number, dim>::create(fctName));
  }
  void set_orthofract_permeability (LuaFunctionHandle fct)
  {
    set_orthofract_permeability (make_sp(new LuaUserData<number, dim> (fct)));
  }
#endif
  
  void set_fract_molecular_diffusion(SmartPtr<CplUserData<number, dim> > user)
  {
    m_imFractDiffusion.set_data(user);
  }
  void set_fract_molecular_diffusion(number val)
  {
    set_fract_molecular_diffusion(make_sp(new ConstUserNumber<dim>(val)));
  }
#ifdef UG_FOR_LUA
  void set_fract_molecular_diffusion(const char* fctName)
  {
    set_fract_molecular_diffusion(LuaUserDataFactory<number, dim>::create(fctName));
  }
  void set_fract_molecular_diffusion (LuaFunctionHandle fct)
  {
    set_fract_molecular_diffusion (make_sp(new LuaUserData<number, dim> (fct)));
  }
#endif
  
  void set_orthofract_molecular_diffusion(SmartPtr<CplUserData<number, dim> > user)
  {
    m_imOrthoDiffusion.set_data(user);
  }
  void set_orthofract_molecular_diffusion(number val)
  {
    set_orthofract_molecular_diffusion(make_sp(new ConstUserNumber<dim>(val)));
  }
#ifdef UG_FOR_LUA
  void set_orthofract_molecular_diffusion(const char* fctName)
  {
    set_orthofract_molecular_diffusion(LuaUserDataFactory<number, dim>::create(fctName));
  }
  void set_orthofract_molecular_diffusion (LuaFunctionHandle fct)
  {
    set_orthofract_molecular_diffusion (make_sp(new LuaUserData<number, dim> (fct)));
  }
#endif
  
  void set_sss_manager(SmartPtr<DDSingularSourcesAndSinks<dim> > sss_mngr) {m_sss_mngr = sss_mngr;}
 
  SmartPtr<DDSingularSourcesAndSinks<dim> > sss_manager() {return m_sss_mngr;}
 
  void set_volume_stabilization(SmartPtr<FlowVolStabData<TDomain> > volStabData)
  {
    m_imOldDensityCo.set_data (volStabData.valid()? m_imDensityCo.user_data() : SPNULL);
    m_spVolStabData = volStabData;
  }
 
  SmartPtr<CplUserData<number, dim> > brine() {return m_exBrine;}
 
  SmartPtr<CplUserData<number, dim> > pressure() {return m_exPressure;}
 
  SmartPtr<CplUserData<MathVector<dim>, dim> > DarcyVel() {return m_exFractDarcyVel;}
 
  SmartPtr<CplUserData<number, dim> > OrthoDarcyVel() {return m_exOrthoFractDarcyVel;}
 
//---- Local discretization interface: ----
protected:
  
  virtual void prepare_setting
  (
    const std::vector<LFEID> & vLfeID,
    bool bNonRegular
  );
 
//---- Assembling functions: ----
protected:
  
  virtual void prep_assemble_loop();
 
  template <typename TElem>
  void prepare_element_loop(ReferenceObjectID roid, int si);
 
  template <typename TElem>
  void prepare_element(const LocalVector& u, GridObject* elem, ReferenceObjectID roid, const position_type vCornerCoords[]);
 
  template <typename TElem>
  void finish_element_loop();
 
  template <typename TElem>
  void ass_JA_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  void ass_JM_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  void ass_dA_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  void ass_dM_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  void ass_rhs_elem(LocalVector& d, GridObject* elem, const position_type vCornerCoords[]);
  
  template <typename TElem>
  void ass_sss_dA_elem(LocalVector& d, const LocalVector& u, TElem* pElem, size_t co, number intensity, number concentration);
  
  template <typename TElem>
  void ass_sss_JA_elem(LocalMatrix& J, const LocalVector& u, TElem* pElem, size_t co, number intensity, number concentration);
  
  virtual bool requests_local_time_series()
  {
    return m_spVolStabData.valid (); // the stabilization requires the old solution
  }
    
//---- Assembling functions for bulk medium: ----
private:
 
  template <typename TElem>
  inline void bulk_prepare_element_loop(ReferenceObjectID roid, int si);
 
  template <typename TElem>
  inline void bulk_prepare_element(const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  inline void bulk_ass_JA_elem(LocalMatrix& J, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  inline void bulk_ass_JM_elem(LocalMatrix& J, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  inline void bulk_ass_dA_elem(LocalVector& d, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  inline void bulk_ass_dM_elem(LocalVector& d, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  inline void bulk_ass_rhs_elem(LocalVector& d, TElem* elem, const position_type vCornerCoords[]);
  
//---- Assembling functions for fractures: ----
private:
 
  template <typename TElem>
  inline void fract_prepare_element_loop(ReferenceObjectID roid, int si);
 
  template <typename TElem>
  inline void fract_prepare_element(const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  inline void fract_ass_JA_elem(LocalMatrix& J, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  inline void fract_bulk_ass_JA_elem(LocalMatrix& J, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  inline void fract_ass_JM_elem(LocalMatrix& J, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  inline void fract_ass_dA_elem(LocalVector& d, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  inline void fract_bulk_ass_dA_elem(LocalVector& d, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  inline void fract_ass_dM_elem(LocalVector& d, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
 
  template <typename TElem>
  inline void fract_ass_rhs_elem(LocalVector& d, TElem* elem, const position_type vCornerCoords[]);
 
//---- Computation of the Darcy velocity: ----
protected:
  
  template<typename TFVGeom, typename TConsGravity>
  inline void compute_ip_Darcy_velocity
  (
    MathVector<dim>& Vel, 
    size_t ip, 
    const TFVGeom& bulk_geo, 
    TConsGravity& ConsGravityMethod, 
    MathVector<TFVGeom::dim> vConsGravity[], 
    number vPressure[], 
    number Viscosity 
  );
 
  template<typename TFVGeom, typename TConsGravity, size_t maxCorners>
  void compute_J_ip_Darcy_velocity
  (
    MathVector<dim>& Vel, 
    MathVector<dim> Vel_c[], 
    MathVector<dim> Vel_p[], 
    size_t ip, 
    const TFVGeom& geo, 
    TConsGravity& ConsGravityMethod, 
    MathVector<TFVGeom::dim> vConsGravity_c[][maxCorners], 
    number vPressure[], 
    number Viscosity, 
    number Viscosity_c[] 
  );
  
//---- Export parameters: ----
private:
  template <typename TElem>
  void ex_brine
  (
    number vValue[],
    const MathVector<dim> vGlobIP[],
    number time, int si,
    const LocalVector& u,
    GridObject* elem,
    const MathVector<dim> vCornerCoords[],
    const MathVector<dim> vLocIP[],
    const size_t nip,
    bool bDeriv,
    std::vector<std::vector<number> > vvvDeriv[]
  );
  
  template <typename TElem>
  void ex_pressure
  (
    number vValue[],
    const MathVector<dim> vGlobIP[],
    number time, int si,
    const LocalVector& u,
    GridObject* elem,
    const MathVector<dim> vCornerCoords[],
    const MathVector<dim> vLocIP[],
    const size_t nip,
    bool bDeriv,
    std::vector<std::vector<number> > vvvDeriv[]
  );
  
  template <typename TElem>
  void ex_darcy_fract
  (
    MathVector<dim> vValue[],
    const MathVector<dim> vGlobIP[],
    number time, int si,
    const LocalVector& u,
    GridObject* elem,
    const MathVector<dim> vCornerCoords[],
    const MathVector<dim> vLocIP[],
    const size_t nip,
    bool bDeriv,
    std::vector<std::vector<MathVector<dim> > > vvvDeriv[]
  );
  
  template <typename TElem>
  void ex_darcy_ortho_fract
  (
    number vValue[],
    const MathVector<dim> vGlobIP[],
    number time, int si,
    const LocalVector& u,
    GridObject* elem,
    const MathVector<dim> vCornerCoords[],
    const MathVector<dim> vLocIP[],
    const size_t nip,
    bool bDeriv,
    std::vector<std::vector<number> > vvvDeriv[]
  );
  
//---- Registration of the template functions: ----
private:
  
  struct RegisterLocalDiscr {
      RegisterLocalDiscr(this_type * pThis) : m_pThis(pThis){}
      this_type * m_pThis;
      template< typename TElem > void operator() (TElem &)
      {m_pThis->register_loc_discr_func<TElem> ();}
  };
 
  template <typename TElem>
  void register_loc_discr_func ();
 
private:
  
//  Global parameters (for all subdomains):
  DataImport<number, dim> m_imDensityCo; 
  DataImport<number, dim> m_imDensityIP; 
  DataImport<number, dim> m_imViscosityCo; 
  DataImport<number, dim> m_imViscosityIP; 
  DataImport<MathVector<dim>, dim> m_imConstGravity; 
  DataImport<number, dim> m_imOldDensityCo; 
  
//  Specific parameters of subdomains:
  DataImport<number, dim> m_imPorosity; 
  
//  Specific parameters for full-dimensional subdomains:
  DataImport<MathMatrix<dim, dim>, dim> m_imPermeability; 
  DataImport<MathMatrix<dim, dim>, dim> m_imDiffusion; 
  
//  Specific parameters for low-dimensional subdomains (fractures):
  DataImport<number, dim> m_imAperture; 
  DataImport<number, dim> m_imFractPermeability; 
  DataImport<number, dim> m_imOrthoPermeability; 
  DataImport<number, dim> m_imFractDiffusion; 
  DataImport<number, dim> m_imOrthoDiffusion; 
  
  SmartPtr<DDSingularSourcesAndSinks<dim> > m_sss_mngr;
  
//  Parameters of the numerical methods of the discretization:
  SmartPtr<fract_manager_type> m_spFractManager; 
  bool m_BoussinesqFlow; 
  bool m_BoussinesqTransport; 
  SmartPtr<IConvectionShapes<dim> > m_spUpwind; 
  
  DataImport<MathVector<dim>, dim> m_imFractDarcyVelIP;
 
  SmartPtr<DataExport<number, dim> > m_exBrine;
 
  SmartPtr<DataExport<number, dim> > m_exPressure;
 
  SmartPtr<DataExport<MathVector<dim>, dim> > m_exFractDarcyVel;
  
  SmartPtr<DataExport<number, dim> > m_exOrthoFractDarcyVel;
 
  SmartPtr<FlowVolStabData<TDomain> > m_spVolStabData;
  bool m_bVolStabDataActive;
 
  MathVector<dim> m_Gravity; 
  
//  Temporary data used in the assembling
  bool m_isFracture; 
  
  TFractFVGeom * m_pFractGeo; 
  size_t m_numFractCo; 
  side_type * m_innerFractSide; 
  side_type * m_outerFractSide; 
  size_t m_innerFractSideIdx; 
  size_t m_outerFractSideIdx; 
  size_t m_innerSideCo[maxFractSideCorners]; 
  size_t m_outerSideCo[maxFractSideCorners]; 
  size_t m_assCo [2 * maxFractSideCorners]; 
  number m_orthGravity; 
  MathVector<dim> m_unitOuterNormal; 
  MathVector<dim> m_elem_loc_coe; 
  MathVector<dim> m_elem_loc_scvf [TFractFVGeom::maxNumSCVF]; 
};
 
 
} // namespace d3f
} // end namespace ug
 
#include "fract_ddf_fv1_impl.h"
 
#endif /* __H__UG__LIB_DISC__DENSITY_DRIVEN_FLOW__FRACT_DDF_FV1_ */
 
/* End of File */