Commit f8159599 authored by Lukas Riedel's avatar Lukas Riedel

Add solute concentration type toggle to Dirichlet BC

The parameter `concentration_type` for a Dirichlet BC of the Transport
model sets if the given concentration is a concentration in the water
phase or the total volume concentration.

* Add enum switch SoluteConcentration to *any* Dirichlet BC.
* Read parameter from boundary condition file in Transport mode.
* Adapt Transport FV local operator to read the setting of a Dirichlet
  BC.
* Add check for solute concentration toggle to BC unit test.
parent ff047eec
......@@ -6,6 +6,15 @@
namespace Dune {
namespace Dorie {
/// Toggle for concentration represented by Dirichlet BC
/** Used for transport simulation only.
*/
enum class SoluteConcentration {
total, //!< Total concentration in the volume \f$ C_t \f$.
water_phase, //!< Concentration in water phase \f$ C_w \f$.
none //!< No type, used for Richards simulation.
};
/// Class representing a Dirichlet boundary condition
/**
* \ingroup Boundary
......@@ -17,9 +26,32 @@ protected:
using Base = BoundaryCondition<RF>;
using Time = typename Base::Time;
/// Setting for the type of concentration encoded in this BC value
const SoluteConcentration _solute_concentration;
public:
/// Re-use base constructor
using Base::Base;
/// Construct a static boundary condition
DirichletBoundaryCondition (
const std::string_view name,
const TimeInterval<Time> interval,
const RF value,
const SoluteConcentration concentration=SoluteConcentration::none)
:
Base(name, interval, value),
_solute_concentration(concentration)
{ }
/// Construct an interpolated boundary condition
DirichletBoundaryCondition (
const std::string_view name,
const TimeInterval<Time> interval,
const RF value_begin,
const std::optional<RF> value_end,
const SoluteConcentration concentration=SoluteConcentration::none)
:
Base(name, interval, value_begin, value_end),
_solute_concentration(concentration)
{ }
/// Destroy this object
virtual ~DirichletBoundaryCondition () override = default;
......@@ -32,6 +64,12 @@ public:
/// The type string of this boundary condition class
static constexpr std::string_view type_str = "dirichlet";
/// Return the quantity encoded in this Dirichlet value
SoluteConcentration concentration_type () const
{
return _solute_concentration;
}
};
} // namespace Dorie
......
......@@ -165,6 +165,32 @@ private:
return std::make_pair(v_begin, v_end);
}
/// Return the concentration type for solute transport Dirichlet BCs
static SoluteConcentration get_solute_c_type (
const YAML::Node& config,
const std::shared_ptr<spdlog::logger> log
)
{
if (not config) {
log->error("Node 'concentration_type' must be given for "
"solute transport boundary conditions");
DUNE_THROW(IOError, "BC 'concentration_type' not specified");
}
const auto type_str = config.as<std::string>();
if (type_str == "water_phase") {
return SoluteConcentration::water_phase;
}
else if (type_str == "total") {
return SoluteConcentration::total;
}
else {
log->error("Unknown solute concentration type: {}",
type_str);
DUNE_THROW(IOError, "Unknown 'concentration_type' in BC");
}
}
public:
/// Create a boundary condition from a YAML node
......@@ -211,9 +237,21 @@ public:
// create the object
if (type == DirichletBoundaryCondition<RF>::type_str) {
return std::make_shared<DirichletBoundaryCondition<RF>>(
name, t_int, v_begin, v_end
);
// check solute concentration type for transport sim
if constexpr (mode == BCMMode::transport) {
const auto solute_c_type = get_solute_c_type(
config["concentration_type"],
log
);
return std::make_shared<DirichletBoundaryCondition<RF>>(
name, t_int, v_begin, v_end, solute_c_type
);
}
else {
return std::make_shared<DirichletBoundaryCondition<RF>>(
name, t_int, v_begin, v_end
);
}
}
else if (type == OutflowBoundaryCondition<RF>::type_str) {
return std::make_shared<OutflowBoundaryCondition<RF>>(
......@@ -223,8 +261,11 @@ public:
else if (type == NeumannBoundaryCondition<RF>::type_str) {
// optional parameter: "horizontal_projection"
bool project = false;
if (const auto node_project = config["horizontal_projection"]) {
project = node_project.as<bool>();
if constexpr (mode != BCMMode::transport) {
const auto node_project = config["horizontal_projection"];
if (node_project) {
project = node_project.as<bool>();
}
}
return std::make_shared<NeumannBoundaryCondition<RF>>(
......
......@@ -14,6 +14,8 @@
#include <dune/dorie/common/typedefs.hh>
#include <dune/dorie/common/boundary_condition/dirichlet.hh>
namespace Dune {
namespace Dorie {
namespace Operator {
......@@ -118,11 +120,9 @@ public:
*/
TransportFVSpatialOperator(
std::shared_ptr<const Parameter> param,
std::shared_ptr<const Boundary> boundary,
DirichletMode::Type dirichlet_mode = DirichletMode::SoluteConcentration
std::shared_ptr<const Boundary> boundary
) : _param(param)
, _boundary(boundary)
, _dirichlet_mode(dirichlet_mode)
{
assert(_gf_water_flux);
assert(_gf_water_content);
......@@ -372,14 +372,16 @@ public:
auto g = bc_value;
// Convert input to total solute if needed
if (_dirichlet_mode == DirichletMode::TotalSolute)
const auto& bc_dirichlet
= dynamic_cast<const DirichletBoundaryCondition<RF>&>(*bc);
if (bc_dirichlet.concentration_type() == SoluteConcentration::total)
{
if (Dune::FloatCmp::gt(water_content_i,WaterContent{0.}))
g /=water_content_i;
else
g = 0.;
}
// Compute the effective hydrodynamic dispersion
auto hydr_disp_i = hydr_disp_func_i(water_flux,water_content_i);
......@@ -451,7 +453,6 @@ private:
const std::shared_ptr<const Boundary> _boundary;
std::shared_ptr<GFWaterFlux> _gf_water_flux;
std::shared_ptr<GFWaterContent> _gf_water_content;
const DirichletMode::Type _dirichlet_mode;
double _time;
};
......
......@@ -64,7 +64,8 @@ void test_neumann_bc (const std::shared_ptr<BC> bc,
/// Test the Dirichlet boundary condition, after specs of the YAML input file
void test_dirichlet_bc (const std::shared_ptr<BC> bc,
const typename Traits::RF time)
const typename Traits::RF time,
const bool transport_mode)
{
using namespace Dune::FloatCmp;
assert(ge(time, 0.5) and le(time, 1.5));
......@@ -72,6 +73,20 @@ void test_dirichlet_bc (const std::shared_ptr<BC> bc,
assert(eq(bc->evaluate(1.0), 0.5));
assert(eq(bc->evaluate(1.5), 1.0));
// perform a cast
using namespace Dune::Dorie;
const auto& bc_dirichlet = dynamic_cast<
const DirichletBoundaryCondition<double>&>(*bc);
// case: not transport
if (transport_mode) {
assert(bc_dirichlet.concentration_type()
== SoluteConcentration::water_phase);
}
else {
assert(bc_dirichlet.concentration_type()
== SoluteConcentration::none);
}
// check interval testing
try {
const auto log = Dune::Dorie::get_logger(Dune::Dorie::log_base);
......@@ -98,7 +113,8 @@ void test_outflow_bc (const std::shared_ptr<BC> bc,
* specialized to the type and respective values of the BC
*/
void test_boundary_condition (const std::shared_ptr<BC> bc,
const typename Traits::RF time)
const typename Traits::RF time,
const bool transport_mode)
{
using namespace Dune::Dorie;
const auto log = get_logger(log_base);
......@@ -111,7 +127,7 @@ void test_boundary_condition (const std::shared_ptr<BC> bc,
test_default_bc(bc, time);
}
else if (type == Operator::BCType::Dirichlet) {
test_dirichlet_bc(bc, time);
test_dirichlet_bc(bc, time, transport_mode);
}
else if (type == Operator::BCType::Neumann) {
test_neumann_bc(bc, time);
......@@ -129,7 +145,8 @@ void test_boundary_condition (const std::shared_ptr<BC> bc,
template<typename BCManager>
void test_bcs_on_grid (const std::shared_ptr<typename Traits::Grid> grid,
const BCManager& boundary,
const typename Traits::RF time)
const typename Traits::RF time,
const bool transport_mode)
{
const auto gv = grid->leafGridView();
......@@ -145,7 +162,7 @@ void test_bcs_on_grid (const std::shared_ptr<typename Traits::Grid> grid,
continue;
const auto bc = boundary.bc(is);
test_boundary_condition(bc, time);
test_boundary_condition(bc, time, transport_mode);
}
}
}
......@@ -198,7 +215,7 @@ int main(int argc, char** argv)
log->info("Setting global time to: {}", time);
boundary.set_time(time);
max_dt = boundary.max_time_step(time);
test_bcs_on_grid(grid, boundary, time);
test_bcs_on_grid(grid, boundary, time, false);
}
// check behavior when reaching end of final BC intervals
......@@ -207,6 +224,15 @@ int main(int argc, char** argv)
assert(test_set_time_throw(boundary));
log->info("...everything went as expected!");
log->info("Initializing the boundary manager in Transport mode");
Dune::Dorie::BoundaryConditionManager<
Traits, Dune::Dorie::BCMMode::transport>
bc_transport (inifile, index_map);
double time = 1.5;
bc_transport.set_time(time);
test_bcs_on_grid(grid, bc_transport, time, true);
return 0;
}
catch (Dune::Exception &e){
......
......@@ -14,6 +14,7 @@ lower:
type: Dirichlet
value: [0.0, 1.0]
time: [0.5, 1.5]
concentration_type: water_phase
left:
index: 2
......
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