added the bare bones of the root water-uptake functionality

dune/dorie/model/richards/flow_source.hh

@@ -7,7 +7,7 @@
 #include <vector>
 #include <memory>
 #include <optional>
-#include <yaml-cpp/yaml.h>
+    #include <yaml-cpp/yaml.h>
 #include <dune/common/parametertree.hh>
 #include <dune/dorie/model/richards/flow_source/feddes_root.hh>
 #include <dune/dorie/common/boundary_condition/manager.hh>
@@ -89,14 +89,17 @@ namespace Dorie{
             for (const auto&& e : plant_node){
                 //auto name = e.first.as<std::string>();
                 const auto& values = e.second;
-                plant_vector.emplace_back(PlantType(index, plant_bc));
+                plant_vector.emplace_back(PlantType(index, plant_bc, values));
                 index++;
                 _log->info("Initialized the vector containing the plants");
             }
         }
-
+        /// Sets time for the plant boundary condition and signals new timestep to the plants.
         void set_time(TF const & time){
             plant_bc->set_time(time);
+            for (auto & p : plant_vector){
+                p.ready_for_new_timestep(time);
+            }
         }
         std::optional<TF> max_time_step(TF const & time) const{
             if (_use_roots)
@@ -111,8 +114,9 @@ namespace Dorie{
         {
             if (_use_roots){
                 RF sink=0.0;
+                Domain const &xGlobal = elem.geometry().center();
                 for (auto & p : plant_vector){
-                    sink+=p.calc_root_water_uptake();
+                    sink+=p.calc_root_water_uptake(xGlobal, matric_head);
                 }
                 return sink;
             }

dune/dorie/model/richards/flow_source/feddes_root.hh

@@ -3,6 +3,7 @@
 
 #include <string>
 #include <memory>
+#include <yaml-cpp/yaml.h>
 #include <dune/dorie/common/boundary_condition/boundary_condition.hh>
 
 namespace Dune{
@@ -14,25 +15,101 @@ namespace Dorie{
 template<typename Traits, typename BCM>
 class Plant{
     typedef typename Traits::RangeField         RF;
+    typedef typename Traits::Domain             Domain;
     /// Alias of the (base) type of boundary conditions
     using BC = BoundaryCondition<RF>;
 
     std::shared_ptr<BCM> boundary_condition_vec;
+    std::string bc_type;
+    RF bc_value;
+
     const int index;
-    
+    const int root_dimension;
+    const RF h1, h2, h3 ,h4;
+    const RF depth1, depth2;
+    const bool growth_enabled;
+    bool calculate_new_root_distribution;
 public:
-    Plant(int index, std::shared_ptr<BCM> bc_plant)
+    Plant(int index, std::shared_ptr<BCM> bc_plant, YAML::Node const & param_yaml)
     : boundary_condition_vec(bc_plant)
     , index(index)
+    , root_dimension((Traits::dim < param_yaml["dimension"].as<int>())? Traits::dim : param_yaml["dimension"].as<int>())
+    , h1(param_yaml["alpha"]["h1"].as<RF>())
+    , h2(param_yaml["alpha"]["h2"].as<RF>())
+    , h3(param_yaml["alpha"]["h3"].as<RF>())
+    , h4(param_yaml["alpha"]["h4"].as<RF>())
+    , depth1(param_yaml["beta"]["const"].as<RF>())
+    , depth2(param_yaml["beta"]["lin"].as<RF>())
+    , growth_enabled(param_yaml["growth_enabled"].as<int>())
     {}
 
-    auto determine_bc_type(){
-        auto type = boundary_condition_vec->get_boundary_condition(index)->type();
-        return type;
+    /// Checks for the type of the boundary condition.
+    /// The check of the BC Type is nessesary because the calculation of the potential transpiration depends on it.
+
+    void set_bc(RF const & time) {
+        bc_type = boundary_condition_vec->get_boundary_condition(index)->get_type_str();
+        bc_value = boundary_condition_vec->get_boundary_condition(index)->evaluate(time);
+    }
+
+    /// Sets the state needed for a new timestep.
+    /// This includes determining the possibly new boundary condition 
+    /// and checking wether or not the root distribution has to be recalculated.
+    
+    void ready_for_new_timestep(RF const & time){
+        set_bc(time);
+        if(growth_enabled)
+            calculate_new_root_distribution=true;
+    }
+
+    /// calculation of the root water uptake.
+    ///
+    /// The calculation of the root water uptake depends on the water_stress_function() , the root_density() and the potential_Transpiration.
+    /// @attention Normalization has to be considered for different boundary conditions.
+    /// @param xGlobal global location
+    /// @param h_m  local matric head 
+
+    RF calc_root_water_uptake(Domain const &xGlobal, RF const &h_m) const {
+        RF sink_value = -1 * water_stress_factor(h_m) *root_density(xGlobal) * potential_transpiration(); // / root_density_norm 
+        return sink_value;
+    }
+    
+    /// Calculation of the so called water-stress-factor the water uptake ability depending on the matric head.
+    ///
+    /// Calculation of the water-stress-factor is based on the pice wise linear mode proposed by Feddes in 1978.
+    /// @param h_m matric head
+
+    RF water_stress_factor(RF const & h_m) const {
+        if ((h_m < h1) && (h_m > h2)){
+            return (h_m - h1) / (h2 - h1);
+        }
+        if ((h_m <= h2) && (h_m >= h3)){
+            return 1.0;
+        }
+        if ((h_m < h3) && (h_m > h4)){
+            return (h_m - h4) / (h3 - h4);
+        }
+        else{ // ((h_m<=h4)||(h_m>=h1))
+            return 0.0;
+        }
+    }
+
+    /// Calculates root density distribution
+    RF root_density(Domain const &xGlobal) const {
+        RF density;
+        if (xGlobal[0]<depth1)
+            density=1;
+        else if (xGlobal[0]<depth2)
+            density=1-(xGlobal[0]-depth1)/(depth2-depth1);
+        if (root_dimension>1){
+            // do some other stuff as well
+        }
+        return density;
     }
 
-    RF calc_root_water_uptake() const {
-        return 0.0;
+    /// calculates the potential Transpiration
+    RF potential_transpiration() const {
+        // has to include bc_type and possibly surface area
+        return bc_value;
     }
 
     // test functions
@@ -42,7 +119,7 @@ public:
     std::string get_bc_type() const {
         return boundary_condition_vec->get_boundary_condition(index)->get_type_str();
     }
-}; // Plant
+}; // end Plant
 
 } //namespace Dorie
 } // namespace Dune

dune/dorie/test/test-flowsource-param.yml

@@ -7,11 +7,10 @@ Plants:
         growth_enabled: 0
         alpha:
             type: feddes
-            feddes:
-                h1: 0.0
-                h2: -0.1
-                h3: -4.1
-                h4: -150
+            h1: 0.0
+            h2: -0.1
+            h3: -4.1
+            h4: -150
         beta:
             type: constlin
             constlin:
@@ -24,11 +23,10 @@ Plants:
         growth_enabled: 0
         alpha:
             type: feddes
-            feddes:
-                h1: 0.0
-                h2: -0.1
-                h3: -4.1
-                h4: -150
+            h1: 0.0
+            h2: -0.1
+            h3: -4.1
+            h4: -150
         beta:
             type: constlin
             constlin:
@@ -41,11 +39,10 @@ Plants:
         growth_enabled: 0
         alpha:
             type: feddes
-            feddes:
-                h1: 0.0
-                h2: -0.1
-                h3: -4.1
-                h4: -150
+            h1: 0.0
+            h2: -0.1
+            h3: -4.1
+            h4: -150
         beta:
             type: constlin
             constlin:

dune/dorie/test/test-flowsource-richards.cc

@@ -80,7 +80,6 @@ TEST_F(SourceTest, BoundaryTest){
 }
 
 TEST_F(SourceTest, BoundaryTypeTest){
-    EXPECT_EQ(sourceclass->get_plant(3)->determine_bc_type(),Dune::Dorie::Operator::BCType::none);
     EXPECT_EQ(sourceclass->get_plant(1)->get_bc_type(),"transpiration_per_surface_area");
     EXPECT_EQ(sourceclass->get_plant(2)->get_bc_type(),"transpiration_per_surface_area");
     EXPECT_EQ(sourceclass->get_plant(3)->get_bc_type(),"transpiration_per_surface_area");