Skip to content
Snippets Groups Projects
Commit f6557b76 authored by Carsten Gräser's avatar Carsten Gräser
Browse files

[examples] Use a hierarchical estimator in adaptive example

parent dbfcfc37
Branches
Tags
1 merge request!258Add constraints support for dune-functions basis
Hide whitespace changes
Inline Side-by-side
examples/poisson-adaptive.cc
+ 227
− 91
View file @ f6557b76
......@@ -4,6 +4,8 @@
#include <vector>
#include <cmath>
#include <ranges>
#include <tuple>
#include <dune/common/bitsetvector.hh>
#include <dune/common/version.hh>
......@@ -22,16 +24,13 @@
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/hierarchicallagrangebasis.hh>
#include <dune/functions/functionspacebases/boundarydofs.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#if DUNE_VERSION_GT(DUNE_VTK, 2, 9)
#include <dune/vtk/vtkwriter.hh>
#include <dune/vtk/datacollectors/lagrangedatacollector.hh>
#else
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#endif
#include <dune/fufem/logger.hh>
#include <dune/fufem/parallel/elementcoloring.hh>
......@@ -45,30 +44,126 @@
using namespace Dune;
template<class BilinearForm, class Residual, class ExtensionBasis, class Constraints, class IsConstraints, class Vector, class... Args>
auto hierarchicalErrorEstimator(BilinearForm b, Residual residual, const ExtensionBasis& extension, const Constraints& constraints, const IsConstraints& isConstrained, Vector& defect, Args... args)
{
using namespace Dune::Fufem::Forms;
auto d = trialFunction(extension);
auto v = testFunction(extension);
auto r = Vector();
auto B = Dune::BCRSMatrix<double>();
auto globalAssembler = GlobalAssembler(extension, extension, 0, std::forward<Args>(args)...);
globalAssembler.assembleFunctional(r, residual(v), constraints);
globalAssembler.assembleOperator(B, b(d,v), constraints);
auto eta = std::vector<double>(extension.size(), 0.0);
double etaSum = 0;
for(auto k : Dune::range(extension.size()))
{
if (not isConstrained[k])
defect[k] = r[k]/B[k][k];
eta[k] = defect[k]*defect[k]*B[k][k];
etaSum += eta[k];
}
constraints.interpolate(defect);
return std::make_tuple(std::sqrt(etaSum), eta);
}
// Guess a value of theta for Doerfler marking.
// Doerfler's marking strategy uses a parameter \theta from the interval (0,1)
// to control the locallity of refinement, where larger values lead to more local
// refinement. He proposes values from the interval [0.25, .5] to balance locallity
// per refinement step with the number of needed refinement steps. In practice
// these values all work well. However, if we are very close to matching the global
// tolerance, a more local refinement might be sufficient to match the tolerance
// in the next step with less unknowns. Given the current estimanted error, the target
// tolerance and the expected local error reduction this function computes a guess for
// the value theta that suffices to match the tolerance.
double guessTheta(double tol, double error, double expectedReduction, double theta_min = 0.25, double theta_max = 0.5)
{
// The value beta is the fraction of the squared local error estimators that
// we need to refine such that we match the tolerance if the local errors
// reduce as expected after refinement.
double beta = (1.-(tol/error)*(tol/error))/(1.-expectedReduction*expectedReduction);
// For err>=tol the value beta should be non-negative.
// But we want to avoid rounding issues if error almost matches tol.
beta = std::max(beta,0.0);
// The paper by Doerfer uses a parameter theta
// with beta=(1-theta)^2
double theta = 1. - std::sqrt(beta);
template<int dim>
auto createUniformCubeGrid()
// Project into desired range of values for theta
return std::max(std::min(theta, theta_max), theta_min);
}
// Mark grid using the strategy proposed by Doerfler.
template<class Grid, class ExtensionBasis, class LocalIndicators>
std::size_t doerflerMarking(Grid& grid, const ExtensionBasis& extensionBasis, const LocalIndicators& eta, double error, double theta)
{
using Grid = Dune::YaspGrid<dim>;
Dune::FieldVector<double,dim> l(1.0);
std::array<int,dim> elements = {{2, 2}};
return std::make_unique<Grid>(l, elements);
auto etaSorted = eta;
std::sort(etaSorted.begin(), etaSorted.end(), std::greater{});
double etaRefinedTarget = (1-theta)*(1-theta)*error*error;
double etaRefined = 0;
double localTol = 0;
double localThreshold = 0;
for(auto eta_k : etaSorted)
{
etaRefined += eta_k;
localTol = eta_k;
if (etaRefined >= etaRefinedTarget)
break;
}
// Cosmetic hack:
// If the error is symmetric we may still get non-symmetric refinement due
// to round-off errors. By decreasing the tolerance slightly (here by 1%),
// it is more likely that we refine contributions that are the same up to
// round-off errors.
localTol *= 0.99;
std::size_t marked = 0;
auto localView = extensionBasis.localView();
for(const auto& element : elements(extensionBasis.gridView()))
{
localView.bind(element);
for(auto i : Dune::range(localView.size()))
{
auto k = localView.index(i);
if (eta[k] >= localTol)
{
grid.mark(1, element);
++marked;
}
}
}
return marked;
}
#if HAVE_DUNE_UGGRID
auto createMixedGrid()
auto createGrid()
{
using Grid = Dune::UGGrid<2>;
Dune::GridFactory<Grid> factory;
for(unsigned int k : Dune::range(8))
factory.insertVertex({0.5*(k%3), 0.5*(k/3)});
factory.insertElement(Dune::GeometryTypes::cube(2), {0, 1, 3, 4});
factory.insertElement(Dune::GeometryTypes::cube(2), {1, 2, 4, 5});
factory.insertElement(Dune::GeometryTypes::simplex(2), {0, 1, 3});
factory.insertElement(Dune::GeometryTypes::simplex(2), {1, 4, 3});
factory.insertElement(Dune::GeometryTypes::simplex(2), {1, 2, 4});
factory.insertElement(Dune::GeometryTypes::simplex(2), {2, 5, 4});
factory.insertElement(Dune::GeometryTypes::simplex(2), {3, 4, 6});
factory.insertElement(Dune::GeometryTypes::simplex(2), {4, 7, 6});
return std::unique_ptr<Grid>{factory.createGrid()};
}
#endif
int main (int argc, char *argv[]) try
{
......@@ -93,17 +188,15 @@ int main (int argc, char *argv[]) try
log("MPI initialized");
constexpr auto order = 2;
constexpr auto order = 1;
constexpr auto extensionOrder = 2;
///////////////////////////////////
// Generate the grid
///////////////////////////////////
#if HAVE_DUNE_UGGRID
auto gridPtr = createMixedGrid();
#else
auto gridPtr = createUniformCubeGrid<2>();
#endif
auto gridPtr = createGrid();
auto& grid = *gridPtr;
using Grid = std::decay_t<decltype(grid)>;
......@@ -122,9 +215,13 @@ int main (int argc, char *argv[]) try
auto dirichletIndicatorFunction = [](auto x) { return true; };
auto rightHandSide = [] (const auto& x) { return 10;};
auto dirichletValues = [](const auto& x){ return 0; };
// auto rightHandSide = [] (const auto& x) { return 0;};
// auto dirichletValues = [](const auto& x){
// return (std::fabs(x[0])<1e-5)*(x[1]<.25)*x[1]*(.25-x[1])*20.;
// };
bool refined = true;
while(refined)
std::size_t refStep = 0;
while(true)
{
auto gridView = grid.leafGridView();
......@@ -246,102 +343,141 @@ int main (int argc, char *argv[]) try
// Solve linear system
// *********************************************
// Technicality: turn the matrix into a linear operator
MatrixAdapter<Matrix,Vector,Vector> op(stiffnessMatrix);
{
// Technicality: turn the matrix into a linear operator
MatrixAdapter<Matrix,Vector,Vector> op(stiffnessMatrix);
// FastAMG is not working for non-blocked matrices in 2.10
#if DUNE_VERSION_GTE(DUNE_ISTL, 2, 11)
auto preconditioner = makeAMGPreconditioner<Vector>(stiffnessMatrix);
#else
// Sequential incomplete LU decomposition as the preconditioner
auto preconditioner = SeqILDL<Matrix,Vector,Vector>(stiffnessMatrix,1.0);
#endif
// FastAMG is not working for non-blocked matrices in 2.10
#if DUNE_VERSION_GTE(DUNE_ISTL, 2, 11)
auto preconditioner = makeAMGPreconditioner<Vector>(stiffnessMatrix);
#else
// Sequential incomplete LU decomposition as the preconditioner
auto preconditioner = SeqILDL<Matrix,Vector,Vector>(stiffnessMatrix,1.0);
#endif
log("Preconditioner created");
log("Preconditioner created");
// Preconditioned conjugate-gradient solver
CGSolver<Vector> cg(op,
preconditioner, // preconditioner
1e-4, // desired residual reduction factor
1000, // maximum number of iterations
2); // verbosity of the solver
// Preconditioned conjugate-gradient solver
CGSolver<Vector> cg(op,
preconditioner, // preconditioner
1e-5, // desired residual reduction factor
1000, // maximum number of iterations
2); // verbosity of the solver
log("Solver created");
log("Solver created");
// Object storing some statistics about the solving process
InverseOperatorResult statistics;
// Object storing some statistics about the solving process
InverseOperatorResult statistics;
// Solve!
cg.apply(sol, rhs, statistics);
// Solve!
cg.apply(sol, rhs, statistics);
log("Linear system solved");
log("Linear system solved");
constraints.interpolate(sol);
log("Constrained solution interpolated");
constraints.interpolate(sol);
log("Constrained solution interpolated");
}
// *********************************************
// Evaluate refinement indicator and adapt grid
// *********************************************
refined = false;
auto tol = 1e-2;
{
using namespace ::Dune::Fufem::Forms;
// A very crude refinement indicator:
// Refine until the local H^1 norm ||grad(u)||_L^2(e) < tol
auto tol = 1e-2;
auto u = bindToCoefficients(trialFunction(basis), sol);
auto constBasis = Dune::Functions::LagrangeBasis<GridView,0>(gridView);
auto extensionBasis = Dune::Functions::HierarchicalLagrangeBasis<GridView,extensionOrder>(gridView);
auto extensionConstraints = Dune::Fufem::makeContinuityConstraints<BitVector>(extensionBasis);
log("Constraints for extension space computed");
auto indicator = std::vector<double>();
auto v = testFunction(constBasis);
auto globalAssembler = GlobalAssembler(constBasis, constBasis, 0);
globalAssembler.assembleFunctional(indicator, integrate(dot(grad(u),grad(u))*v));
log("Refinement indicator evaluated");
auto defect = Vector();
Dune::Functions::istlVectorBackend(defect).resize(extensionBasis);
Dune::Functions::istlVectorBackend(defect) = 0;
auto localView = constBasis.localView();
for(const auto& element : elements(basis.gridView()))
// Mark all constrained and lower order DOFs to be ignored
auto ignore = extensionConstraints.isConstrained();
{
localView.bind(element);
if (std::sqrt(indicator[localView.index(0)]) > tol)
auto localView = extensionBasis.localView();
for(const auto& element : elements(extensionBasis.gridView()))
{
refined = true;
grid.mark(1, element);
localView.bind(element);
const auto& localCoefficients = localView.tree().finiteElement().localCoefficients();
for(auto i : Dune::range(localCoefficients.size()))
if (localCoefficients.localKey(i).codim() == Grid::dimension)
ignore[localView.index(i)] = true;
}
}
}
log("Grid marked for refinement");
if (refined)
{
// Interpolate all non-ignored boundary DOFs
auto extensionIsBoundary = std::vector<bool>();
Dune::Fufem::markBoundaryPatchDofs(dirichletPatch, extensionBasis, extensionIsBoundary);
for(auto i: Dune::range(ignore.size()))
if (ignore[i])
extensionIsBoundary[i] = false;
interpolate(extensionBasis, defect, dirichletValues, extensionIsBoundary);
log("Boundary DOFs for extension space marked and interpolated");
// Mark all boundary DOFs for being ignored, too.
for(auto i: Dune::range(ignore.size()))
if (extensionIsBoundary[i])
ignore[i] = true;
// Define defect problem assemblers
auto f = Coefficient(Dune::Functions::makeGridViewFunction(rightHandSide, basis.gridView()));
auto a = [&](auto u, auto v) {
return integrate(dot(grad(u),grad(v)));
};
auto residual = [&](auto v) {
return integrate(f*v - dot(grad(u),grad(v)));
};
// Compute local error estimates
#if DUNE_VERSION_GT(DUNE_FUNCTIONS, 2, 9)
auto [error, eta] = hierarchicalErrorEstimator(a, residual, extensionBasis, extensionConstraints, ignore, defect, std::cref(gridViewPartition), threadCount);
#else
auto [error, eta] = hierarchicalErrorEstimator(a, residual, extensionBasis, extensionConstraints, ignore, defect);
#endif
log(Dune::formatString("Estimated total error : % 12.5e", error));
// *********************************************
// Write solution
// *********************************************
{
auto d = bindToCoefficients(trialFunction(extensionBasis), defect);
using Dune::Vtk::LagrangeDataCollector;
using Dune::Vtk::UnstructuredGridWriter;
using Dune::VTK::FieldInfo;
auto vtkWriter = UnstructuredGridWriter(LagrangeDataCollector(basis.gridView(), extensionOrder));
vtkWriter.addPointData(u, FieldInfo("sol", FieldInfo::Type::scalar, 1));
vtkWriter.addPointData(d, FieldInfo("defect", FieldInfo::Type::scalar, 1));
vtkWriter.write(Dune::formatString("poisson-adaptive-%03d", refStep));
log("Solution written to vtk file");
}
if (error <= tol)
break;
double expectedReduction = .5;
double theta = guessTheta(tol, error, expectedReduction);
log(Dune::formatString("Computed guess for refinement parameter theta : % 12.5e", theta));
doerflerMarking(grid, extensionBasis, eta, error, theta);
log("Grid marked for refinement");
grid.preAdapt();
grid.adapt();
grid.postAdapt();
log("Grid adapted");
}
else
{
using namespace ::Dune::Fufem::Forms;
auto u = bindToCoefficients(trialFunction(basis), sol);
#if DUNE_VERSION_GT(DUNE_VTK, 2, 9)
using Dune::Vtk::LagrangeDataCollector;
using Dune::Vtk::UnstructuredGridWriter;
auto vtkWriter = UnstructuredGridWriter(LagrangeDataCollector(basis.gridView(), order));
vtkWriter.addPointData(u, VTK::FieldInfo("sol", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.write("poisson-adaptive");
#else
// We need to use the SubsamplingVTKWriter, because dune-vtk's LagrangeDataCollector
// does not work with mixed grids in 2.9.
SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(2));
vtkWriter.addVertexData(u, VTK::FieldInfo("sol", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.write("poisson-adaptive");
#endif
log("Solution written to vtk file");
log(Dune::formatString("Grid adapted. New grid has %d levels.", grid.maxLevel()));
refStep++;
}
}
}
}
// Error handling
catch (Exception& e) {
std::cout << e.what() << std::endl;
}
catch (Exception& e) {
std::cout << e.what() << std::endl;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment