Commit debf828e authored by Lukas Riedel's avatar Lukas Riedel

Merge branch '95-modify-richardssimulation-to-have-the-base-class-94' into 'master'

Resolve "Modify RichardsSimulation to have the base class #94."

Closes #90 and #95

See merge request !89
parents cd7e2bc7 5286f3fd
......@@ -49,11 +49,19 @@
solution vector.
* `GradientFluxAdapter` was reimplemented and renamed `WaterFluxAdapter`.
* Every grid function adapter has its own file and are gathered in the
subdirectory [dune/dorie/solver/adapters](dune/dorie/solver/adapters)
subdirectory [dune/dorie/solver/adapters](dune/dorie/solver/adapters).
* `VTKAdapters` are now managed with `shared_ptr` instead of references.
* `OutputWriter` class is deprecated in favor of an minimal extension of the
usual `VTKSequenceWriter` for grid functions called
`GridFunctionVTKSequenceWriter`.
* `RichardsSimulation` now models the `SimulationBase` class (See details in
https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/issues/94). For this,
several methods were renamed with basically no change in behaviour.
* `RichardsSimulation` and `TimeController` are able to receive a suggestion for
the next timestep.
* Static const `bool`s for adaptivity and output writing in `BaseTraits` were
removed in favor of run-time policies (See details in
https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/issues/94).
* Swich to the stable release branch `releases/2.6` of `dune-testtools`.
* The `build` CI stage now explicitly builds system tests and unit tests. The
job `build:main` has been removed. The `build:debug` jobs now build the
......
......@@ -43,19 +43,15 @@ using Sim = Dune::Dorie::RichardsSimulation<Traits>;
template<int dim, int order>
using Simplex = Dune::Dorie::RichardsSimulationTraits<Dune::Dorie::BaseTraits<Dune::UGGrid<dim>,
Dune::GeometryType::BasicType::simplex,true,false>,order>;
template<int dim, int order>
using SimplexAdaptive = Dune::Dorie::RichardsSimulationTraits<Dune::Dorie::BaseTraits<Dune::UGGrid<dim>,
Dune::GeometryType::BasicType::simplex,true,true>,order>;
Dune::GeometryType::BasicType::simplex>,order>;
template<int dim, int order>
using Cube = Dune::Dorie::RichardsSimulationTraits<Dune::Dorie::BaseTraits<Dune::YaspGrid<dim>,
Dune::GeometryType::BasicType::cube,true,false>,order>;
Dune::GeometryType::BasicType::cube>,order>;
template<int dim, int order>
using CubeAdaptive = Dune::Dorie::RichardsSimulationTraits<Dune::Dorie::BaseTraits<Dune::UGGrid<dim>,
Dune::GeometryType::BasicType::cube,true,true>,order>;
Dune::GeometryType::BasicType::cube>,order>;
/// Main Program Function: Initialize parameter file reader, build grid and call Richards Solver.
/** As simplex and rectangular grids demand different FiniteElementMaps, the program calls different functions for these tasks.
......@@ -89,21 +85,21 @@ int main(int argc, char** argv)
Dune::ParameterTreeParser ptreeparser;
ptreeparser.readINITree(inifilename,inifile);
// Allow for a debugger as gdb or lldb to hook into the process, even if run in parallel
// (by manually setting the variable i to a nonzero value)
bool debugMode = inifile.get<bool>("misc.debugMode");
if (debugMode)
{
int i = 0;
char hostname[256];
gethostname(hostname, sizeof(hostname));
if(helper.rank()==0)
std::cout << "Debug mode activated. Use your debugger to set the variable 'i' to a value > 0 in each process" << std::endl;
printf("PID %d on %s ready for attachment\n", getpid(), hostname);
fflush(stdout);
while (0 == i)
sleep(5);
}
// Allow for a debugger as gdb or lldb to hook into the process, even if run in parallel
// (by manually setting the variable i to a nonzero value)
bool debugMode = inifile.get<bool>("misc.debugMode");
if (debugMode)
{
int i = 0;
char hostname[256];
gethostname(hostname, sizeof(hostname));
if(helper.rank()==0)
std::cout << "Debug mode activated. Use your debugger to set the variable 'i' to a value > 0 in each process" << std::endl;
printf("PID %d on %s ready for attachment\n", getpid(), hostname);
fflush(stdout);
while (0 == i)
sleep(5);
}
// Read necessary variables
const std::string gtype = inifile.get<std::string>("grid.gridType");
......@@ -113,11 +109,17 @@ int main(int argc, char** argv)
const std::string outputPath = inifile.get<std::string>("output.outputPath");
const bool adaptivity = inifile.get<bool>("adaptivity.useAdaptivity");
Dune::Dorie::AdaptivityPolicy adapt_policy =
Dune::Dorie::AdaptivityPolicy::None;
if (adaptivity) {
adapt_policy = Dune::Dorie::AdaptivityPolicy::WaterFlux;
}
// Attempt to create output directory
mkdir(outputPath.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
int result = access(outputPath.c_str(), W_OK);
if (result != 0)
DUNE_THROW(Dune::IOError,"Output folder " << outputPath << " not writable");
int result = access(outputPath.c_str(), W_OK);
if (result != 0)
DUNE_THROW(Dune::IOError,"Output folder " << outputPath << " not writable");
if (helper.rank()==0){
std::cout << "INPUT FILE: " << inifilename << std::endl;
......@@ -133,47 +135,27 @@ int main(int argc, char** argv)
{
if (gtype == "gmsh"){
auto grid = Dune::Dorie::build_grid_gmsh<Dune::UGGrid<2>>(inifile,helper);
if(adaptivity){
switch(FEorder){
case 1:{
Sim<SimplexAdaptive<2,1>> sim(helper,grid,inifile);
sim.run();
break;
}
case 2:{
Sim<SimplexAdaptive<2,2>> sim(helper,grid,inifile);
sim.run();
break;
}
case 3:{
Sim<SimplexAdaptive<2,3>> sim(helper,grid,inifile);
sim.run();
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
switch(FEorder){
case 1:{
Sim<Simplex<2,1>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
}
else{ // no adaptivity
switch(FEorder){
case 1:{
Sim<Simplex<2,1>> sim(helper,grid,inifile);
sim.run();
break;
}
case 2:{
Sim<Simplex<2,2>> sim(helper,grid,inifile);
sim.run();
break;
}
case 3:{
Sim<Simplex<2,3>> sim(helper,grid,inifile);
sim.run();
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
case 2:{
Sim<Simplex<2,2>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 3:{
Sim<Simplex<2,3>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
}
}
else if (gtype == "rectangular"){
......@@ -182,16 +164,19 @@ int main(int argc, char** argv)
switch(FEorder){
case 1:{
Sim<CubeAdaptive<2,1>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 2:{
Sim<CubeAdaptive<2,2>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 3:{
Sim<CubeAdaptive<2,3>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
......@@ -230,47 +215,27 @@ int main(int argc, char** argv)
{
if (gtype == "gmsh"){
auto grid = Dune::Dorie::build_grid_gmsh<Dune::UGGrid<3>>(inifile,helper);
if(adaptivity){
switch(FEorder){
case 1:{
Sim<SimplexAdaptive<3,1>> sim(helper,grid,inifile);
sim.run();
break;
}
case 2:{
Sim<SimplexAdaptive<3,2>> sim(helper,grid,inifile);
sim.run();
break;
}
case 3:{
Sim<SimplexAdaptive<3,3>> sim(helper,grid,inifile);
sim.run();
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
switch(FEorder){
case 1:{
Sim<Simplex<3,1>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
}
else{ // no adaptivity
switch(FEorder){
case 1:{
Sim<Simplex<3,1>> sim(helper,grid,inifile);
sim.run();
break;
}
case 2:{
Sim<Simplex<3,2>> sim(helper,grid,inifile);
sim.run();
break;
}
case 3:{
Sim<Simplex<3,3>> sim(helper,grid,inifile);
sim.run();
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
case 2:{
Sim<Simplex<3,2>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 3:{
Sim<Simplex<3,3>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
}
}
else if (gtype == "rectangular"){
......@@ -279,16 +244,19 @@ int main(int argc, char** argv)
switch(FEorder){
case 1:{
Sim<CubeAdaptive<3,1>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 2:{
Sim<CubeAdaptive<3,2>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 3:{
Sim<CubeAdaptive<3,3>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
......
......@@ -8,9 +8,8 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube,true,true>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex,true,true>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex,true,false>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,1>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,9 +8,8 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube,true,true>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex,true,true>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex,true,false>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,2>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,9 +8,8 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube,true,true>,3>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex,true,true>,3>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex,true,false>,3>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,3>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,3>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,9 +8,8 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube,true,true>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex,true,true>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex,true,false>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,1>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,9 +8,8 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube,true,true>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex,true,true>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex,true,false>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,2>>;
} // namespace Dorie
} // namespace Dune
......@@ -8,9 +8,8 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube,true,true>,3>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex,true,true>,3>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex,true,false>,3>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,3>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,3>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,7 +8,7 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube,true,false>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,1>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,7 +8,7 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube,true,false>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,2>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,7 +8,7 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube,true,false>,3>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,3>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,7 +8,7 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube,true,false>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,1>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,7 +8,7 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube,true,false>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,2>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,7 +8,7 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube,true,false>,3>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,3>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -38,8 +38,6 @@ public:
virtual void mark_grid (Grid& grid, GV& gv, GFS& gfs, const Param& param, const Boundary& boundary, const RF time, U& u) { }
/// Do nothing.
/** \return Boolean if operators have to be updated
*/
virtual void adapt_grid (Grid& grid, GFS& gfs, U& u) { }
};
......@@ -257,7 +255,27 @@ template<typename Traits, typename GFS, typename Param, typename Boundary, typen
struct AdaptivityHandlerFactory
{
private:
static constexpr bool enabled = Traits::adapt_grid;
/*---------------------------------------------------------------------*//**
* @brief Helper to assert wheter a grid is adaptable.
*
* @tparam GridType Type of the grid
*/
template<class GridType>
struct isAdaptable : std::false_type {};
// The only known adaptable grid for dorie is UGGrid
/*---------------------------------------------------------------------*//**
* @brief Helper to assert wheter a grid is adaptable. UGGrid
* specialization
*
* @tparam dim Dimension of the UGGrid
*/
template<int dim>
struct isAdaptable<Dune::UGGrid<dim>> : std::true_type {};
using Grid = typename Traits::Grid;
using AHB = AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U,order>;
......@@ -267,10 +285,14 @@ private:
public:
/// Create the factory, taking the parameters for building an AdaptivityHandler
/**
* \param _inifile Parameter file parser
* \param _grid Grid to adapt (reference is not saved)
static constexpr bool enabled = isAdaptable<Grid>::value;
/*---------------------------------------------------------------------*//**
* @brief Create the factory, taking the parameters for building an
* AdaptivityHandler
*
* @param _inifile Parameter file parser
* @param _grid Grid to adapt (reference is not saved)
*/
AdaptivityHandlerFactory (Dune::ParameterTree& _inifile, Grid& _grid) :
inifile(_inifile),
......@@ -290,9 +312,6 @@ public:
template<bool active = enabled>
std::enable_if_t<active,std::unique_ptr<AHB>> create ()
{
static_assert(std::is_same<Grid,Dune::UGGrid<2>>::value ||
std::is_same<Grid,Dune::UGGrid<3>>::value,
"Adaptivity is only implemented for UGGrid!");
std::unique_ptr<AHB> p;
p = std::make_unique<AdaptivityHandler<Traits,GFS,Param,Boundary,U,order>>(inifile,grid);
return p;
......
......@@ -61,7 +61,7 @@ RichardsSimulation<Traits>::RichardsSimulation (Dune::MPIHelper& _helper, std::s
operator_setup();
// --- Utility Class Setup --- //
if constexpr (Traits::write_output)
if (output_policy() != OutputPolicy::None)
{
const int subsamling_lvl = _inifile.get<int>("output.subsamplingLevel", 0);
const auto subsamling_intervals = Dune::refinementLevels(subsamling_lvl);
......@@ -74,15 +74,18 @@ RichardsSimulation<Traits>::RichardsSimulation (Dune::MPIHelper& _helper, std::s
"./");
}
AdaptivityHandlerFactory adaptivity_fac(inifile,*grid);
adaptivity = adaptivity_fac.create();
if constexpr (AdaptivityHandlerFactory::enabled)
{
AdaptivityHandlerFactory adaptivity_fac(inifile,*grid);
adaptivity = adaptivity_fac.create();
}
if(verbose>1 && helper.rank()==0)
std::cout << "::: setup time " << std::setw(12) << std::setprecision(4) << std::scientific << timer.elapsed() << std::endl;
}
template<typename Traits>
void RichardsSimulation<Traits>::operator_setup ()
void RichardsSimulation<Traits>::operator_setup()
{
Dune::PDELab::constraints(*this->gfs,*this->cc,false);
// --- Grid Operators ---
......@@ -128,77 +131,78 @@ 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_succeed = false;
while(not step_succeed)
{
std::shared_ptr<U> unext = std::make_shared<U>(*u);
pdesolver->setMaxIterations(controller->getIterations());
if (not parmetis_warning)
dwarn.push(false);
osm->apply(t, dt, *u, *unext);
if (not parmetis_warning)
dwarn.pop();
u = unext;
}
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);
const RF t = controller->getTime();
const RF dt = controller->getDT();
bool exception = false;
const bool solver_warnings = verbose > 0 && helper.rank() == 0 ?
true : false;
try
{
std::shared_ptr<U> unext = std::make_shared<U>(*u);
pdesolver->setMaxIterations(controller->getIterations());
if (not solver_warnings)
dwarn.push(false);
osm->apply(time_before, dt, *u, *unext);
if (not solver_warnings)
dwarn.pop();
u = unext;
}
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);
}
// saving last times for adaptivity
time_before = t;
dt_before = dt;
// controller reacts to outcome of solution
step_succeed = controller->validate(exception);
}
// controller reacts to outcome of solution
return controller->validate(exception);
if (this->output_policy() == OutputPolicy::EndOfStep)
write_data();
}
template<typename Traits>
void RichardsSimulation<Traits>::run ()
void RichardsSimulation<Traits>::mark_grid()