Commit 25baeadc authored by Lukas Riedel's avatar Lukas Riedel 🎧
Browse files

refactored Adaptivity: added base and factory classes, added new to Simulation

parent a2bf27ec
......@@ -30,18 +30,46 @@ public:
using ESTGO = Dune::PDELab::GridOperator<typename SimTraits::GFS::Type,typename GFS::Type,ESTLOP,MBE,RF,RF,RF,NoTrafo,NoTrafo>;
};
/// Adaptivity base class. Does nothing.
template<typename Traits>
class AdaptivityHandlerBase
{
private:
using Grid = typename Traits::Grid;
using GV = typename Traits::GV;
using SimTraits = Dune::Dorie::SimulationTraits<Traits>;
using GFS = typename SimTraits::GFS::Type;
using Param = typename SimTraits::Parameters;
using U = typename SimTraits::U;
public:
AdaptivityHandlerBase () = default;
virtual ~AdaptivityHandlerBase () = default;
/// Do nothing.
/** \return Boolean if operators have to be updated
*/
virtual bool adapt_grid (Grid& grid, GV& gv, GFS& gfs, Param& param, U& uold, U& unew) { return false; }
};
/// Specialized SimulationBase class providing functions and members for adaptive grid refinement
template<typename Traits>
class AdaptivityHandler : public virtual SimulationBase<Traits>
class AdaptivityHandler : public AdaptivityHandlerBase<Traits>
{
private:
static constexpr int order = Traits::fem_order;
using RF = typename Traits::RF;
using Grid = typename Traits::Grid;
using GV = typename Traits::GV;
using SimTraits = Dune::Dorie::SimulationTraits<Traits>;
using AdaptTraits = Dune::Dorie::AdaptivityTraits<Traits,SimTraits>;
using GFS = typename SimTraits::GFS::Type;
using Param = typename SimTraits::Parameters;
using U = typename SimTraits::U;
using AGFS = typename AdaptTraits::GFS::Type;
using ESTLOP = typename AdaptTraits::ESTLOP;
using MBE = typename AdaptTraits::MBE;
......@@ -50,47 +78,52 @@ private:
private:
const Dune::ParameterTree& inifile; //!< Parameter file parser
const int maxLevel; //!< Maximum local grid refinement level
const int minLevel; //!< Minimum local grid refinement level
const std::string strategy; //!< Refinement strategy
const double alpha; //!< Refinement factor
const double beta; //!< Coarsening factor
const float adaptivityThreshold; //!< Global error threshold below which no refinement is applied
const int verbose; //!< output verbosity of this object
public:
/// Initialize members from config file parameters.
/**
* \param _helper Dune MPI instance controller
* \param _grid Shared pointer to the grid
* \param _inifile Dune parameter file parser
*/
AdaptivityHandler (Dune::MPIHelper& _helper, std::shared_ptr<Grid> _grid, Dune::ParameterTree& _inifile)
: SimulationBase<Traits>(_helper,_grid,_inifile),
maxLevel(_inifile.get<int>("adaptivity.maxLevel")),
minLevel(_inifile.get<int>("adaptivity.minLevel")),
strategy(_inifile.get<std::string>("adaptivity.markingStrategy")),
alpha(_inifile.get<float>("adaptivity.refinementFraction")),
beta(_inifile.get<float>("adaptivity.coarseningFraction")),
adaptivityThreshold(_inifile.get<float>("adaptivity.threshold"))
AdaptivityHandler (Dune::ParameterTree& _inifile, Grid& grid)
: AdaptivityHandlerBase<Traits>(),
inifile(_inifile),
maxLevel(inifile.get<int>("adaptivity.maxLevel")),
minLevel(inifile.get<int>("adaptivity.minLevel")),
strategy(inifile.get<std::string>("adaptivity.markingStrategy")),
alpha(inifile.get<float>("adaptivity.refinementFraction")),
beta(inifile.get<float>("adaptivity.coarseningFraction")),
adaptivityThreshold(inifile.get<float>("adaptivity.threshold")),
verbose(inifile.get<int>("output.verbose"))
{
if(strategy!="elementFraction" && strategy!="errorFraction" && strategy!="threshold" && strategy!="targetTolerance")
DUNE_THROW(Dune::IOError,"Valid values of adaptivity.markingStrategy: elementFraction, errorFraction, threshold, targetTolerance");
// QkDG space cannot handle triagles arising from closure of adapted grid
if(Traits::GridGeometryType == Dune::GeometryType::BasicType::cube){
this->grid->setClosureType(Traits::Grid::ClosureType::NONE);
grid.setClosureType(Traits::Grid::ClosureType::NONE);
}
}
/// Adapt the grid according to the estimated error and the strategy applied
void adapt_grid ()
/// Adapt the grid according to the estimated flux error and the strategy applied
/**
* \return Boolean if grid operators have to be set up again (true if grid changed)
*/
bool adapt_grid (Grid& grid, GV& gv, GFS& gfs, Param& param, U& uold, U& unew) override
{
Dune::Timer timer2;
float t2;
double maxeta(0.0);
int multiRefinementCounter = 0;
if(this->verbose>1){
if(verbose>1){
std::cout << "GRID ADAPTION [" << strategy << "]:" << std::endl;
}
......@@ -98,7 +131,7 @@ public:
Dune::Timer timer3;
float t_setup,t_est,t_sqrt,t_strtgy,t_mark,t_adapt;
if(this->verbose>1 /*&& strategy == "targetTolerance"*/){
if(verbose>1 /*&& strategy == "targetTolerance"*/){
std::cout << " Refinement Step " << multiRefinementCounter << ": ";
}
......@@ -106,16 +139,16 @@ public:
double eta_beta(0.0);
// Compute error estimate eta
AGFS p0gfs = AdaptTraits::GFS::create(this->gv);
ESTLOP estlop(this->inifile,*this->param);
AGFS p0gfs = AdaptTraits::GFS::create(gv);
ESTLOP estlop(inifile,param);
MBE mbe(0.0);
ESTGO estgo(*this->gfs,p0gfs,estlop,mbe);
ESTGO estgo(gfs,p0gfs,estlop,mbe);
U0 eta(p0gfs,0.0);
t_setup = timer3.elapsed();
timer3.reset();
estgo.residual(*this->unew,eta);
estgo.residual(uold,eta);
t_est = timer3.elapsed();
timer3.reset();
......@@ -126,7 +159,7 @@ public:
native(eta)[i] = sqrt(native(eta)[i]); // eta contains squares
if (native(eta)[i] > maxeta) maxeta = native(eta)[i];
}
if (this->verbose>1) {
if (verbose>1) {
std::cout << "Max Local Error: " << maxeta << " " << std::endl;
if (maxeta < adaptivityThreshold) {
std::cout << " Max local error smaller than threshold " << adaptivityThreshold << ". Skipping grid refinement." << std::endl;
......@@ -138,16 +171,16 @@ public:
// Stop adaptivity if target error is reached
if ((maxeta < alpha) && strategy == "targetTolerance") {
if(this->verbose>1)
if(verbose>1)
std::cout << " Target error " << alpha << " reached!" << std::endl;
break;
}
// Apply marking strategy
if (strategy == "elementFraction")
Dune::PDELab::element_fraction( eta, alpha, beta, eta_alpha, eta_beta, this->verbose-1 );
Dune::PDELab::element_fraction( eta, alpha, beta, eta_alpha, eta_beta, verbose-1 );
else if (strategy == "errorFraction")
Dune::PDELab::error_fraction( eta, alpha, beta, eta_alpha, eta_beta, this->verbose-1 );
Dune::PDELab::error_fraction( eta, alpha, beta, eta_alpha, eta_beta, verbose-1 );
else{ //((strategy == "threshold") || (strategy == "targetTolerance")) {
eta_alpha = alpha;
eta_beta = beta; }
......@@ -159,12 +192,12 @@ public:
if (maxeta < adaptivityThreshold)
eta_alpha = maxeta + 1; // Only refine elements with error > maxerror + 1
Dune::PDELab::mark_grid( *this->grid, eta, eta_alpha, eta_beta, minLevel, maxLevel, this->verbose-1);
Dune::PDELab::mark_grid( grid, eta, eta_alpha, eta_beta, minLevel, maxLevel, verbose-1);
t_mark = timer3.elapsed();
timer3.reset();
Dune::PDELab::adapt_grid( *this->grid, *this->gfs, *this->uold, *this->unew, (this->order) * 2 );
Dune::PDELab::adapt_grid( grid, gfs, uold, unew, (order) * 2 );
t_adapt = timer3.elapsed();
timer3.reset();
......@@ -184,7 +217,7 @@ public:
} while (strategy == "targetTolerance");
if (this->verbose>1) {
if (verbose>1) {
t2 = timer2.elapsed();
ios_base_all_saver restore(std::cout);
std::cout << std::setprecision(4) << std::scientific;
......@@ -192,10 +225,52 @@ public:
std::cout << "::: total adaptation time" << std::setw(12) << t2 << std::endl;
}
// Update Grid Operators
this->operator_setup();
return true;
}
};
/// Factory class for creating an AdaptivityHandler based on Traits
/** \tparam Traits Dune::Dorie::BaseTraits defining the Simulation
* The constexpr boolean adapt_grid in Dune::Dorie::BaseTraits will
* define which create() function is enabled at compile time.
*/
template<typename Traits>
struct AdaptivityHandlerFactory
{
private:
static constexpr bool enabled = Traits::adapt_grid;
using Grid = typename Traits::Grid;
Dune::ParameterTree& inifile;
Grid& grid;
public:
AdaptivityHandlerFactory (Dune::ParameterTree& _inifile, Grid& _grid) :
inifile(_inifile),
grid(_grid)
{
static_assert(std::is_same<Grid,Dune::UGGrid<2>>::value ||
std::is_same<Grid,Dune::UGGrid<3>>::value,
"Adaptivity is only implemented for UGGrid!");
}
/// Create a dummy AdaptivityHandler
template<bool active = enabled>
std::unique_ptr<AdaptivityHandlerBase<Traits>> create (typename std::enable_if_t<!active,int> = 0)
{
return std::make_unique<AdaptivityHandlerBase<Traits>>();
}
/// Create a true AdaptivityHandler
template<bool active = enabled>
std::unique_ptr<AdaptivityHandlerBase<Traits>> create (typename std::enable_if_t<active,int> = 0)
{
std::unique_ptr<AdaptivityHandlerBase<Traits>> p;
p = std::make_unique<AdaptivityHandler<Traits>>(inifile,grid);
return p;
}
};
} // namespace Dorie
......
......@@ -68,8 +68,8 @@
#include "util.hh"
#include "output.hh"
#include "simulation_base.hh"
#include "adaptivity.hh"
#include "simulation_base.hh"
#include "knofu.hh"
#include "simulation.hh"
......@@ -98,9 +98,10 @@ int main(int argc, char **argv)
const std::string outputPath = inifile.get<std::string>("output.outputPath");
const bool adaptivity = inifile.get<bool>("adaptivity.useAdaptivity");
auto yasp_grid = Dune::Dorie::build_grid_cube<Dune::YaspGrid<2>>(inifile,helper);
auto ug_grid = Dune::Dorie::build_grid_cube<Dune::UGGrid<2>>(inifile,helper);
Dune::Dorie::SimulationBase<Dune::Dorie::BaseTraits<Dune::UGGrid,Dune::GeometryType::BasicType::cube,2,1,true>> sim(helper,ug_grid,inifile);
Dune::Dorie::SimulationBase<Dune::Dorie::BaseTraits<Dune::UGGrid,Dune::GeometryType::BasicType::cube,2,1,true,false>> sim(helper,ug_grid,inifile);
sim.run();
/*
......
......@@ -74,6 +74,7 @@ protected:
std::unique_ptr<U> unew;
std::unique_ptr<OutputWriterBase<Traits>> output;
std::unique_ptr<AdaptivityHandlerBase<Traits>> adaptivity;
const int verbose;
......@@ -84,7 +85,12 @@ public:
{
output->write_vtk_output(controller->getTime());
while(controller->doStep()) {
compute_time_step();
if(!compute_time_step()){
continue;
}
if(adaptivity->adapt_grid(*grid, gv, *gfs, *param, *uold, *unew)){ // reset operators if grid changes
operator_setup();
}
output->write_vtk_output(controller->getTime());
}
}
......@@ -134,10 +140,13 @@ public:
// --- Operator Setup --- //
operator_setup();
// --- Utility Class Setup --- //
output = Dune::Dorie::OutputWriterFactory<Traits>::create(inifile,gv,*gfs,*param,*unew);
auto adaptivity_fac = Dune::Dorie::AdaptivityHandlerFactory<Traits>(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;
output = Dune::Dorie::OutputWriterFactory<Traits>::create(inifile,gv,*gfs,*param,*unew);
}
/// Create the Grid Operators, Solvers, and Time Step Operator
......
......@@ -92,12 +92,13 @@ public:
};
/// Traits struct defining basic types based on grid an geometry
template< template<int> class GridType, Dune::GeometryType::BasicType GeometryType, int dimension, int order , bool output>
template< template<int> class GridType, Dune::GeometryType::BasicType GeometryType, int dimension, int order , bool output, bool adaptivity>
struct BaseTraits
{
static constexpr int dim = dimension;
static constexpr int fem_order = order;
static constexpr bool write_output = output;
static constexpr bool adapt_grid = adaptivity;
static constexpr Dune::GeometryType::BasicType GridGeometryType = GeometryType;
using RF = double;
......
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