param_van_genuchten.hh 4.04 KB
Newer Older
Dion Haefner's avatar
Dion Haefner committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_DORIE_PARAMETERS_VAN_GENUCHTEN_HH
#define DUNE_DORIE_PARAMETERS_VAN_GENUCHTEN_HH

namespace Dune {
  namespace Dorie {

    /// Mualem-van Genuchten parameterization class
    template<typename T, class Interpolator>
    class MualemVanGenuchten: public ParameterizationBase<T,Interpolator>
    {
    private:
      typedef ParameterizationBase<T,Interpolator> PB;
      typedef typename T::RangeField RF;
      typedef typename T::Domain Domain;
      typedef Parameter<T> Par;
      Par thetaR, thetaS, alpha, n, tau, k0;
      std::vector<Par*> parameterArray;

    public:

      explicit MualemVanGenuchten(const ParameterTree& fieldConfig_, const Dune::MPIHelper& helper_,
        Interpolator& interpolator_, H5File& h5file_, const int verbose_)
      : PB(fieldConfig_, "vanGenuchten", parameterArray, helper_, interpolator_, h5file_, verbose_)
25 26
      , thetaR("theta_r", 0.0, 1.0), thetaS("theta_s", 0.0, 1.0), alpha("alpha", -20.0, 0.0)
      , n("n", 1.0, 100.0), tau("tau", -100.0, 1.0), k0("k0", 0.0, 1.0)
Dion Haefner's avatar
Dion Haefner committed
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
      {
        parameterArray.push_back(& thetaR);
        parameterArray.push_back(& thetaS);
        parameterArray.push_back(& alpha);
        parameterArray.push_back(& n);
        parameterArray.push_back(& tau);
        parameterArray.push_back(& k0);
        PB::read_parameters();
      }

	/******************************************
	*  PARAMETERIZATION QUERIES
	*  According to Mualem & van Genuchten
	******************************************/

		  /// Return water saturation for given matric head at current position
		  /** \f[ \Theta(h) = \left( 1 + \left( \alpha h \right)^n \right)^{-m} \f]
       * with \f$m = 1 - 1/n\f$
		   * \param h Matric head
		   * \return Soil water saturation
		   */
48
		  RF saturation (const RF& h, const Domain& pos) const override
Dion Haefner's avatar
Dion Haefner committed
49
		  {
50 51 52 53 54
				if (h >= 0.)
					return 1.;
	      const RF n_x = PB::interpolate(n,pos);
	      const RF m_x = 1. - 1./n_x;
				return std::pow(1 + std::pow(PB::interpolate(alpha,pos) * h, n_x),-m_x);
Dion Haefner's avatar
Dion Haefner committed
55 56 57 58 59 60 61
		  }

		  /// Return volumetric water content for given matric head at current position
		  /** \f[ \theta(h) = \Theta(h) (\theta_s - \theta_r) + \theta_r \f]
		   * \param h Matric head
		   * \return Soil Water Content
		   */
62
		  RF waterContent(const RF& h, const Domain& pos) const override
Dion Haefner's avatar
Dion Haefner committed
63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
		  {
  			return saturation(h,pos) * (PB::interpolate(thetaS,pos) - PB::interpolate(thetaR,pos)) + PB::interpolate(thetaR,pos);
		  }

		  /// Return saturated conductivity at current position
		  /** \f[ K = K_0 \f]
		   * \return Saturated conductivity
		   */
		  RF K(const Domain& pos) const
		  {
			  return PB::interpolate(k0,pos);
		  }

		  /// Return conductivity factor for given water saturation at current position
		  /** Saturated conductivity K is multiplied with factor condFactor < 1 to receive saturation-dependent conductivity.
       * \f[ K_\text{fac} = \Theta^\tau \left( \Theta - \left( 1 - \Theta^{1/m} \right)^m \right)^2 \f]
       * with \f$m = 1 - 1/n\f$
		   * \param sat Soil water saturation
		   * \return Conductivity Factor
		   */
83
		  RF condFactor(const RF& sat, const Domain& pos) const override
Dion Haefner's avatar
Dion Haefner committed
84
		  {
85
        const RF m_x = 1. - 1./PB::interpolate(n,pos);
Dion Haefner's avatar
Dion Haefner committed
86 87 88
			  return std::pow(sat,PB::interpolate(tau,pos))
         * std::pow(1 - std::pow(1 - std::pow(sat,1/m_x),m_x),2);
		  }
89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106


		  /// Return matric head for given water content at current position
		  /**
		   *  \param waterContent Soil water content
		   *  \param pos Position in global coordinates
		   *  \return Matric head
		   *  \f[
						h_m = \alpha^{-1} * \left[ \Theta^{-1/m} -1 \right]^{1/n}
		   		]
		   */
		  RF matricHead(const RF waterContent, const Domain& pos) const override
		  {
		  	const RF n_x = PB::interpolate(n,pos);
		  	const RF m_x = 1. - 1./n_x;
		  	const RF sat_x = (waterContent - PB::interpolate(thetaR,pos)) / (PB::interpolate(thetaS,pos) - PB::interpolate(thetaR,pos));
		  	return std::pow(PB::interpolate(alpha,pos),-1) * std::pow( std::pow(sat_x,-1./m_x)-1. , 1./n_x);
		  }
Dion Haefner's avatar
Dion Haefner committed
107 108 109 110 111
    };
  }
}

#endif // DUNE_DORIE_PARAMETERS_VAN_GENUCHTEN_HH