rivers.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
 * 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
 */
 
/*
 * Contribution of the surface rivers represented as sets of 1d segments in 3d.
 */
#ifndef __H__UG__PLUGINS__D3F__RIVERS__
#define __H__UG__PLUGINS__D3F__RIVERS__
 
#include <vector>
#include <limits>
#include <string>
 
// ug4 headers
#include "common/common.h"
#include "common/math/misc/math_constants.h"
#include "lib_grid/grid_objects/grid_dim_traits.h"
#include "lib_disc/function_spaces/grid_function.h"
 
// d3f headers
#include "level_set_pos.h"
#include "segment_network.h"
 
namespace ug {
namespace d3f {
 
template <int dim>
inline number cone_mollifier
(
  number d, 
  number smooth_len 
);
 
template <>
inline number cone_mollifier<1>
(
  number d, 
  number smooth_len 
)
{
  d /= smooth_len;
  if (d >= 1) return 0;
  return (1 - d) / (smooth_len / 2);
}
 
template <>
inline number cone_mollifier<2>
(
  number d, 
  number smooth_len 
)
{
  d /= smooth_len;
  if (d >= 1) return 0;
  return (1 - d) / (PI * smooth_len * smooth_len / 3);
}
 
template <typename TGridFunction>
class VertexLowDimRecharge
{
public:
 
  typedef TGridFunction grid_func_type;
  
  typedef typename TGridFunction::domain_type domain_type;
  
  typedef typename TGridFunction::algebra_type algebra_type;
  
  typedef typename domain_type::grid_type grid_type;
  
  static const int dim = TGridFunction::dim;
  
  static const int mol_dim = dim - 1;
  
  typedef typename grid_dim_traits<dim>::side_type g_surf_elem_t;
  
  typedef typename geometry_traits<g_surf_elem_t>::iterator g_surf_iter_t;
  
  typedef lg_ntree<dim-1, dim, g_surf_elem_t> top_tracer_tree_t;
  typedef RayElemIntersectionRecord<g_surf_elem_t*> top_intersection_record_t;
  
public:
 
  VertexLowDimRecharge
  (
    number smoothLen, 
    SmartPtr<grid_func_type> spLSF 
  )
  : m_density (998.2), m_viscosity (1e-3),
    m_perm (1), m_e (1), m_gravity (- 9.81),
    m_smoothLen (smoothLen),
    m_restrict_depth (false),
    m_max_depth_diff (0),
    m_no_inflow (false),
    m_lsPosZ (spLSF),
    m_is_constant (false), m_constant_flux (0)
  {}
  
  number smooth_len () const {return m_smoothLen;}
  
  void set_permeability (number k) {m_perm = k;}
  
  void set_bottom_thickness (number e) {m_e = e;}
  
  void set_density (number rho) {m_density = rho;}
  
  void set_viscosity (number mu) {m_viscosity = mu;}
  
  void set_gravity (number g) {m_gravity = g;}
  
  void restrict_depth (number max_depth_diff)
  {
    m_restrict_depth = true;
    m_max_depth_diff = max_depth_diff;
  }
  
  bool restrict_depth () const {return m_restrict_depth;}
  
  number max_depth_diff () const {return m_max_depth_diff;}
  
  void set_no_inflow (bool b) {m_no_inflow = b;}
  
  bool no_inflow () const {return m_no_inflow;}
  
  void set_relative_to (const char* subsetNames, number top_tolerance, int top_grid_level);
  
  void set_relative_to (const char* subsetNames, number top_tolerance)
  {
    set_relative_to (subsetNames, top_tolerance, -1);
  }
  
  void set_relative_to (const char* subsetNames)
  {
    set_relative_to (subsetNames, 1e-6, -1);
  }
  
  void set_constant_flux (number flux)
  {
    m_is_constant = true;
    m_constant_flux = flux;
  }
  
  void compute
  (
    SmartPtr<grid_func_type> spRecharge 
  );
  
  void simple_compute
  (SmartPtr<grid_func_type> spRecharge); 
 
 
protected:
 
  void prepare_ls_height () {m_lsPosZ.reinit ();}
  
  number get_ls_height_at (MathVector<dim-1>& x)
  {
    number fs_depth;
    if (! m_lsPosZ.get_height_at (x, fs_depth))
      fs_depth = - std::numeric_limits<number>::max (); // the segment is not covered by the free surface
    return fs_depth;
  }
  
  bool ls_depth_is_relative () {return m_spTopTracerTree.valid ();}
  
  size_t min_top_depth (MathVector<dim-1>& x, number& min_depth)
  {
    MathVector<dim> X;
    for (size_t i = 0; i < dim - 1; i++) X[i] = x[i];
    X[dim-1] = 0;
    size_t n = get_top_z (X);
    if (n != 0)
    {
    // get the minimum z:
      min_depth = m_topTracePoints[0][dim-1];
      for (size_t i = 1; i < m_topTracePoints.size (); i++)
        if (min_depth > m_topTracePoints[i][dim-1])
          min_depth = m_topTracePoints[i][dim-1];
    }
    return n;
  }
  
  number object_recharge
  (
    number h_s, 
    number h_fs, 
    number w 
  ) const
  {
    if (! m_is_constant)
      return (- m_perm * m_density * m_gravity * (h_s - h_fs) / m_e / m_viscosity) * w;
    else
      return m_constant_flux;
  }
  
  size_t get_top_z (const MathVector<dim> & over);
  
  virtual void compute_depths () = 0;
  
  virtual number smoothed_recharge_at
  (
    MathVector<dim> & x 
  ) const = 0;
  
  virtual number simple_smoothed_recharge_at
  (
    MathVector<dim> & x 
  ) const = 0;
 
protected:
  
  number m_density;
  
  number m_viscosity;
  
  number m_perm;
  
  number m_e;
  
  number m_gravity;
  
  number m_smoothLen;
  
  bool m_restrict_depth;
  
  number m_max_depth_diff;
  
  bool m_no_inflow;
  
  LSPositionZ<grid_func_type> m_lsPosZ;
  
  SmartPtr<top_tracer_tree_t> m_spTopTracerTree;
  
  std::vector<top_intersection_record_t> m_topIntersectionRecords;
  
  std::vector<MathVector<dim> > m_topTracePoints;
  
  number m_top_tolerance;
  
//--- Some debugging facilities
  
  bool m_is_constant;
  number m_constant_flux;
};
 
template <typename TGridFunction>
class VertexWellRecharge
: public VertexLowDimRecharge<TGridFunction>
{
  typedef VertexWellRecharge<TGridFunction> this_type;
 
public:
 
  typedef TGridFunction grid_func_type;
  
  typedef typename TGridFunction::domain_type domain_type;
  
  typedef typename TGridFunction::algebra_type algebra_type;
  
  typedef typename domain_type::grid_type grid_type;
  
  static const int dim = TGridFunction::dim;
  
private:
  
  struct well_recharge_data_t
  {
    MathVector<dim-1> x; 
    number well_depth_spec; 
    
    number fs_depth; 
    number well_depth; 
    
    well_recharge_data_t
    (
      MathVector<dim-1> the_x,
      number the_well_depth_spec
    )
    : x (the_x), well_depth_spec (the_well_depth_spec)
    {}
  };
 
public:
 
  VertexWellRecharge
  (
    number smoothLen, 
    SmartPtr<grid_func_type> spLSF 
  )
  : VertexLowDimRecharge<TGridFunction> (smoothLen, spLSF),
    m_width (1)
  {}
  
  virtual ~VertexWellRecharge () {}
  
  void set_width (number riverWidth) {m_width = riverWidth;}
  
  void add_well (MathVector<dim-1>& x, number depth)
  {
    m_wells.push_back (well_recharge_data_t (x, depth));
  }
  
  void add_well (std::vector<number> coord, number depth)
  {
    if (coord.size () != (size_t) (dim - 1))
      UG_THROW ("VertexWellRecharge::add_well: Wrong number of coordinates specified.\n");
    MathVector<dim-1> x;
    for (size_t i = 0; i < dim-1; i++) x[i] = coord[i];
    m_wells.push_back (well_recharge_data_t (x, depth));
  }
  
private:
  
  virtual void compute_depths ();
  
  virtual number smoothed_recharge_at
  (
    MathVector<dim> & x 
  ) const;
  
  virtual number simple_smoothed_recharge_at
  (
    MathVector<dim> & x 
  ) const;
 
private:
 
  number m_width;
  
  std::vector<well_recharge_data_t> m_wells;
};
 
template <typename TGridFunction>
class VertexRiverRecharge
: public VertexLowDimRecharge<TGridFunction>
{
  typedef VertexRiverRecharge<TGridFunction> this_type;
 
public:
 
  typedef TGridFunction grid_func_type;
  
  typedef typename TGridFunction::domain_type domain_type;
  
  typedef typename TGridFunction::algebra_type algebra_type;
  
  typedef typename domain_type::grid_type grid_type;
  
  static const int dim = TGridFunction::dim;
  
private:
 
  typedef ANumber depth_attachment_type;
 
  typedef Grid::EdgeAttachmentAccessor<depth_attachment_type> depth_accessor_type;
  
public:
  
  VertexRiverRecharge
  (
    SmartPtr<SegmentNetwork> spRiverNetwork, 
    number smoothLen, 
    SmartPtr<grid_func_type> spLSF 
  )
  : VertexLowDimRecharge<TGridFunction> (smoothLen, spLSF),
    m_spRiverNetwork (spRiverNetwork),
    m_bNodalFSDepth (false),
    m_bCheckFSCovering (true), m_bCheckTopCovering (false),
    m_bPrintTotalRecharge (false)
  {
    Grid & riverGrid = m_spRiverNetwork->grid ();
    
    riverGrid.attach_to_edges (m_aFSDepth);
    m_aaFSDepth.access (riverGrid, m_aFSDepth);
    
    riverGrid.attach_to_edges (m_aRiverDepth);
    m_aaRiverDepth.access (riverGrid, m_aRiverDepth);
    
    riverGrid.attach_to_edges (m_aRiverWidth);
    m_aaRiverWidth.access (riverGrid, m_aRiverWidth);
    
    set_river_width (1); // a "very default" value
  }
  
  virtual ~VertexRiverRecharge ()
  {
    Grid & riverGrid = m_spRiverNetwork->grid ();
    riverGrid.detach_from_edges (m_aFSDepth);
    riverGrid.detach_from_edges (m_aRiverDepth);
    riverGrid.detach_from_edges (m_aRiverWidth);
  }
  
  void set_river_width (number riverWidth);
  
  void set_river_width (number riverWidth, int ssi);
  
  void set_river_width (number riverWidth, const char * ss_name);
  
  void use_nodal_fs_height (bool v) {m_bNodalFSDepth = v;}
  
  void check_fs_covering (bool v) {m_bCheckFSCovering = v;}
  
  void check_top_covering (bool v) {m_bCheckTopCovering = v;}
  
  void print_total_recharge (bool v) {m_bPrintTotalRecharge = v;}
  
  void print_total_recharge (const char * f_name)
  {
    m_bPrintTotalRecharge = true;
    m_sTotalRechargeFile = f_name;
  }
 
  void file_output_prefix (const char * s) {m_sOutputPrefix = s;}
  
  void reset_river_depth ();
  
  void save_river_depth (const char * file_name);
  
  SmartPtr<SegmentNetwork> river_network () const {return m_spRiverNetwork;}
  
  number fs_depth_for (const Edge * seg) const {return m_aaFSDepth [seg];}
  
  number river_depth_for (const Edge * seg) const {return m_aaRiverDepth [seg];}
 
  number river_width_for (const Edge * seg) const {return m_aaRiverWidth [seg];}
 
private:
 
  void compute_fs_depths_at_nodes ();
 
  void compute_fs_depths_at_segments ();
  
  virtual void compute_depths ();
  
  virtual number smoothed_recharge_at
  (
    MathVector<dim> & x 
  ) const;
  
  virtual number simple_smoothed_recharge_at
  (
    MathVector<dim> & x 
  ) const;
  
  void do_print_total_recharge ();
 
private:
  
  SmartPtr<SegmentNetwork> m_spRiverNetwork;
  
  depth_attachment_type m_aFSDepth;
  
  depth_accessor_type m_aaFSDepth;
  
  depth_attachment_type m_aRiverDepth;
  
  depth_accessor_type m_aaRiverDepth;
  
  depth_attachment_type m_aRiverWidth;
  
  depth_accessor_type m_aaRiverWidth;
  
  bool m_bNodalFSDepth;
  
  bool m_bCheckFSCovering;
  
  bool m_bCheckTopCovering;
  
  bool m_bPrintTotalRecharge;
  
  std::string m_sTotalRechargeFile;
  std::string m_sOutputPrefix;
};
 
} // namespace d3f
} // end namespace ug
 
#include "rivers_impl.h"
 
#endif // __H__UG__PLUGINS__D3F__RIVERS__
/* End of File */