Commit 56d5f506 authored by Lukas Riedel's avatar Lukas Riedel

Merge branch...

Merge branch '90-split-adaptivityhandler-adapt_grid-method-into-two-parts-mark_grid-and-adapt_grid' into 'master'

Resolve "Split AdaptivityHandler::adapt_grid() method into two parts: mark_grid and adapt_grid"

Closes #90

See merge request !91
parents fedee2f2 6b88bd07
......@@ -27,6 +27,9 @@
* `RichardsSimulation` now has its own `RichardsSimulationTraits` derived from
`BaseTraits`, which defines all its member types. `BaseTraits` now have
reduced content and are intended to be shared between models/simulations.
* Grid adaptation now is done in two steps: (1) mark entities of the grid to
refine/coarse and (2) adapt the grid and project the degrees of freedom on
the new grid.
* Data structures for storing and accessing parameters. The new class
`FlowParameters` maps parameter sets and scaling factors to every grid cell
on the coarsest grid level, and only stores one set of parameters for
......
......@@ -380,14 +380,10 @@ adding an empty line, make text **bold** or ``monospaced``.
Coarsening occurs for all elements that carry an error smaller
than β.
**targetTolerance**: All elements with a local error η > α are being refined as often
as necessary until the requirement η < α is met for all grid cells.
Coarsening occurs for all elements that carry an error smaller
than β.
</definition>
<values> elementFraction, errorFraction, targetTolerance, threshold </values>
<values> elementFraction, errorFraction, threshold </values>
<suggestion> elementFraction </suggestion>
<comment> elementFraction, errorFraction, targetTolerance, threshold </comment>
<comment> elementFraction, errorFraction, threshold </comment>
</parameter>
<parameter name="refinementFraction">
......
......@@ -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& u) { }
/// 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& u) { return false; }
virtual void adapt_grid (Grid& grid, GFS& gfs, U& u) { }
};
/// 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){
......@@ -107,16 +112,18 @@ public:
}
}
/// Adapt the grid according to the estimated flux error and the strategy applied
/**
* \param grid Grid to adapt
* \param gv Leaf grid view of the grid
/// Mark grid cells for adaptation.
/** Marking is based on the estimated flux error and the strategy applied.
* \param grid Grid to mark (and later adapt)
* \param gv Current grid view
* \param gfs Solution GridFunctionSpace
* \param param Parameter handler
* \param param Parameter interface
* \param boundary Boundary condition interface
* \param time The time for which boundary conditions are queried inside
* the error estimator operator.
* \param u 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,
......@@ -125,116 +132,118 @@ public:
const RF time,
U& u) override
{
Dune::Timer timer2;
float t2;
double maxeta(0.0);
int multiRefinementCounter = 0;
Dune::Timer timer_total;
Dune::Timer timer_seq;
float t_setup, t_sqrt, t_est, t_strtgy, t_mark;
if(verbose>1){
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;
double eta_alpha(0.0);
double eta_beta(0.0);
double maxeta(0.0);
if(verbose>1)
std::cout << " Refinement Step " << multiRefinementCounter << ": ";
// 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 = timer_seq.elapsed();
timer_seq.reset();
// Compute error estimate eta
estgo.residual(u,eta);
t_est = timer_seq.elapsed();
timer_seq.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 = timer_seq.elapsed();
timer_seq.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 = timer_seq.elapsed();
timer_seq.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;
}
double eta_alpha(0.0);
double eta_beta(0.0);
Dune::PDELab::mark_grid(grid, eta, eta_alpha, eta_beta, minLevel, maxLevel, verbose-1);
// 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_mark = timer_seq.elapsed();
timer_seq.reset();
t_setup = timer3.elapsed();
timer3.reset();
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;
}
// Compute error estimate eta
estgo.residual(u,eta);
total_time += timer_total.elapsed();
}
t_est = timer3.elapsed();
timer3.reset();
/// Adapt the grid, and update grid function space and solution.
/** This assumes that the grid has been marked for refinement. Otherwise,
* calling this function has no effect.
* After adapting the grid, the grid function space is updated and the
* solution vector is interpolated onto the new grid.
* \param grid The (marked) grid to adapt
* \param gfs The solution grid function space
* \param u The solution vector to be interpolated
*/
void adapt_grid (
Grid& grid,
GFS& gfs,
U& u) override
{
Dune::Timer timer;
// 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); }
);
Dune::PDELab::adapt_grid(grid, gfs, u, 2*order);
// 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, u, 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");
float t_adapt = timer.elapsed();
total_time += t_adapt;
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.;
}
};
......
......@@ -184,13 +184,16 @@ void RichardsSimulation<Traits>::run ()
if constexpr (Traits::write_output)
write_data(controller->getTime());
if(controller->doStep()
&& adaptivity->adapt_grid(*grid, gv, *gfs, *fparam, *fboundary,
time+dt, // need "old" boundary condition
*u))
{
operator_setup(); // reset operators if grid changes
}
if( controller->doStep() )
{
// need "old" boundary condition
adaptivity->mark_grid(*grid, gv, *gfs, *fparam, *fboundary, time+dt, *u);
adaptivity->adapt_grid(*grid, *gfs, *u);
// reset operators always
operator_setup();
}
}
}
......
......@@ -284,13 +284,15 @@ public:
{
continue;
}
if (this->adaptivity->adapt_grid(
*this->grid, this->gv, *this->gfs, *this->fparam, *this->fboundary,
time, *this->u))
{
// reset operators if grid changes
this->operator_setup();
}
this->adaptivity->mark_grid(*this->grid, this->gv, *this->gfs,
*this->fparam, *this->fboundary, time,
*this->u);
this->adaptivity->adapt_grid( *this->grid, *this->gfs, *this->u);
// reset operators always
this->operator_setup();
// compute new water content
const auto wc_new = water_content(*this->u);
......
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