Commit a87eeb54 authored by Santiago Ospina's avatar Santiago Ospina

Move changes on richards to another branch

Signed-off-by: default avatarSantiago Ospina <saospina@hugo.iwr.uni-heidelberg.de>
parent b9997219
......@@ -5,8 +5,6 @@ namespace Dorie{
template<typename Traits>
RichardsSimulation<Traits>::RichardsSimulation (Dune::MPIHelper& _helper, std::shared_ptr<Grid> _grid, Dune::ParameterTree& _inifile) :
Dune::Dorie::SimulationBase(Traits::write_output ? OutputPolicy::EndOfStep : OutputPolicy::None,
Traits::adapt_grid ? AdaptivityPolicy::WaterFlux : AdaptivityPolicy::None),
helper(_helper), grid(_grid), inifile(_inifile), gv(grid->leafGridView()),
mbe_slop(estimate_mbe_entries<typename MBE::size_type>(Traits::dim,Traits::GridGeometryType)),
mbe_tlop(1),
......@@ -117,78 +115,64 @@ void RichardsSimulation<Traits>::operator_setup ()
}
template<typename Traits>
void RichardsSimulation<Traits>::step ()
bool RichardsSimulation<Traits>::compute_time_step ()
{
bool step_succed = false;
while(controller->doStep() & not step_succed)
{
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
{
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
step_succed = 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);
}
if (this->outputPolicy() == OutputPolicy::EndOfStep)
writeData();
}
template<typename Traits>
void RichardsSimulation<Traits>::markGrid ()
{
// TODO: split the current adapt_grid method
// controller reacts to outcome of solution
return controller->validate(exception);
}
template<typename Traits>
void RichardsSimulation<Traits>::adaptGrid ()
void RichardsSimulation<Traits>::run ()
{
double time = currentTime();
double eps = inifile.get<double>("time.minTimestep")/10.;
if constexpr (Traits::adapt_grid)
adaptivity->adapt_grid(*grid, gv, *gfs, *param, *fboundary,
time-eps, // need "old" boundary condition TODO: check that this is right
*uold, *unew);
else
DUNE_THROW(NotImplemented,"Dorie does not implement adapticivity for this simulation object!");
}
output->write_vtk_output(controller->getTime());
template<typename Traits>
void RichardsSimulation<Traits>::postAdaptGrid ()
{
operator_setup();
}
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>::writeData ()
{
output->write_vtk_output(controller->getTime());
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
}
}
}
} // namespace Dorie
......
......@@ -13,7 +13,6 @@
#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"
......@@ -127,7 +126,7 @@ struct RichardsSimulationTraits : public BaseTraits
* for Richards Simulation
*/
template<typename RichardsSimulationTraits>
class RichardsSimulation : public SimulationBase
class RichardsSimulation
{
public:
using Traits = RichardsSimulationTraits;
......@@ -246,8 +245,14 @@ 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.
......@@ -259,74 +264,28 @@ public:
RichardsSimulation (Dune::MPIHelper& _helper, std::shared_ptr<Grid> _grid, Dune::ParameterTree& _inifile);
/*------------------------------------------------------------------------*//**
* @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 step () override;
/**
* @brief Writes the data.
*/
void virtual writeData() override;
/**
* @brief Mark the grid in order to improve the current model.
*/
virtual void markGrid() override;
/**
* @brief Adapt the grid together it every dependency of the
* grid (e.g. solution vector and grid function spaces).
* @brief Create the Grid Operators, Solvers, and Time Step Operator.
* This function must be called after changes to grid or GFS!
*/
virtual void adaptGrid() override;
void operator_setup ();
/**
* @brief Setup operators to fit the new grid.
/// 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
*/
virtual void postAdaptGrid() override;
/**
* @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();};
/**
* @brief Suggest a time step to the model.
*
* @param[in] suggestion for the internal time step of the model.
*/
void suggestTimestep(double dt) final
{/*TODO: ensures that controller will not provide a dt higher than suggested dt*/};
protected:
/*------------------------------------------------------------------------*//**
* @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. Catch errors in the Newton and the
* Linear Solver.
*
* @throws Dune::Exception Fatal error occurs during computation
*
* @return The time step.
*/
void operator_setup ();
bool compute_time_step ();
};
......
......@@ -16,18 +16,17 @@ 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
const R tBegin; //!< simulation time limits
R tEnd; //!< simulation time limits
const int itmax, itmin; //!< newton iteration 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:
......@@ -44,7 +43,6 @@ 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"))
......@@ -76,10 +74,6 @@ 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,12 +213,8 @@ public:
TestSimulation(
Dune::MPIHelper &_helper,
std::shared_ptr<Grid> _grid,
Dune::ParameterTree &_inifile
) : RichardsSimulation<Traits>(_helper, _grid, _inifile)
, acc(0.)
, acc_square(0.)
, wc_old(0.)
, flux_int(0.)
Dune::ParameterTree &_inifile)
: RichardsSimulation<Traits>(_helper, _grid, _inifile)
{ }
/// Compute the water flux into a certain domain
......@@ -260,47 +256,53 @@ public:
return std::accumulate(eta_n.begin(), eta_n.end(), 0.0);
}
void preStep() override
{
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;
}
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 ()
/// Run the simulation
/** \return Error quantity
*/
double run_test()
{
acc = 0.0; // accumulated deviation
acc_square = 0.0; // accumulated squared deviation
double acc{0.0}; // accumulated deviation
double acc_square{0.0}; // accumulated squared deviation
this->run();
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