Commit 30a82ba0 authored by Santiago Ospina's avatar Santiago Ospina

First support for RichardsSimulation with the new base class

Signed-off-by: default avatarSantiago Ospina <saospina@hugo.iwr.uni-heidelberg.de>
parent 2e17b93c
......@@ -34,11 +34,6 @@ public:
: SimulationBase(OutputPolicy::None,AdaptivityPolicy::None)
{}
/**
* @brief Writes the data.
*/
void virtual writeData();
/**
* @brief Sets the output policy.
*
......@@ -68,16 +63,46 @@ public:
*/
AdaptivityPolicy adaptivityPolicy() {return _adapt_policy;}
/**
* @brief Writes the data.
*/
void virtual writeData()
{
if (_output_policy != OutputPolicy::None)
DUNE_THROW(Dune::IOError,"Simulation can't write data!");
}
/**
* @brief Mark the grid in order to improve the current model.
*/
virtual void markGrid();
virtual void markGrid()
{
DUNE_THROW(Dune::NotImplemented,"Simulation can't mark the grid for adaptivity!");
}
/**
* @brief Operations before adaptation of the grid
*/
virtual void preAdaptGrid() {};
/**
* @brief Adapt the grid together it every dependency of the
* grid (e.g. solution vector and grid function spaces).
*/
virtual void adapt();
virtual void adaptGrid()
{
if (_adapt_policy != AdaptivityPolicy::None)
DUNE_THROW(Dune::NotImplemented,"Simulation can't be adapted!");
else
DUNE_THROW(Dune::InvalidStateException,
"Simulation uses adaptivity policy None. Method adapt() must not be called!");
}
/**
* @brief Operations after adaptation of the grid
*/
virtual void postAdaptGrid() {};
/**
* @brief Method that provides the begin time of the model.
......@@ -107,59 +132,48 @@ public:
*/
virtual void suggestTimestep(double dt) {};
/**
* @brief Operations before step of the simulation.
*/
virtual void preStep() {};
/**
* @brief Performs one steps in direction to endTime()
*/
virtual void step() = 0;
/**
* @brief Operations after step of the simulation.
*/
virtual void postStep() {};
/**
* @brief Runs the model performing steps until currentTime() equals endTime()
*/
virtual void run();
virtual void run()
{
if (outputPolicy() != OutputPolicy::None)
writeData();
while( Dune::FloatCmp::lt(currentTime(),endTime()) )
{
step();
if (adaptivityPolicy() != AdaptivityPolicy::None)
if ( Dune::FloatCmp::lt(currentTime(),endTime()) )
{
markGrid();
preAdaptGrid();
adaptGrid();
postAdaptGrid();
}
}
}
private:
OutputPolicy _output_policy;
AdaptivityPolicy _adapt_policy;
};
void SimulationBase::writeData()
{
if (_output_policy != OutputPolicy::None)
DUNE_THROW(Dune::IOError,"Simulation can't write data!");
}
void SimulationBase::adapt()
{
if (_adapt_policy != AdaptivityPolicy::None)
DUNE_THROW(Dune::NotImplemented,"Simulation can't be adapted!");
else
DUNE_THROW(Dune::InvalidStateException,
"Simulation uses adaptivity policy None. Method adapt() must not be called!");
}
void SimulationBase::markGrid()
{
DUNE_THROW(Dune::NotImplemented,"Simulation can't mark the grid for adaptivity!");
}
void SimulationBase::run()
{
if (outputPolicy() != OutputPolicy::None)
writeData();
while( Dune::FloatCmp::lt(currentTime(),endTime()) )
{
step();
if (adaptivityPolicy() != AdaptivityPolicy::None)
if ( Dune::FloatCmp::lt(currentTime(),endTime()) )
{
markGrid();
adapt();
}
}
}
} // namespace Dorie
} // namespace Dune
......
......@@ -115,64 +115,69 @@ void RichardsSimulation<Traits>::operator_setup ()
}
template<typename Traits>
bool RichardsSimulation<Traits>::compute_time_step ()
void RichardsSimulation<Traits>::step ()
{
const RF t = controller->getTime();
const RF dt = controller->getDT();
bool exception = false;
const bool parmetis_warning = verbose > 0 && helper.rank() == 0 ?
true : false;
try
bool step_succed = false;
while(controller->doStep() & not step_succed)
{
pdesolver->setMaxIterations(controller->getIterations());
const RF t = controller->getTime();
const RF dt = controller->getDT();
bool exception = false;
const bool parmetis_warning = verbose > 0 && helper.rank() == 0 ?
true : false;
if (not parmetis_warning)
dwarn.push(false);
osm->apply(t, dt, *uold, *unew);
if (not parmetis_warning)
dwarn.pop();
try
{
pdesolver->setMaxIterations(controller->getIterations());
*uold = *unew;
}
catch (Dune::PDELab::NewtonError &e){
exception = true;
std::cerr << " ERROR in Newton PDESolver: " << e << std::endl;
}
catch (Dune::ISTLError &e){
exception = true;
std::cerr << " ERROR in Linear Matrix Solver: " << e << std::endl;
}
catch (Dune::Exception &e){
std::cerr << "FATAL ERROR occured!" << std::endl;
DUNE_THROW(Dune::Exception,e);
}
if (not parmetis_warning)
dwarn.push(false);
osm->apply(t, dt, *uold, *unew);
if (not parmetis_warning)
dwarn.pop();
// controller reacts to outcome of solution
return controller->validate(exception);
*uold = *unew;
}
catch (Dune::PDELab::NewtonError &e){
exception = true;
std::cerr << " ERROR in Newton PDESolver: " << e << std::endl;
}
catch (Dune::ISTLError &e){
exception = true;
std::cerr << " ERROR in Linear Matrix Solver: " << e << std::endl;
}
catch (Dune::Exception &e){
std::cerr << "FATAL ERROR occured!" << std::endl;
DUNE_THROW(Dune::Exception,e);
}
// controller reacts to outcome of solution
step_succed = controller->validate(exception);
}
if (this->outputPolicy() == OutputPolicy::EndOfStep)
output->write_vtk_output(controller->getTime());
}
template<typename Traits>
void RichardsSimulation<Traits>::run ()
void RichardsSimulation<Traits>::markGrid ()
{
output->write_vtk_output(controller->getTime());
// TODO: split the current adapt_grid method
}
while(controller->doStep()) {
const auto time = controller->getTime();
const auto dt = controller->getDT();
if(!compute_time_step()){
continue;
}
output->write_vtk_output(controller->getTime());
template<typename Traits>
void RichardsSimulation<Traits>::adaptGrid ()
{
double time = currentTime();
double eps = inifile.get<double>("time.minTimestep")/10.;
adaptivity->adapt_grid(*grid, gv, *gfs, *param, *fboundary,
time-eps, // need "old" boundary condition TODO: check that this is right
*uold, *unew);
}
if(controller->doStep()
&& adaptivity->adapt_grid(*grid, gv, *gfs, *param, *fboundary,
time+dt, // need "old" boundary condition
*uold, *unew))
{
operator_setup(); // reset operators if grid changes
}
}
template<typename Traits>
void RichardsSimulation<Traits>::postAdaptGrid ()
{
operator_setup();
}
} // namespace Dorie
......
......@@ -13,6 +13,7 @@
#include <dune/pdelab/gridoperator/onestep.hh>
#include <dune/pdelab/newton/newton.hh>
#include "base_simulation.hh"
#include "util.hh"
#include "output.hh"
#include "adaptivity.hh"
......@@ -126,7 +127,7 @@ struct RichardsSimulationTraits : public BaseTraits
* for Richards Simulation
*/
template<typename RichardsSimulationTraits>
class RichardsSimulation
class RichardsSimulation : public SimulationBase
{
public:
using Traits = RichardsSimulationTraits;
......@@ -245,14 +246,8 @@ protected:
std::unique_ptr<AdaptivityHandler> adaptivity;
const int verbose;
public:
/**
* @brief Execute the simulation until tEnd is reached.
*/
void run ();
/*------------------------------------------------------------------------*//**
* @brief Construct basic simulation. Prepare setup for basic
* simulation.
......@@ -264,28 +259,59 @@ public:
RichardsSimulation (Dune::MPIHelper& _helper, std::shared_ptr<Grid> _grid, Dune::ParameterTree& _inifile);
/*------------------------------------------------------------------------*//**
* @brief Create the Grid Operators, Solvers, and Time Step Operator.
* This function must be called after changes to grid or GFS!
* @brief Compute a time step. Models the time step requirement of
* SimulationBase
*
* @throws Dune::Exception Fatal error occurs during computation
*
* @return The time step.
*/
void operator_setup ();
void step ();
/// Compute a time step
/** Catch errors in the Newton and the Linear Solver.
* \throw Dune::Exception Fatal error occurs during computation
* \return True if computation succeeded
/**
* @brief Mark the grid in order to improve the current model.
*/
virtual void markGrid();
/**
* @brief Adapt the grid together it every dependency of the
* grid (e.g. solution vector and grid function spaces).
*/
virtual void adaptGrid();
/**
* @brief Setup operators to fit the new grid.
*/
virtual void postAdaptGrid();
/**
* @brief Method that provides the begin time of the model.
*
* @return Begin time of the model.
*/
double beginTime() final {return controller->getBeginTime();};
/**
* @brief Method that provides the end time of the model.
*
* @return End time of the model.
*/
double endTime() final {return controller->getEndTime();};
/**
* @brief Method that provides the current time of the model.
*
* @return Current time of the model.
*/
double currentTime() final {return controller->getTime();};
protected:
/*------------------------------------------------------------------------*//**
* @brief Compute a time step. Catch errors in the Newton and the
* Linear Solver.
*
* @throws Dune::Exception Fatal error occurs during computation
*
* @return The time step.
* @brief Create the Grid Operators, Solvers, and Time Step Operator.
* This function must be called after changes to grid or GFS!
*/
bool compute_time_step ();
void operator_setup ();
};
......
......@@ -16,17 +16,18 @@ class CalculationController
{
private:
R time; //!< current time
R dt; //!< current time step
int it; //!< current number of allowed iterations for Newton solver
int verbose; //!< verbosity level
R time; //!< current time
R dt; //!< current time step
int it; //!< current number of allowed iterations for Newton solver
int verbose; //!< verbosity level
const R dtmax, dtmin; //!< time step limits
const R eps; //!< error margin
const R dtinc, dtdec; //!< time step mutliplicators
R tEnd; //!< simulation time limits
const int itmax, itmin; //!< newton iteration limits
const R tBegin; //!< simulation time limits
R tEnd; //!< simulation time limits
const int itmax, itmin; //!< newton iteration limits
const TSC& tsc; //!< Class with Function getNextTimeStamp()
const TSC& tsc; //!< Class with Function getNextTimeStamp()
public:
......@@ -43,6 +44,7 @@ public:
, eps(dtmin/10.0)
, dtinc(config.get<R>("time.timestepIncreaseFactor"))
, dtdec(config.get<R>("time.timestepDecreaseFactor"))
, tBegin(config.get<R>("time.start"))
, tEnd(config.get<R>("time.end"))
, itmax(config.get<int>("time.maxIterations"))
, itmin(config.get<int>("time.minIterations"))
......@@ -74,6 +76,10 @@ public:
inline R getTime () const { return time; }
/// Return timestep for next calculation
inline R getDT () const { return dt; }
/// Return begin time
inline R getBeginTime () const { return tBegin; }
/// Return end time
inline R getEndTime () const { return tEnd; }
/// Return maximum number of Newton iterations allowed for next calculation
inline int getIterations () const { return it; }
/// Return boolean whether new time step shall be computed
......
......@@ -213,8 +213,12 @@ public:
TestSimulation(
Dune::MPIHelper &_helper,
std::shared_ptr<Grid> _grid,
Dune::ParameterTree &_inifile)
: RichardsSimulation<Traits>(_helper, _grid, _inifile)
Dune::ParameterTree &_inifile
) : RichardsSimulation<Traits>(_helper, _grid, _inifile)
, acc(0.)
, acc_square(0.)
, wc_old(0.)
, flux_int(0.)
{ }
/// Compute the water flux into a certain domain
......@@ -256,53 +260,93 @@ public:
return std::accumulate(eta_n.begin(), eta_n.end(), 0.0);
}
/// Run the simulation
/** \return Error quantity
*/
double run_test()
void preStep() override
{
double acc{0.0}; // accumulated deviation
double acc_square{0.0}; // accumulated squared deviation
const auto time = this->controller->getTime();
const auto dt = this->controller->getDT();
wc_old = water_content(*this->uold);
// compute expected integrated flux (water content change)
flux_int = water_flux(*this->uold, time) * dt;
}
while (this->controller->doStep())
{
const auto time = this->controller->getTime();
const auto dt = this->controller->getDT();
// compute current water content
const auto wc_old = water_content(*this->uold);
// compute expected integrated flux (water content change)
const auto flux_int = water_flux(*this->uold, time) * dt;
// compute time step
if (!this->compute_time_step())
{
continue;
}
if (this->adaptivity->adapt_grid(
*this->grid, this->gv, *this->gfs, *this->param, *this->fboundary,
time, *this->uold, *this->unew))
{
// reset operators if grid changes
this->operator_setup();
}
// compute new water content
const auto wc_new = water_content(*this->unew);
const auto deviation = (wc_new - wc_old) - flux_int;
acc += deviation;
acc_square += deviation * deviation;
std::cout << "wc_old: " << wc_old << std::endl;
std::cout << "wc_new: " << wc_new << std::endl;
std::cout << "integrated flux: " << flux_int << std::endl;
std::cout << "deviation: " << deviation << std::endl;
}
void postStep() override
{
// compute new water content
const auto wc_new = water_content(*this->unew);
const auto deviation = (wc_new - wc_old) - flux_int;
acc += deviation;
acc_square += deviation * deviation;
std::cout << "wc_old: " << wc_old << std::endl;
std::cout << "wc_new: " << wc_new << std::endl;
std::cout << "integrated flux: " << flux_int << std::endl;
std::cout << "deviation: " << deviation << std::endl;
}
double run_test ()
{
acc = 0.0; // accumulated deviation
acc_square = 0.0; // accumulated squared deviation
this->run();
std::cout << "total deviation: " << acc << std::endl;
std::cout << "total normalized deviation: " << std::sqrt(acc_square) << std::endl;
return acc;
}
/// Run the simulation
/** \return Error quantity
*/
// double run_test()
// {
// double acc{0.0}; // accumulated deviation
// double acc_square{0.0}; // accumulated squared deviation
// while (this->controller->doStep())
// {
// const auto time = this->controller->getTime();
// const auto dt = this->controller->getDT();
// // compute current water content
// const auto wc_old = water_content(*this->uold);
// // compute expected integrated flux (water content change)
// const auto flux_int = water_flux(*this->uold, time) * dt;
// // compute time step
// if (!this->compute_time_step())
// {
// continue;
// }
// if (this->adaptivity->adapt_grid(
// *this->grid, this->gv, *this->gfs, *this->param, *this->fboundary,
// time, *this->uold, *this->unew))
// {
// // reset operators if grid changes
// this->operator_setup();
// }
// // compute new water content
// const auto wc_new = water_content(*this->unew);
// const auto deviation = (wc_new - wc_old) - flux_int;
// acc += deviation;
// acc_square += deviation * deviation;
// std::cout << "wc_old: " << wc_old << std::endl;
// std::cout << "wc_new: " << wc_new << std::endl;
// std::cout << "integrated flux: " << flux_int << std::endl;
// std::cout << "deviation: " << deviation << std::endl;
// }
// std::cout << "total deviation: " << acc << std::endl;
// std::cout << "total normalized deviation: " << std::sqrt(acc_square) << std::endl;
// return acc;
// }
private:
double acc; // accumulated deviation
double acc_square; // accumulated squared deviation
double wc_old;
double flux_int;
};
} // namespace Dorie
......
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