Use water content instead of saturation for transport solver

Signed-off-by: Santiago Ospina De Los Ríos's avatarSantiago Ospina <santiago.ospina@iup.uni-heidelberg.de>
parent b0386663
......@@ -22,7 +22,7 @@
`build-cmake/test/coverage`. The total coverage is reported to GitLab.
* New class [`SimulationTransport`](dune/dorie/interface/transport_simulation.hh)
which solves the transport equation with a finite volume scheme for give
grid functions of saturation and velocity.
grid functions of water content and velocity.
### Changed
* `Simulation` is renamed `RichardsSimulation` and moved to
......
......@@ -82,15 +82,15 @@ void TransportSimulation<Traits>::operator_setup ()
if (!gf_water_flux)
DUNE_THROW(Dune::InvalidStateException,
"Pointer to gf_water_flux is invalid!");
if (!gf_saturation)
if (!gf_water_content)
DUNE_THROW(Dune::InvalidStateException,
"Pointer to gf_saturation is invalid!");
"Pointer to gf_water_content is invalid!");
// --- Local Operators ---
double diff_coeff = inifile.get<double>("parameters.molecularDiffusion");
slop = std::make_unique<SLOP>(*sboundary,gf_water_flux,gf_saturation,diff_coeff);
slop = std::make_unique<SLOP>(*sboundary,gf_water_flux,gf_water_content,diff_coeff);
tlop = std::make_unique<TLOP>(gf_saturation);
tlop = std::make_unique<TLOP>(gf_water_content);
Dune::PDELab::constraints(*this->gfs,*this->cc,false);
// --- Grid Operators ---
......@@ -158,11 +158,11 @@ void TransportSimulation<Traits>::step ()
<< water_flux_interval.end << "]\n" <<
"Transport time interval [" << t << "," << t+dt << "]");
if (Dune::FloatCmp::gt(saturation_interval.begin,t+dt))
if (Dune::FloatCmp::gt(water_content_interval.begin,t+dt))
DUNE_THROW(Dune::InvalidStateException,
"Saturation is in invalid state:\n"
<< "Saturation time domain [" << saturation_interval.begin << ","
<< saturation_interval.end << "]\n" <<
"Water Content is in invalid state:\n"
<< "Water Content time domain [" << water_content_interval.begin << ","
<< water_content_interval.end << "]\n" <<
"Transport time interval [" << t << "," << t+dt << "]");
if (Dune::FloatCmp::lt(water_flux_interval.end,t+dt))
......@@ -172,11 +172,11 @@ void TransportSimulation<Traits>::step ()
<< water_flux_interval.end << "]\n" <<
"Transport time interval [" << t << "," << t+dt << "]");
if (Dune::FloatCmp::lt(saturation_interval.end,t+dt))
if (Dune::FloatCmp::lt(water_content_interval.end,t+dt))
DUNE_THROW(Dune::InvalidStateException,
"Water Saturation is in invalid state:\n"
<< "Saturation time domain [" << saturation_interval.begin << ","
<< saturation_interval.end << "]\n" <<
"Water Water Content is in invalid state:\n"
<< "Water Content time domain [" << water_content_interval.begin << ","
<< water_content_interval.end << "]\n" <<
"Transport time interval [" << t << "," << t+dt << "]");
const bool solver_warnings = verbose > 0 && helper.rank() == 0 ?
......@@ -215,8 +215,8 @@ template<typename Traits>
void TransportSimulation<Traits>::update_adapters () const
{
c_w = std::make_shared<GFSolute>(gfs,u);
assert(gf_saturation);
c_total = std::make_shared<GFTotalSolute>(*c_w,*gf_saturation);
assert(gf_water_content);
c_total = std::make_shared<GFTotalSolute>(*c_w,*gf_water_content);
}
template<typename Traits>
......
......@@ -34,12 +34,12 @@ namespace Dorie{
* class extends BaseTraits to be used for the DG implementation of
* the Transport simulation
*
* @tparam BaseTraits Traits defining domain and range field
* properties of the simulation.
* @tparam GFWaterFluxType Grid Function type of the water flux
* @tparam GFSaturationType Grid Function type of the saturation
* @tparam BaseTraits Traits defining domain and range field
* properties of the simulation.
* @tparam GFWaterFluxType Grid Function type of the water flux
* @tparam GFWaterContentType Grid Function type of the water content
*/
template<class BaseTraits, class GFWaterFluxType, class GFSaturationType>
template<class BaseTraits, class GFWaterFluxType, class GFWaterContentType>
struct TransportSimulationTraits : public BaseTraits
{
static constexpr int dim = BaseTraits::dim;
......@@ -54,9 +54,9 @@ struct TransportSimulationTraits : public BaseTraits
"Water Flux Grid Function has wrong range");
using GFWaterFlux = GFWaterFluxType;
/// (Inestationary) Grid Function Type for water content
static_assert(GFSaturationType::Traits::dimRange == 1,
"Saturation Grid Function has wrong range");
using GFSaturation = GFSaturationType;
static_assert(GFWaterContentType::Traits::dimRange == 1,
"Water Content Grid Function has wrong range");
using GFWaterContent = GFWaterContentType;
/// GFS Helper
using GFSHelper = Dune::Dorie::LinearSolverGridFunctionSpaceHelper<GV,RF,
......@@ -76,9 +76,9 @@ struct TransportSimulationTraits : public BaseTraits
/// Local spatial operator
using SLOP = Dune::Dorie::Operator::TransportFVSpatialOperator<SoluteBoundary,
GFWaterFlux,
GFSaturation>;
GFWaterContent>;
/// Local temporal operator
using TLOP = Dune::Dorie::Operator::TransportFVTemporalOperator<GFSaturation>;
using TLOP = Dune::Dorie::Operator::TransportFVTemporalOperator<GFWaterContent>;
/// Time controller
using CalculationController = Dune::Dorie::CalculationController<RF,SoluteBoundary>;
......@@ -106,7 +106,7 @@ struct TransportSimulationTraits : public BaseTraits
/// Solute
using GFSolute = Dune::PDELab::DiscreteGridFunction<GFS,U>;
using GFTotalSolute = Dune::PDELab::ProductGridFunctionAdapter<GFSolute,GFSaturation>;
using GFTotalSolute = Dune::PDELab::ProductGridFunctionAdapter<GFSolute,GFWaterContent>;
// -- Utility Class Definitions -- //
/// VTK Output writer base class
......@@ -146,7 +146,7 @@ protected:
/// (Inestationary) Grid Function Type for water flux
using GFWaterFlux = typename Traits::GFWaterFlux;
/// (Inestationary) Grid Function Type for water content
using GFSaturation = typename Traits::GFSaturation;
using GFWaterContent = typename Traits::GFWaterContent;
/// GFS Helper
using GFSHelper = typename Traits::GFSHelper;
......@@ -228,12 +228,12 @@ protected:
mutable std::shared_ptr<GFSolute> c_w;
mutable std::shared_ptr<GFTotalSolute> c_total;
std::shared_ptr<GFWaterFlux> gf_water_flux;
std::shared_ptr<GFSaturation> gf_saturation;
std::shared_ptr<GFWaterContent> gf_water_content;
TimeInterval water_flux_interval = TimeInterval{std::numeric_limits<TimeField>::min(),
std::numeric_limits<TimeField>::min()};
TimeInterval saturation_interval = TimeInterval{std::numeric_limits<TimeField>::min(),
std::numeric_limits<TimeField>::min()};
TimeInterval water_flux_interval = TimeInterval{std::numeric_limits<TimeField>::lowest(),
std::numeric_limits<TimeField>::lowest()};
TimeInterval water_content_interval = TimeInterval{std::numeric_limits<TimeField>::lowest(),
std::numeric_limits<TimeField>::lowest()};
std::unique_ptr<Writer> vtkwriter;
......@@ -313,21 +313,22 @@ public:
}
/*-----------------------------------------------------------------------*//**
* @brief Sets a function for the saturation for a given time interval.
* @brief Sets a function for the water content for a given time
* interval.
*
* @param[in] _time_interval The time interval.
* @param[in] _gf_saturation The grid function saturation
* @param[in] _time_interval The time interval.
* @param[in] _gf_water_content The grid function water content
*/
void set_saturation(TimeInterval _time_interval,
std::shared_ptr<GFSaturation> _gf_saturation)
void set_water_content(TimeInterval _time_interval,
std::shared_ptr<GFWaterContent> _gf_water_content)
{
if (_gf_saturation)
if (_gf_water_content)
{
saturation_interval = _time_interval;
gf_saturation = _gf_saturation;
water_content_interval = _time_interval;
gf_water_content = _gf_water_content;
} else
DUNE_THROW(Dune::InvalidStateException,
"Pointer to gf_saturation is invalid!");
"Pointer to gf_water_content is invalid!");
}
/*-----------------------------------------------------------------------*//**
......
......@@ -110,20 +110,20 @@ int main(int argc, char** argv) {
typename Traits::Scalar sat;
flux[dim-1] = inifile.get<double>("_constantWater.flux",-5.5e-8);
sat = inifile.get<double>("_constantWater.saturation",1.);
sat = inifile.get<double>("_constantWater.waterContent",1.);
using GFFlux = typename Traits::GFWaterFlux;
using GFSat = typename Traits::GFSaturation;
using GFSat = typename Traits::GFWaterContent;
auto grid = grid_mapper->get_grid();
auto gf_flux = std::make_shared<GFFlux>(grid->leafGridView(),flux);
auto gf_sat = std::make_shared<GFSat>(grid->leafGridView(),sat);
auto gf_theta = std::make_shared<GFSat>(grid->leafGridView(),sat);
typename Traits::TimeInterval interval{sim.begin_time(),sim.end_time()};
sim.set_water_flux(interval,gf_flux);
sim.set_saturation(interval,gf_sat);
sim.set_water_content(interval,gf_theta);
sim.run();
} else {
......@@ -141,23 +141,23 @@ int main(int argc, char** argv) {
Sim<Traits> sim(inifile_transport,grid_mapper,helper);
typename Traits::Vector flux;
typename Traits::Scalar sat;
typename Traits::Scalar theta;
flux[dim-1] = inifile.get<double>("_constantWater.flux",-5.5e-8);
sat = inifile.get<double>("_constantWater.saturation",1.);
theta = inifile.get<double>("_constantWater.waterContent",1.);
auto grid = grid_mapper->get_grid();
using GFFlux = typename Traits::GFWaterFlux;
using GFSat = typename Traits::GFSaturation;
using GFSat = typename Traits::GFWaterContent;
auto gf_flux = std::make_shared<GFFlux>(grid->leafGridView(),flux);
auto gf_sat = std::make_shared<GFSat>(grid->leafGridView(),sat);
auto gf_theta = std::make_shared<GFSat>(grid->leafGridView(),theta);
typename Traits::TimeInterval interval{sim.begin_time(),sim.end_time()};
sim.set_water_flux(interval,gf_flux);
sim.set_saturation(interval,gf_sat);
sim.set_water_content(interval,gf_theta);
sim.run();
} else {
......
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