The TS-GitLab will have to shut down towards the end of the year — please think about migrating your projects to GitLab.com or GitHub.
(This is still a very early message, meant to keep you informed. There will be more urgent ones in the future.)

Commit dea7b801 authored by Simon Lüdke's avatar Simon Lüdke

Read in Boundary Condition from plant

parent 83ff29d1
......@@ -16,7 +16,7 @@
#include <dune/dorie/common/boundary_condition/dirichlet.hh>
#include <dune/dorie/common/boundary_condition/neumann.hh>
#include <dune/dorie/common/boundary_condition/outflow.hh>
#include <dune/dorie/common/boundary_condition/plantroot.hh>
#include <dune/dorie/common/boundary_condition/feddes.hh>
namespace Dune {
namespace Dorie {
......@@ -316,7 +316,7 @@ public:
}
else if(type == FeddesBoundaryCondition<RF>::type_str){
// all transpiration fluxes can be interpreted as a Neumann-BC ove Feddes Boundary therefore is very simmilar
return std::make_shared<FeddesBoundaryCondition<RF>>(name, type, t_int, v_begin, v_end, false);
return std::make_shared<FeddesBoundaryCondition<RF>>(name, t_int, v_begin, v_end, false);
}
// should never reach this part...
else {
......
......@@ -292,16 +292,7 @@ private:
BCList conditions_list;
for (const auto& condition : values["conditions"]) {
const auto name = condition.first.as<std::string>();
_log->trace(" Creating boundary condition. "
"Mapping index: {}, Name: {}",
index, name);
conditions_list.push_back(_bcf.create(name,
condition.second,
time_end));
}
if (conditions_list.empty()) {
_log->error("Failed to read boundary conditions. "
_log->trace(" Creating bouboundary_index_map boundary conditions. "
"Boundary: {}, Index: {}",
key.as<std::string>(), index);
DUNE_THROW(IOError, "Error reading conditions associated "
......@@ -493,6 +484,13 @@ private:
cfg_default,
t_end);
}
public:
/// Retrive the BoundaryCondition as a shared pointer given it's index.
std::shared_ptr<BC> get_by_index(int const & index){
return _bc_cache[index];
}
};
} // namespace Dorie
......
#ifndef DUNE_DORIE_COMMON_BOUNDARY_CONDITION_PLANT_HH
#define DUNE_DORIE_COMMON_BOUNDARY_CONDITION_PLANT_HH
#include <dune/dorie/common/boundary_condition/boundary_condition.hh>
namespace Dune {
namespace Dorie {
/// Class representing a Neumann boundary condition (fixed flux)
/** In the \ref RichardsModel model, this class typically models irrigation and
* evaporation from soil surfaces and therefore stores an additional switch
* which can be used to consider the surface tilt when evaluating the boundary
* condition.
*
* If #_horizontal_projection is set to true, the local operator scales the
* boundary condition \f$ j \f$ with
* \f[
* j = j_N \left| \mathbf{\hat{n}}_F \cdot \mathbf{\hat{g}}_F \right|
* \f]
* where the vectors are the unit normal vector of the boundary face, and
* the unit vector along the gravitational force, respectively.
*
* #_horizontal_projection is ignored in the \ref TransportModel model.
*
* \ingroup Boundary
*/
template<typename RF>
class FeddesBoundaryCondition : public BoundaryCondition<RF>
{
protected:
using Base = BoundaryCondition<RF>;
using Time = typename Base::Time;
/// Value of the surface adjustment switch
const bool _horizontal_projection;
public:
/// The type string of this boundary condition class
static constexpr std::string_view type_str = "feddes";
/// Construct a static boundary condition
FeddesBoundaryCondition (const std::string_view name,
const std::string_view type,
const TimeInterval<Time> interval,
const RF value,
const bool horizontal_projection)
:
Base(name, interval, value),
_horizontal_projection(horizontal_projection),
type_str(type)
{ }
/// Construct an interpolated boundary condition
FeddesBoundaryCondition (const std::string_view name,
const std::string_view type,
const TimeInterval<Time> interval,
const RF value_begin,
const std::optional<RF> value_end,
const bool horizontal_projection)
:
Base(name, interval, value_begin, value_end),
_horizontal_projection(horizontal_projection),
type_str(type)
{ }
/// Destroy this object
virtual ~FeddesBoundaryCondition () override = default;
/// Return the type of this boundary condition
virtual Operator::BCType type () const override
{
return Operator::BCType::Feddes;
}
};
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_COMMON_BOUNDARY_CONDITION_PLANT_HH
......@@ -28,6 +28,7 @@ namespace Dorie{
typedef typename Traits::TimeField TF;
typedef typename Traits::Domain Domain;
typedef typename Traits::Element Element;
typedef typename Traits::GridView GridView;
using BCM = Dune::Dorie::BoundaryConditionManager<Traits, BCMMode::feddes_source>;
using PlantType = Plant<Traits,BCM>;
enum {dim=Traits::dim};
......@@ -38,19 +39,23 @@ namespace Dorie{
const bool _use_roots;
std::vector<PlantType> plant_vector;
std::shared_ptr<BCM> plant_bc;
public:
FlowSource(const Dune::ParameterTree& config)
FlowSource(const Dune::ParameterTree& config, std::shared_ptr<GridView> grid_view)
: _config(config)
, _plant_param_filepath(config.get<std::string>("FlowSource.PlantParameter"))
, _use_roots(config.get<bool>("FlowSource.use_roots"))
{
if(_use_roots){
initialize();
initialize(grid_view);
}
}
void initialize(){
/// Reads the data from the _plant_param_filepath and splits the yaml node
/// into a part for the bc manager
/// and a part for the plant initialization
void initialize(std::shared_ptr<GridView> grid_view){
// read in the yaml file
YAML::Node plant_param_yaml;
try{
......@@ -60,6 +65,7 @@ namespace Dorie{
_log->error("unsuccesfully read yaml file {}", _plant_param_filepath);
return;
}
TF end_time = _config.get<TF>("richards.time.end");
// create boundary condition index map
int number_of_plants = plant_param_yaml["NumberOfPlants"].as<int>();
......@@ -70,10 +76,14 @@ namespace Dorie{
// initialize the plants
YAML::Node plant_node = plant_param_yaml["Plants"];
initialize_plants (plant_node, number_of_plants);
initialize_plants (plant_node, number_of_plants, grid_view, end_time);
}
/// initialize the boundary condition manager for plants
/** \param bc_node Yaml Node containing the information about the plant boundarys. Is passed to the bc manager.
* \param number_of_plants
*/
void initialize_bc_manager(YAML::Node const & bc_node, int const & number_of_plants){
_log->info("Initialize plant-BCM ...");
std::vector<int> plant_bc_index_map;
......@@ -81,19 +91,25 @@ namespace Dorie{
for (int index=1; index<number_of_plants+1; index++){
plant_bc_index_map.emplace_back(index);
}
plant_bc = std::make_shared<BCM>(_config,plant_bc_index_map,_log,bc_node);
plant_bc = std::make_shared<BCM>(_config, plant_bc_index_map, _log, bc_node);
_log->info("Initialized the plant-boundary manager");
}
void initialize_plants(YAML::Node const & plant_node, int const & number_of_plants){
/// takes care of plant initialization
/** \param plant_node Yaml node containing the information of the plants
* \param number_of_plants Number of Plants that should be initialized
* \param grid_view Shared Pointer to a grid view
* \param end_time Time of simulation end.
*/
void initialize_plants(YAML::Node const & plant_node, int number_of_plants, std::shared_ptr<GridView> grid_view, TF end_time){
_log->info("Initialize the plant vector ...");
plant_vector.reserve(number_of_plants);
int index=1;
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, values));
plant_vector.emplace_back(PlantType(index, plant_bc, values, grid_view, end_time));
index++;
_log->info("Initialized the vector containing the plants");
}
......
......@@ -3,8 +3,10 @@
#include <string>
#include <memory>
#include <array>
#include <yaml-cpp/yaml.h>
#include <dune/dorie/common/boundary_condition/boundary_condition.hh>
#include <dune/dorie/common/logging.hh>
namespace Dune{
namespace Dorie{
......@@ -14,12 +16,18 @@ namespace Dorie{
///
template<typename Traits, typename BCM>
class Plant{
typedef typename Traits::TimeField TF; //unit for time normaly double
typedef typename Traits::RangeField RF;
typedef typename Traits::Domain Domain;
typedef typename Traits::Vector Vector;
typedef typename Traits::GridView GridView;
/// Alias of the (base) type of boundary conditions
using BC = BoundaryCondition<RF>;
std::shared_ptr<BCM> boundary_condition_vec;
/// Alias of the Intersection Type
using Intersection = typename Traits::Intersection;
std::shared_ptr<spdlog::logger> _log{get_logger("base")};
std::shared_ptr<BCM> boundary_condition_list;
std::string bc_type;
RF bc_value;
......@@ -28,27 +36,66 @@ class Plant{
const RF h1, h2, h3 ,h4;
const RF depth1, depth2;
const bool growth_enabled;
const TF end_time;
bool calculate_new_root_distribution;
Vector location {0.01,0};
RF surface_area;
public:
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>())
{}
Plant(int index, std::shared_ptr<BCM> bc_plant, YAML::Node const & param_yaml, std::shared_ptr<GridView> grid_view, TF end_time)
: boundary_condition_list(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>())
, end_time(end_time)
{
initialize(grid_view);
}
/// finds Intersection with the top BoundaryCondition at Seed-location and saves it as a pointer
/// this intersection is used to access the boundary condition and the surface area of one cell
void initialize(std::shared_ptr<GridView> grid_view) {
//change that if 3 dim is needed
for (const auto & element : elements(*grid_view)){
if (not element.hasBoundaryIntersections())
continue;
const auto geo_type = element.type();
const auto &ref = Dune::ReferenceElements<typename Traits::DF, Traits::dim>::general(geo_type);
auto pos_local = element.geometry().local(location);
if (not ref.checkInside(pos_local))
continue;
for (const auto &intersec : intersections(*grid_view, element)){
if (intersec.boundary()){
// if the intersection is an boundary it is saved in the pointer
// check if the found boundary intersection has Feddes Boundary condition
if (intersec.boundarySegmentIndex()==1){ // 1 is the index for the upper Boundary
surface_area = intersec.geometry().volume();
return;
}
}
}
}
_log->error("Seed loc has no intersection with a boundary");
}
/// 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);
if (time < end_time)
{
boundary_condition_list->set_time(time);
}
// Index-1 because plant index 1.. while bc index 0..
bc_value = boundary_condition_list->get_by_index(index-1)->evaluate(time);
}
/// Sets the state needed for a new timestep.
......@@ -107,9 +154,10 @@ public:
}
/// calculates the potential Transpiration
/// if unstructured grids are used this area calculation has to be changed
/// if different types of BC are used functionality has to be added
RF potential_transpiration() const {
// has to include bc_type and possibly surface area
return bc_value;
return bc_value*surface_area;
}
/* // test functions
......@@ -124,4 +172,7 @@ public:
} //namespace Dorie
} // namespace Dune
#endif
\ No newline at end of file
#endif
......@@ -57,7 +57,7 @@ ModelRichards<Traits>::ModelRichards (
// --- Operator Helper Classes ---
const auto& boundary_map = _grid_creator.boundary_index_map();
fboundary = std::make_shared<FlowBoundary>(inifile, boundary_map, this->_log);
fsource = std::make_shared<FlowSource>(inifile);
fsource = std::make_shared<FlowSource>(inifile, std::make_shared<GV>(gv) );
finitial = FlowInitialFactory::create(inifile, gv, fparam, this->_log);
// Read local operator settings
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment