Commit de44584d authored by Santiago Ospina's avatar Santiago Ospina

Merge branch...

Merge branch '98-implement-a-simulation-object-for-transport-transportsimulation' into 101-couple-the-transportsimulation-with-richardssimulation
parents f76510b5 16ecee35
......@@ -106,7 +106,11 @@ adding an empty line, make text **bold** or ``monospaced``.
</category>
<category name="transport.initial">
<!-- TODO! -->
<parameter name="value">
<definition> Initial constant value over the whole domain. </definition>
<values> float </values>
<suggestion> 0 </suggestion>
</parameter>
</category>
<category name="transport.time">
......@@ -174,10 +178,10 @@ adding an empty line, make text **bold** or ``monospaced``.
</category>
<category name="transport.parameters">
<parameter name="molecularDifussion">
<definition> Molecular difussion of the fluid phase </definition>
<parameter name="molecularDiffusion">
<definition> Molecular diffusion of the fluid phase </definition>
<values> float &gt; 0 </values>
<comment> 2E-9 </comment>
<suggestion> 2E-9 </suggestion>
</parameter>
</category>
......@@ -186,7 +190,16 @@ adding an empty line, make text **bold** or ``monospaced``.
<definition> Numerical scheme to perform time steps in the simulation
</definition>
<values> string </values>
<suggestion> explicit_euler, implicit_euler, heun, shu3, rk4, alex2, fractional, alex3 </suggestion>
<suggestion> explicit_euler, implicit_euler, heun, shu3, rk4, alex2, fractional, alex3
</suggestion>
</parameter>
<parameter name="courant">
<definition> Courant number for explicit methods. Explicit methods have a
maximum time step according to the CFL condition.
</definition>
<values> float &gt; 0 </values>
<suggestion> 0.8 </suggestion>
</parameter>
</category>
</dorie>
......@@ -29,9 +29,9 @@ CoupledSimulation<Traits>::CoupledSimulation(
GFWaterFluxPtr waterFlux = _richards->get_waterflux();
auto time_begin = _richards->begin_time();
// minimal interval to allow printing of initial condition
TimeInterval time_domain = TimeInterval(time_begin,time_begin+1e-10);
TimeInterval time_domain = TimeInterval{time_begin,time_begin+1e-10};
_transport->set_waterflux(time_domain,waterFlux);
_transport->set_water_flux(time_domain,waterFlux);
_transport->set_saturation(time_domain,saturation);
}
......@@ -58,7 +58,7 @@ void CoupledSimulation<Traits>::step()
GFSaturationPtr saturation = _richards->get_saturation();
GFWaterFluxPtr waterFlux = _richards->get_waterflux();
TimeInterval time_domain = TimeInterval(time_begin,time_end);
TimeInterval time_domain = TimeInterval{time_begin,time_end};
auto _current_time_transport = _transport->current_time();
auto _current_time_richards = _richards->current_time();
......@@ -69,7 +69,7 @@ void CoupledSimulation<Traits>::step()
return;
}
_transport->set_waterflux(time_domain,waterFlux);
_transport->set_water_flux(time_domain,waterFlux);
_transport->set_saturation(time_domain,saturation);
while (Dune::FloatCmp::lt(_transport->current_time(),_richards->current_time()))
......
......@@ -284,9 +284,9 @@ public:
/*------------------------------------------------------------------------*//**
* @brief Construct Richards simulation.
*
* @param _grid_mapper Mapper for grid and volume/boundary data
* @param _helper Dune MPI instance controller
* @param _inifile Dune parameter file parser
* @param _inifile Dune parameter file parser
* @param _grid_mapper Mapper for grid and volume/boundary data
* @param _helper Dune MPI instance controller
*/
RichardsSimulation (
Dune::ParameterTree& _inifile,
......@@ -343,7 +343,7 @@ public:
* @brief Suggest a time step to the model.
*
* @param[in] dt Suggestion for the internal time step of the model. The
* new internal time step shall not be bigger than dt
* new internal time step shall not be bigger than dt.
*/
void suggest_timestep(double dt) final {controller->suggest_timestep(dt);}
......
......@@ -18,7 +18,6 @@ TransportSimulation<Traits>::TransportSimulation(
, output_type(_inifile.get<bool>("output.asciiVtk") ? Dune::VTK::OutputType::ascii : Dune::VTK::OutputType::base64)
{
Dune::Timer timer;
// const int verbose_low = helper.rank() == 0 ? verbose : 0;
// --- Grid Function Space ---
gfs = std::make_shared<GFS>(GFSHelper::create(gv));
......@@ -28,15 +27,15 @@ TransportSimulation<Traits>::TransportSimulation(
// --- Operator Helper Classes ---
sboundary = std::make_unique<SoluteBoundary>(inifile);
// ssource = std::make_unique<SoluteSource>(inifile,*param);
sinitial = std::make_unique<SoluteInitial>(gv,0.);
// --- Solution Vectors and Initial Condition ---
u = std::make_shared<U>(*gfs,.0);
const RF ini_value = inifile.get<RF>("initial.value",0.);
u = std::make_shared<U>(*gfs,ini_value);
Dune::PDELab::interpolate(*sinitial,*gfs,*u);
controller = std::make_unique<CalculationController>(inifile,*sboundary,helper);
_dt = controller->getDT();
std::string ts_method = _inifile.get<std::string>("numerics.timestepMethod");
if (ts_method== "explicit_euler")
......@@ -80,16 +79,18 @@ TransportSimulation<Traits>::TransportSimulation(
template<typename Traits>
void TransportSimulation<Traits>::operator_setup()
{
if (!_gfWaterFlux)
DUNE_THROW(Dune::Exception," ");
if (!_gfSaturation)
DUNE_THROW(Dune::Exception," ");
if (!gf_water_flux)
DUNE_THROW(Dune::InvalidStateException,
"Pointer to gf_water_flux is invalid!");
if (!gf_saturation)
DUNE_THROW(Dune::InvalidStateException,
"Pointer to gf_saturation is invalid!");
// --- Local Operators ---
double D_eff = inifile.get<double>("parameters.molecularDifussion");
slop = std::make_unique<SLOP>(*sboundary,_gfWaterFlux,_gfSaturation,D_eff /*inifile,*param,gv,*fboundary,*fsource*/);
double diff_coeff = inifile.get<double>("parameters.molecularDiffusion");
slop = std::make_unique<SLOP>(*sboundary,gf_water_flux,gf_saturation,diff_coeff);
tlop = std::make_unique<TLOP>(_gfSaturation/*inifile,*param,gv*/);
tlop = std::make_unique<TLOP>(gf_saturation);
Dune::PDELab::constraints(*this->gfs,*this->cc,false);
// --- Grid Operators ---
......@@ -101,6 +102,8 @@ void TransportSimulation<Traits>::operator_setup()
ls = std::make_unique<LS>();
pdesolver = std::make_unique<PDESOLVER>(*igo,*ls,1e-6);
// --- Time Step Operators ---
osm = std::make_unique<OSM>(*ts_param,*igo,*pdesolver);
osm->setVerbosityLevel(verbose+1);
......@@ -131,30 +134,50 @@ void TransportSimulation<Traits>::step()
operator_setup();
bool succeed = false;
bool exception = false;
// Suggest times step for explicit methods.
// This assumes that timestep is never increased on failure!
if (not ts_param->implicit()) {
double courant = inifile.get<double>("numerics.courant");
assert(courant>0.);
const auto cfl = Dune::Dorie::cfl_condition(*gf_water_flux);
controller->suggest_timestep(courant*cfl);
}
while (not succeed)
{
// suggest a new dt
this->suggest_timestep(controller->getDT());
const RF t = controller->getTime();
// Obtain time variables
const auto t = controller->getTime();
const auto dt = controller->getDT();
// Check whether time intervals of grid functions are valid
if (Dune::FloatCmp::gt(_waterFluxInterval.begin(),t+_dt))
DUNE_THROW(Dune::Exception,"Water Flux is in invalid state:\n" <<
"Water Flux time domain [" << _waterFluxInterval.begin() << "," << _waterFluxInterval.end() << "]\n" <<
"Transport time interval [" << t << "," << t+_dt << "]");
if (Dune::FloatCmp::gt(_saturationInterval.begin(),t+_dt))
DUNE_THROW(Dune::Exception,"Saturation is in invalid state:\n" <<
"Saturation time domain [" << _saturationInterval.begin() << "," << _saturationInterval.end() << "]\n" <<
"Transport time interval [" << t << "," << t+_dt << "]");
if (Dune::FloatCmp::lt(_waterFluxInterval.end(),t+_dt))
DUNE_THROW(Dune::Exception,"Water Flux is in invalid state:\n" <<
"Water Flux time domain [" << _waterFluxInterval.begin() << "," << _waterFluxInterval.end() << "]\n" <<
"Transport time interval [" << t << "," << t+_dt << "]");
if (Dune::FloatCmp::lt(_saturationInterval.end(),t+_dt))
DUNE_THROW(Dune::Exception,"Water Saturation is in invalid state:\n" <<
"Saturation time domain [" << _saturationInterval.begin() << "," << _saturationInterval.end() << "]\n" <<
"Transport time interval [" << t << "," << t+_dt << "]");
if (Dune::FloatCmp::gt(water_flux_interval.begin,t+dt))
DUNE_THROW(Dune::InvalidStateException,
"Water Flux is in invalid state:\n"
<< "Water Flux time domain [" << water_flux_interval.begin << ","
<< water_flux_interval.end << "]\n" <<
"Transport time interval [" << t << "," << t+dt << "]");
if (Dune::FloatCmp::gt(saturation_interval.begin,t+dt))
DUNE_THROW(Dune::InvalidStateException,
"Saturation is in invalid state:\n"
<< "Saturation time domain [" << saturation_interval.begin << ","
<< saturation_interval.end << "]\n" <<
"Transport time interval [" << t << "," << t+dt << "]");
if (Dune::FloatCmp::lt(water_flux_interval.end,t+dt))
DUNE_THROW(Dune::InvalidStateException,
"Water Flux is in invalid state:\n"
<< "Water Flux time domain [" << water_flux_interval.begin << ","
<< water_flux_interval.end << "]\n" <<
"Transport time interval [" << t << "," << t+dt << "]");
if (Dune::FloatCmp::lt(saturation_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" <<
"Transport time interval [" << t << "," << t+dt << "]");
const bool solver_warnings = verbose > 0 && helper.rank() == 0 ?
true : false;
......@@ -165,13 +188,11 @@ void TransportSimulation<Traits>::step()
if (not solver_warnings)
dwarn.push(false);
osm->apply(t, _dt, *u, *unext);
osm->apply(t, dt, *u, *unext);
if (not solver_warnings)
dwarn.pop();
u = unext;
if (not ts_param->implicit())
this->suggest_timestep(slop->suggestTimestep(_dt));
}
catch (Dune::ISTLError &e){
exception = true;
......@@ -184,9 +205,6 @@ void TransportSimulation<Traits>::step()
succeed = true;
}
// update dt for next step
_dt = controller->getDT();
// controller reacts to outcome of solution
controller->validate(exception);
......@@ -197,9 +215,9 @@ void TransportSimulation<Traits>::step()
template<typename Traits>
void TransportSimulation<Traits>::update_adapters () const
{
Cw = std::make_shared<GFSolute>(gfs,u);
assert(_gfSaturation);
Ctotal = std::make_shared<GFTotalSolute>(*Cw,*_gfSaturation);
c_w = std::make_shared<GFSolute>(gfs,u);
assert(gf_saturation);
c_total = std::make_shared<GFTotalSolute>(*c_w,*gf_saturation);
}
template<typename Traits>
......@@ -211,22 +229,11 @@ void TransportSimulation<Traits>::write_data () const
update_adapters();
if (inifile.get<bool>("output.vertexData"))
{
vtkwriter->template addVertexData<GFSolute>(Cw,"solute");
vtkwriter->template addVertexData<GFTotalSolute>(Ctotal,"solute (total)");
if (inifile.get<bool>("output.debugMode",false))
{
vtkwriter->template addVertexData<GFSaturation>(_gfSaturation,"saturation");
vtkwriter->template addVertexData<GFWaterFlux>(_gfWaterFlux,"flux");
}
vtkwriter->template addVertexData<GFSolute>(c_w,"solute");
vtkwriter->template addVertexData<GFTotalSolute>(c_total,"solute (total)");
} else {
vtkwriter->template addCellData<GFSolute>(Cw,"solute");
vtkwriter->template addCellData<GFTotalSolute>(Ctotal,"solute (total)");
if (inifile.get<bool>("output.debugMode",false))
{
vtkwriter->template addCellData<GFSaturation>(_gfSaturation,"saturation");
vtkwriter->template addCellData<GFWaterFlux>(_gfWaterFlux,"flux");
}
vtkwriter->template addCellData<GFSolute>(c_w,"solute");
vtkwriter->template addCellData<GFTotalSolute>(c_total,"solute (total)");
}
try{
......
......@@ -24,6 +24,7 @@
#include "../solver/util_controller.hh"
#include "../solver/util_writer.hh"
#include "../solver/util_grid_creator.hh"
#include "../solver/util_cfl_condition.hh"
namespace Dune{
namespace Dorie{
......@@ -33,10 +34,10 @@ 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 GFSaturationType Grid Function type of the water content
* @tparam GFWaterFluxType Grid Function type of the water flux
* @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
*/
template<class BaseTraits, class GFWaterFluxType, class GFSaturationType>
struct TransportSimulationTraits : public BaseTraits
......@@ -70,8 +71,6 @@ struct TransportSimulationTraits : public BaseTraits
using GridMapper = Dune::Dorie::GridMapper<typename BaseTraits::Grid>;
/// Class defining boundary conditions
using SoluteBoundary = Dune::Dorie::SoluteBoundary<BaseTraits>;
/// Class defining source terms
// using SoluteSource = Dune::Dorie::SoluteSource<BaseTraits,Parameters>; // TODO
/// Class defining the initial condition
using SoluteInitial = Dune::Dorie::SoluteInitial<BaseTraits>;
/// Local spatial operator
......@@ -120,6 +119,11 @@ struct TransportSimulationTraits : public BaseTraits
* for the transport equation. It can perform time-steps and print
* the solution.
*
* @todo Implement source term.
* @todo Implement more complex initial conditions.
* @todo Implement adaptivity.
* @todo Implement DG local operator.
*
* @tparam TransportSimulationTraits Traits containing the type definitions
* which this class will use.
*/
......@@ -156,8 +160,6 @@ protected:
using GridMapper = typename Traits::GridMapper;
/// Class defining boundary conditions
using SoluteBoundary = typename Traits::SoluteBoundary;
/// Class defining source terms
// using SoluteSource = typename Traits::SoluteSource; // TODO
/// Class defining the initial condition
using SoluteInitial = typename Traits::SoluteInitial;
/// Local spatial operator
......@@ -206,7 +208,6 @@ protected:
std::unique_ptr<CC> cc;
std::unique_ptr<SoluteBoundary> sboundary;
// std::unique_ptr<SoluteSource> ssource;
std::unique_ptr<SoluteInitial> sinitial;
std::unique_ptr<CalculationController> controller;
......@@ -224,31 +225,30 @@ protected:
std::shared_ptr<U> u;
mutable std::shared_ptr<GFSolute> Cw;
mutable std::shared_ptr<GFTotalSolute> Ctotal;
std::shared_ptr<GFWaterFlux> _gfWaterFlux;
std::shared_ptr<GFSaturation> _gfSaturation;
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;
TimeInterval _waterFluxInterval = TimeInterval(std::numeric_limits<TimeField>::min(),
std::numeric_limits<TimeField>::min());
TimeInterval _saturationInterval = TimeInterval(std::numeric_limits<TimeField>::min(),
std::numeric_limits<TimeField>::min());
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()};
std::unique_ptr<Writer> vtkwriter;
const int verbose;
const Dune::VTK::OutputType output_type;
TimeField _dt;
bool operator_setup_flag = true;
public:
/*------------------------------------------------------------------------*//**
/*-----------------------------------------------------------------------*//**
* @brief Construct basic simulation. Prepare setup for basic simulation.
*
* @param _grid_mapper Mapper for grid and volume/boundary data
* @param _helper Dune MPI instance controller
* @param _inifile Dune parameter file parser
* @param _inifile Dune parameter file parser
* @param _grid_mapper Mapper for grid and volume/boundary data
* @param _helper Dune MPI instance controller
*/
TransportSimulation (
Dune::ParameterTree& _inifile,
......@@ -256,54 +256,93 @@ public:
Dune::MPIHelper& _helper
);
// TODO: Define data exchange with richards
void set_waterflux(TimeInterval time_interval, std::shared_ptr<GFWaterFlux> gfWaterFlux)
{
if (!gfWaterFlux)
DUNE_THROW(Dune::InvalidStateException,"Pointer to gfWaterFlux is invalid!");
_waterFluxInterval = time_interval;
_gfWaterFlux = gfWaterFlux;
operator_setup_flag = true;
}
/*------------------------------------------------------------------------*//**
* @brief Compute a time step. Models the time step requirement of
* SimulationBase
*
* @throws Dune::Exception Fatal error occurs during computation
*/
void step () override;
// TODO: Define data exchange with richards
void set_saturation(TimeInterval time_interval, std::shared_ptr<GFSaturation> gfSaturation)
{
if (!gfSaturation)
DUNE_THROW(Dune::InvalidStateException,"Pointer to gfWaterFlux is invalid!");
_saturationInterval = time_interval;
_gfSaturation = gfSaturation;
operator_setup_flag = true;
}
/*-----------------------------------------------------------------------*//**
* @brief Method that provides the begin time of the model.
*
* @return Begin time of the model.
*/
double begin_time() const final {return controller->getBeginTime();}
void step () override;
/*-----------------------------------------------------------------------*//**
* @brief Method that provides the end time of the model.
*
* @return End time of the model.
*/
double end_time() const final {return controller->getEndTime();}
TimeField begin_time() const override
{
return inifile.get<TimeField>("time.start");
}
TimeField end_time() const override
{
return inifile.get<TimeField>("time.end");
}
/*-----------------------------------------------------------------------*//**
* @brief Method that provides the current time of the model.
*
* @return Current time of the model.
*/
double current_time() const final {return controller->getTime();}
TimeField current_time() const override
/*-----------------------------------------------------------------------*//**
* @brief Suggest a time step to the model.
*
* @param[in] dt Suggestion for the internal time step of the model. The
* new internal time step shall not be bigger than dt.
*/
void suggest_timestep(TimeField dt) final {controller->suggest_timestep(dt);}
/*-----------------------------------------------------------------------*//**
* @brief Sets a function for the water flux for a given time interval.
*
* @param[in] _time_interval The time interval.
* @param[in] _gf_water_flux The grid function water flux.
*/
void set_water_flux( TimeInterval _time_interval,
std::shared_ptr<GFWaterFlux> _gf_water_flux)
{
return controller->getTime();
if (_gf_water_flux)
{
water_flux_interval = _time_interval;
gf_water_flux = _gf_water_flux;
operator_setup_flag = true;
} else
DUNE_THROW(Dune::InvalidStateException,
"Pointer to gf_water_flux is invalid!");
}
void suggest_timestep(TimeField dt) override
/*-----------------------------------------------------------------------*//**
* @brief Sets a function for the saturation for a given time interval.
*
* @param[in] _time_interval The time interval.
* @param[in] _gf_saturation The grid function saturation
*/
void set_saturation(TimeInterval _time_interval,
std::shared_ptr<GFSaturation> _gf_saturation)
{
_dt = std::min(_dt,dt);
if (_gf_saturation)
{
saturation_interval = _time_interval;
gf_saturation = _gf_saturation;
operator_setup_flag = true;
} else
DUNE_THROW(Dune::InvalidStateException,
"Pointer to gf_saturation is invalid!");
}
// TODO: Define data exchange with richards
/*-----------------------------------------------------------------------*//**
* @brief Get a function for the solute for the current time of the
* model.
*
* @return A pointer to a grid function of solute.
*/
std::shared_ptr<GFSolute> get_solute()
{
if (!Cw)
DUNE_THROW(Dune::InvalidStateException,"Solute in invalid state!");
return Cw;
if (!c_w)
DUNE_THROW(Dune::InvalidStateException,"Pointer to c_w in invalid state!");
return c_w;
}
protected:
......
......@@ -194,18 +194,10 @@ public:
}
};
template<class TF>
class TimeInterval
{
public:
using TimeField = TF;
TimeInterval(TimeField begin, TimeField end) : _begin(begin), _end(end) {}
const TimeField& begin() const {return _begin;}
TimeField& begin() {return _begin;}
const TimeField& end() const {return _end;}
TimeField& end() {return _end;}
private:
TimeField _begin, _end;
template<typename TF>
struct TimeInterval {
TF begin;
TF end;
};
......
......@@ -10,6 +10,17 @@
namespace Dune{
namespace Dorie{
/**
* @brief Class for volume raviart thomas local finite element.
*
* @todo Find a solution for the geometry type deprecation in dune-grid.
*
* @tparam DF { description }
* @tparam RF { description }
* @tparam k { description }
* @tparam d { description }
* @tparam gt { description }
*/
template<class DF, class RF, int k, int d, Dune::GeometryType::BasicType gt>
class VolumeRaviartThomasLocalFiniteElement
{
......@@ -42,7 +53,9 @@ public:
static constexpr GeometryType type()
{
return GeometryType(gt,d);
DUNE_NO_DEPRECATED_BEGIN
return GeometryType(gt,d);
DUNE_NO_DEPRECATED_END
}
private:
......
......@@ -9,106 +9,106 @@
#include "util_boundary_condition.hh"
namespace Dune{
namespace Dorie{
namespace Dorie{
/// Boundary type and condition value queries
/** This class containts functions that return the type of boundary conditions
* and the values of the boundary condition parameters. Both types of queries
* can be time dependent.
*/
template<typename Traits>
class SoluteBoundary
{
private:
/// Boundary type and condition value queries
/** This class containts functions that return the type of boundary conditions
* and the values of the boundary condition parameters. Both types of queries
* can be time dependent.
*/
template<typename Traits>
class SoluteBoundary
{
private:
typedef typename Traits::RangeField RF;
typedef typename Traits::Domain Domain;
typedef typename Traits::IntersectionDomain ID;
typedef typename Traits::Intersection Intersection;
typedef typename Traits::RangeField RF;
typedef typename Traits::Domain Domain;
typedef typename Traits::IntersectionDomain ID;
typedef typename Traits::Intersection Intersection;
enum {dim = Traits::dim};
enum {dim = Traits::dim};
//! Object for handling the BC data file and function queries
std::unique_ptr<BCReadoutInterface<Traits>> bcDataHandler;
//! Object for handling the BC data file and function queries
std::unique_ptr<BCReadoutInterface<Traits>> bcDataHandler;
public:
public: