richards_source.hh 3.13 KB
Newer Older
Dion Haefner's avatar
Dion Haefner committed
1 2 3 4
// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_DORIE_RICHARDS_SOURCE_HH
#define DUNE_DORIE_RICHARDS_SOURCE_HH

5
#include <dune/common/parametertree.hh>
Dion Haefner's avatar
Dion Haefner committed
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98

namespace Dune{
  namespace Dorie{
    /**
    * @brief Sink and source terms for flow equation
    *
    * This class implements sinks (q<0) and sources (q>0) inside of the grid.
    */
    template<typename Traits, typename Parameter>
    class FlowSource
    {

      private:

        typedef typename Traits::RangeField         RF;
        typedef typename Traits::Domain             Domain;
        typedef typename Traits::Element            Element;

        enum {dim=Traits::dim};

        const Parameter& parameter; //!< Parameterization class
        const ParameterTree& config; //!< Parameter File Parser
        const Domain extensions; //!< Spatial extensions of grid
        const RF eps; //!< Error margin

        // Root extraction
        //const RF pumprate,maxDepth;
        //const std::vector<RF> fullDist;

      public:

        FlowSource(const Dune::ParameterTree& config_, const Parameter& parameter_)
          : parameter(parameter_), config(config_)
          , extensions(config.get<Domain>("grid.extensions"))
          , eps(1e-10)
        {}


        /// Return value of the sink/source at certain position and time
        /** \param elem Element enitity
        * \param x Position in local entity coordinates
        * \param time Time value
        * \param xglobal Global coordinates
        * \return Sink/source term value
        * \todo To do: make head available in q to calculate root water uptake dependent of the potential,
        *  return correct values
        */
        RF q (const Element& elem, const Domain& x, RF time) const
        {
          /*
          if (dim == 2)
            const RF A = extensions[0];
          else if (dim == 3)
            const RF A = extensions[0]*extensions[1];
          else
            return 0.0;
          */

          return 0.0;
        }

      private:

        // Returns the root distribution at a certain depth, for further use in function q
        /*
        RF rootDist (const Element& elem, const Domain& x, RF time) const
        {
          const Domain& xGlobal = elem.geometry().global(x);
          const RF height = extensions[dim-1];
          const RF width = extensions[0];

          const RF rootStartHeight = height - fullDist[0];
          const RF rootEndHeight = height - fullDist[1];
          const RF rootMaxHeight = height - maxDepth;

          if (dim < 2)
            return 0.0;
          if (  xGlobal[dim-1]  > rootStartHeight && xGlobal[0] > eps && xGlobal[0] < width-eps )
            return (xGlobal[dim-1] - height) / fullDist[0];
          else if ( xGlobal[dim-1]  > rootEndHeight && xGlobal[dim-1]  < rootStartHeight && xGlobal[0] > eps && xGlobal[0] < width-eps )
            return 1.0;
          else if ( xGlobal[dim-1]  > rootMaxHeight && xGlobal[dim-1]  < rootEndHeight && xGlobal[0] > eps && xGlobal[0] < width-eps )
            return 1.0 - (rootEndHeight - xGlobal[dim-1]) / (rootEndHeight - rootMaxHeight);
          else
            return 0.0;
        }
        */

    };
  }
}

#endif