Commit 93f81a50 authored by Santiago Ospina's avatar Santiago Ospina

split in grid adaptivity interface

parent 2001ba22
......@@ -34,10 +34,13 @@ public:
AdaptivityHandlerBase () = default;
virtual ~AdaptivityHandlerBase () = default;
/// Do nothing.
virtual void mark_grid (Grid& grid, GV& gv, GFS& gfs, const Param& param, const Boundary& boundary, const RF time, U& uold, U& unew) { }
/// Do nothing.
/** \return Boolean if operators have to be updated
*/
virtual bool adapt_grid (Grid& grid, GV& gv, GFS& gfs, const Param& param, const Boundary& boundary, const RF time, U& uold, U& unew) { return false; }
virtual void adapt_grid (Grid& grid, GFS& gfs, U& uold, U& unew) { }
};
/// Specialized SimulationBase class providing functions and members for adaptive grid refinement
......@@ -79,6 +82,7 @@ private:
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
float total_time; //!< accumulated time of grid marking and pgrid
public:
......@@ -96,10 +100,11 @@ public:
alpha(inifile.get<float>("adaptivity.refinementFraction")),
beta(inifile.get<float>("adaptivity.coarseningFraction")),
adaptivityThreshold(inifile.get<float>("adaptivity.threshold")),
verbose(inifile.get<int>("output.verbose"))
verbose(inifile.get<int>("output.verbose")),
total_time(0.)
{
if(strategy!="elementFraction" && strategy!="errorFraction" && strategy!="threshold" && strategy!="targetTolerance")
DUNE_THROW(Dune::IOError,"Valid values of adaptivity.markingStrategy: elementFraction, errorFraction, threshold, targetTolerance");
if(strategy!="elementFraction" && strategy!="errorFraction" && strategy!="threshold")
DUNE_THROW(Dune::IOError,"Valid values of adaptivity.markingStrategy: elementFraction, errorFraction, threshold");
// QkDG space cannot handle triagles arising from closure of adapted grid
if(Traits::GridGeometryType == Dune::GeometryType::BasicType::cube){
......@@ -117,7 +122,7 @@ public:
* \param unew Solution vector
* \return Boolean if grid operators have to be set up again (true if grid changed)
*/
bool adapt_grid (
void mark_grid (
Grid& grid,
GV& gv,
GFS& gfs,
......@@ -128,7 +133,6 @@ public:
U& unew) override
{
Dune::Timer timer2;
float t2;
double maxeta(0.0);
int multiRefinementCounter = 0;
......@@ -136,107 +140,107 @@ public:
std::cout << "GRID ADAPTION [" << strategy << "]:" << std::endl;
}
do { // Loop for multiple refinements
Dune::Timer timer3;
float t_setup, t_sqrt, t_est, t_strtgy, t_mark, t_adapt;
Dune::Timer timer3;
float t_setup, t_sqrt, t_est, t_strtgy, t_mark;
if(verbose>1)
std::cout << " Refinement Step " << multiRefinementCounter << ": ";
double eta_alpha(0.0);
double eta_beta(0.0);
// set up helper GFS
AGFS p0gfs = AGFSHelper::create(gv);
ESTLOP estlop(inifile, param, boundary);
estlop.setTime(time);
MBE mbe(0);
ESTGO estgo(gfs,p0gfs,estlop,mbe);
U0 eta(p0gfs,0.0);
t_setup = timer3.elapsed();
timer3.reset();
// Compute error estimate eta
estgo.residual(uold,eta);
t_est = timer3.elapsed();
timer3.reset();
// unsquare errors
using Dune::PDELab::Backend::native;
auto& eta_nat = native(eta);
std::transform(eta_nat.begin(), eta_nat.end(), eta_nat.begin(),
[](const auto val) { return std::sqrt(val); }
);
// compute largest error
maxeta = eta.infinity_norm();
if (verbose>1)
std::cout << "Largest Local Error: " << maxeta << " " << std::endl;
t_sqrt = timer3.elapsed();
timer3.reset();
// Apply marking strategy
if (strategy == "elementFraction")
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, verbose-1);
else{ //((strategy == "threshold") || (strategy == "targetTolerance")) {
eta_alpha = alpha;
eta_beta = beta; }
t_strtgy = timer3.elapsed();
timer3.reset();
// Skip refinement if threshold is met, but still allow coarsening
if (maxeta < adaptivityThreshold) {
eta_alpha = maxeta + 1; // Only refine elements with error > maxerror + 1
if (verbose>1)
std::cout << " Max local error smaller than threshold " << adaptivityThreshold << ". Skipping grid refinement." << std::endl;
}
if(verbose>1)
std::cout << " Refinement Step " << multiRefinementCounter << ": ";
Dune::PDELab::mark_grid(grid, eta, eta_alpha, eta_beta, minLevel, maxLevel, verbose-1);
double eta_alpha(0.0);
double eta_beta(0.0);
t_mark = timer3.elapsed();
timer3.reset();
// set up helper GFS
AGFS p0gfs = AGFSHelper::create(gv);
ESTLOP estlop(inifile, param, boundary);
estlop.setTime(time);
MBE mbe(0);
ESTGO estgo(gfs,p0gfs,estlop,mbe);
U0 eta(p0gfs,0.0);
if(this->verbose>3){
ios_base_all_saver restore(std::cout);
std::cout << std::setprecision(4) << std::scientific;
std::cout << "::: setup time" << std::setw(12) << t_setup << std::endl;
std::cout << "::: error estimation time" << std::setw(12) << t_est << std::endl;
std::cout << "::: unsquaring time" << std::setw(12) << t_sqrt << std::endl;
std::cout << "::: strategy application time" << std::setw(12) << t_strtgy << std::endl;
std::cout << "::: grid marking time" << std::setw(12) << t_mark << std::endl;
}
t_setup = timer3.elapsed();
timer3.reset();
total_time += timer2.elapsed();
}
// Compute error estimate eta
estgo.residual(uold,eta);
void adapt_grid (
Grid& grid,
GFS& gfs,
U& uold,
U& unew) override
{
Dune::Timer timer;
t_est = timer3.elapsed();
timer3.reset();
Dune::PDELab::adapt_grid(grid, gfs, uold, unew, 2*order);
// unsquare errors
using Dune::PDELab::Backend::native;
auto& eta_nat = native(eta);
std::transform(eta_nat.begin(), eta_nat.end(), eta_nat.begin(),
[](const auto val) { return std::sqrt(val); }
);
float t_adapt = timer.elapsed();
total_time += t_adapt;
// compute largest error
maxeta = eta.infinity_norm();
if (verbose>1)
std::cout << "Largest Local Error: " << maxeta << " " << std::endl;
t_sqrt = timer3.elapsed();
timer3.reset();
// Stop adaptivity if target error is reached
if ((maxeta < alpha) && strategy == "targetTolerance") {
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, verbose-1);
else if (strategy == "errorFraction")
Dune::PDELab::error_fraction(eta, alpha, beta, eta_alpha, eta_beta, verbose-1);
else{ //((strategy == "threshold") || (strategy == "targetTolerance")) {
eta_alpha = alpha;
eta_beta = beta; }
t_strtgy = timer3.elapsed();
timer3.reset();
// Skip refinement if threshold is met, but still allow coarsening
if (maxeta < adaptivityThreshold) {
eta_alpha = maxeta + 1; // Only refine elements with error > maxerror + 1
if (verbose>1)
std::cout << " Max local error smaller than threshold " << adaptivityThreshold << ". Skipping grid refinement." << std::endl;
}
Dune::PDELab::mark_grid(grid, eta, eta_alpha, eta_beta, minLevel, maxLevel, verbose-1);
t_mark = timer3.elapsed();
timer3.reset();
Dune::PDELab::adapt_grid(grid, gfs, uold, unew, 2*order);
t_adapt = timer3.elapsed();
timer3.reset();
multiRefinementCounter++;
if(this->verbose>3){
ios_base_all_saver restore(std::cout);
std::cout << std::setprecision(4) << std::scientific;
std::cout << "::: setup time" << std::setw(12) << t_setup << std::endl;
std::cout << "::: error estimation time" << std::setw(12) << t_est << std::endl;
std::cout << "::: unsquaring time" << std::setw(12) << t_sqrt << std::endl;
std::cout << "::: strategy application time" << std::setw(12) << t_strtgy << std::endl;
std::cout << "::: grid marking time" << std::setw(12) << t_mark << std::endl;
std::cout << "::: grid adaptation time" << std::setw(12) << t_adapt << std::endl;
}
} while (strategy == "targetTolerance");
if(this->verbose>3)
std::cout << "::: grid adaptation time" << std::setw(12) << t_adapt << std::endl;
if (verbose>1) {
t2 = timer2.elapsed();
total_time += timer.elapsed();
ios_base_all_saver restore(std::cout);
std::cout << std::setprecision(4) << std::scientific;
std::cout << "::: total adaptation time" << std::setw(12) << t2 << std::endl;
std::cout << "::: total adaptation time" << std::setw(12) << total_time << std::endl;
}
return true;
total_time = 0.;
}
};
......
......@@ -165,13 +165,16 @@ void RichardsSimulation<Traits>::run ()
}
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
}
if( controller->doStep() )
{
// need "old" boundary condition
adaptivity->mark_grid(*grid, gv, *gfs, *param, *fboundary, time+dt,
*uold, *unew);
adaptivity->adapt_grid(*grid, *gfs, *uold, *unew);
operator_setup();
}
}
}
......
......@@ -278,13 +278,15 @@ public:
{
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();
}
this->adaptivity->mark_grid(*this->grid, this->gv, *this->gfs,
*this->param, *this->fboundary, time,
*this->uold, *this->unew);
this->adaptivity->adapt_grid( *this->grid, *this->gfs, *this->uold,
*this->unew);
this->operator_setup();
// compute new water content
const auto wc_new = water_content(*this->unew);
......
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