4 #ifndef DUNE_PDELAB_ADAPTIVITY_ADAPTIVITY_HH
5 #define DUNE_PDELAB_ADAPTIVITY_ADAPTIVITY_HH
7 #include<dune/common/exceptions.hh>
13 #include<unordered_map>
14 #include<dune/common/dynmatrix.hh>
15 #include<dune/geometry/quadraturerules.hh>
26 #include<dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
32 template<
typename GFS>
36 typedef typename GFS::Traits::GridView::template Codim<0>::Entity
Cell;
41 typedef std::array<std::size_t,TypeTree::TreeInfo<GFS>::leafCount + 1>
LeafOffsets;
47 assert(leaf_offsets.back() > 0);
54 if (leaf_offsets.back() == 0)
59 std::partial_sum(leaf_offsets.begin(),leaf_offsets.end(),leaf_offsets.begin());
61 assert(leaf_offsets.back() ==
_lfs.size());
78 template<
typename MassMatrices,
typename Cell>
79 struct inverse_mass_matrix_calculator
80 :
public TypeTree::TreeVisitor
81 ,
public TypeTree::DynamicTraversal
84 static const int dim = Cell::Geometry::mydimension;
85 typedef std::size_t size_type;
86 typedef typename MassMatrices::value_type MassMatrix;
87 typedef typename MassMatrix::field_type DF;
88 typedef typename Dune::QuadratureRule<DF,dim>::const_iterator QRIterator;
90 template<
typename GFS,
typename TreePath>
91 void leaf(
const GFS& gfs, TreePath treePath)
93 auto& fem = gfs.finiteElementMap();
95 size_type local_size = fe.localBasis().size();
98 mass_matrix.resize(local_size,local_size);
100 using Range =
typename GFS::Traits::FiniteElementMap::Traits::
101 FiniteElement::Traits::LocalBasisType::Traits::RangeType;
102 std::vector<Range> phi;
103 phi.resize(std::max(phi.size(),local_size));
107 fe.localBasis().evaluateFunction(ip.position(),phi);
108 const DF factor = ip.weight();
110 for (size_type i = 0; i < local_size; ++i)
111 for (size_type j = 0; j < local_size; ++j)
112 mass_matrix[i][j] += phi[i] * phi[j] * factor;
115 mass_matrix.invert();
120 inverse_mass_matrix_calculator(MassMatrices& mass_matrices,
const Cell& element, size_type intorder)
144 template<
class GFS,
class U>
147 using EntitySet =
typename GFS::Traits::EntitySet;
148 using Element =
typename EntitySet::Element;
150 typedef typename U::ElementType DF;
155 typedef std::array<MassMatrix,TypeTree::TreeInfo<GFS>::leafCount>
MassMatrices;
163 , _intorder(intorder)
164 , _inverse_mass_matrices(GlobalGeometryTypeIndex::size(Element::dimension))
174 auto gt =
e.geometry().type();
177 if (inverse_mass_matrices[0].N() > 0)
178 return inverse_mass_matrices;
180 inverse_mass_matrix_calculator<MassMatrices,Element> calculate_mass_matrices(
181 inverse_mass_matrices,
186 TypeTree::applyToTree(_gfs,calculate_mass_matrices);
188 return inverse_mass_matrices;
195 std::vector<MassMatrices> _inverse_mass_matrices;
199 template<
typename GFS,
typename DOFVector,
typename TransferMap>
201 :
public TypeTree::TreeVisitor
202 ,
public TypeTree::DynamicTraversal
210 using IDSet =
typename EntitySet::Traits::GridView::Grid::LocalIdSet;
213 static const int dim = Geometry::mydimension;
214 typedef typename DOFVector::ElementType
RF;
223 using DF =
typename EntitySet::Traits::CoordinateField;
225 template<
typename LFSLeaf,
typename TreePath>
226 void leaf(
const LFSLeaf& leaf_lfs, TreePath treePath)
229 auto& fem = leaf_lfs.gridFunctionSpace().finiteElementMap();
233 using Range =
typename LFSLeaf::Traits::GridFunctionSpace::Traits::FiniteElementMap::
234 Traits::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
238 auto coarse_phi = std::vector<Range>{};
239 auto fine_phi = std::vector<Range>{};
241 auto fine_geometry =
_current.geometry();
242 auto coarse_geometry =
_ancestor.geometry();
247 auto coarse_local = coarse_geometry.local(fine_geometry.global(ip.position()));
249 fe->localBasis().evaluateFunction(ip.position(),fine_phi);
251 fe->localBasis().evaluateFunction(coarse_local,coarse_phi);
252 const DF factor = ip.weight()
253 * fine_geometry.integrationElement(ip.position())
254 / coarse_geometry.integrationElement(coarse_local);
256 auto val = Range{0.0};
257 for (
size_type i = 0; i < fine_phi.size(); ++i)
259 val.axpy(
_u_fine[fine_offset + i],fine_phi[i]);
262 assert(inverse_mass_matrix.M()==coarse_phi.size());
263 for (
size_type i = 0; i < coarse_phi.size(); ++i)
266 for (
size_type j = 0; j < inverse_mass_matrix.M(); ++j)
267 x.axpy(inverse_mass_matrix[i][j],coarse_phi[j]);
268 (*_u_coarse)[coarse_offset + i] += factor * (x * val);
289 size_type max_level =
_lfs.gridFunctionSpace().gridView().grid().maxLevel();
309 for (
const auto& child : descendantElements(
_ancestor,max_level))
326 TypeTree::applyToTree(
_lfs,*
this);
336 TransferMap& transfer_map,
337 std::size_t int_order = 2)
340 ,
_id_set(gfs.gridView().grid().localIdSet())
357 typename DOFVector::template ConstLocalView<LFSCache>
_u_view;
369 template<
typename GFS,
typename DOFVector,
typename CountVector>
371 :
public TypeTree::TreeVisitor
372 ,
public TypeTree::DynamicTraversal
380 using IDSet =
typename EntitySet::Traits::GridView::Grid::LocalIdSet;
383 typedef typename DOFVector::ElementType
RF;
388 using DF =
typename EntitySet::Traits::CoordinateField;
390 template<
typename FiniteElement>
393 using Range =
typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
394 using Domain =
typename FiniteElement::Traits::LocalBasisType::Traits::DomainType;
424 mutable std::vector<Range>
_phi;
430 template<
typename LeafLFS,
typename TreePath>
431 void leaf(
const LeafLFS& leaf_lfs, TreePath treePath)
433 using FiniteElement =
typename LeafLFS::Traits::FiniteElementType;
435 auto& fem = leaf_lfs.gridFunctionSpace().finiteElementMap();
442 _u_tmp.resize(fe.localBasis().size());
444 fe.localInterpolation().interpolate(f,
_u_tmp);
471 TypeTree::applyToTree(
_lfs,*
this);
485 ,
_id_set(gfs.entitySet().gridView().grid().localIdSet())
497 typename DOFVector::template LocalView<LFSCache>
_u_view;
498 typename CountVector::template LocalView<LFSCache>
_uc_view;
522 template<
class Gr
id,
class GFSU,
class U,
class Projection>
525 typedef typename Grid::LeafGridView LeafGridView;
526 typedef typename LeafGridView::template Codim<0>
527 ::template Partition<Dune::Interior_Partition>::Iterator LeafIterator;
528 typedef typename Grid::template Codim<0>::Entity Element;
529 typedef typename Grid::LocalIdSet IDSet;
530 typedef typename IDSet::IdType ID;
533 typedef std::unordered_map<ID,std::vector<typename U::ElementType> >
MapType;
541 : _leaf_offset_cache(gfs)
553 Visitor visitor(gfsu,projection,u,_leaf_offset_cache,transfer_map);
556 for(
const auto& cell : elements(gfsu.entitySet(),Partitions::interior))
565 void replayData(Grid& grid, GFSU& gfsu, Projection& projection, U& u,
const MapType& transfer_map)
567 const IDSet& id_set = grid.localIdSet();
570 CountVector uc(gfsu,0);
573 Visitor visitor(gfsu,u,uc,_leaf_offset_cache);
576 for (
const auto& cell : elements(gfsu.entitySet(),Partitions::interior))
578 Element ancestor = cell;
580 typename MapType::const_iterator map_it;
581 while ((map_it = transfer_map.find(id_set.id(ancestor))) == transfer_map.end())
583 if (!ancestor.hasFather())
585 "transferMap of GridAdaptor didn't contain ancestor of element with id " << id_set.id(ancestor));
586 ancestor = ancestor.father();
589 visitor(cell,ancestor,map_it->second);
593 DOFHandle addHandle1(gfsu,u);
594 gfsu.entitySet().gridView().communicate(addHandle1,
595 Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
597 CountHandle addHandle2(gfsu,uc);
598 gfsu.entitySet().gridView().communicate(addHandle2,
599 Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
602 typename CountVector::iterator ucit = uc.begin();
603 for (
typename U::iterator uit = u.begin(), uend = u.end(); uit != uend; ++uit, ++ucit)
604 (*uit) /= ((*ucit) > 0 ? (*ucit) : 1.0);
631 template<
class Gr
id,
class GFS,
class X>
635 Projection projection(gfs,int_order);
644 grid_adaptor.
backupData(grid,gfs,projection,x1,transferMap1);
654 grid_adaptor.
replayData(grid,gfs,projection,x1,transferMap1);
671 template<
class Gr
id,
class GFS,
class X>
672 void adapt_grid (Grid& grid, GFS& gfs, X& x1, X& x2,
int int_order)
675 Projection projection(gfs,int_order);
684 grid_adaptor.
backupData(grid,gfs,projection,x1,transferMap1);
686 grid_adaptor.
backupData(grid,gfs,projection,x2,transferMap2);
696 grid_adaptor.
replayData(grid,gfs,projection,x1,transferMap1);
698 grid_adaptor.
replayData(grid,gfs,projection,x2,transferMap2);
708 template <
typename G,
typename... X>
709 struct GFSWithVectors
713 using Tuple = std::tuple<X&...>;
715 GFSWithVectors (GFS& gfs,
int integrationOrder, X&... x) :
717 _integrationOrder(integrationOrder),
722 int _integrationOrder;
727 template <
typename Gr
id>
728 void iteratePacks(Grid& grid);
729 template <
typename Grid,
typename X,
typename... XS>
730 void iteratePacks(Grid& grid, X& x, XS&... xs);
735 template<std::size_t I = 0,
typename Grid,
typename X,
typename... XS>
737 iterateTuple(Grid& grid, X& x, XS&... xs)
740 iteratePacks(grid,xs...);
755 template<std::size_t I = 0,
typename Grid,
typename X,
typename... XS>
757 iterateTuple(Grid& grid, X& x, XS&... xs)
760 using GFS =
typename X::GFS;
761 using Tuple =
typename X::Tuple;
762 using V =
typename std::decay<typename std::tuple_element<I,Tuple>::type>::type;
769 Projection projection(x._gfs,x._integrationOrder);
770 GridAdaptor<Grid,GFS,V,Projection> gridAdaptor(x._gfs);
774 gridAdaptor.backupData(grid,x._gfs,projection,std::get<I>(x._tuple),transferMap);
778 iterateTuple<I + 1, Grid, X, XS...>(grid,x,xs...);
782 std::get<I>(x._tuple) = V(x._gfs,0.0);
783 gridAdaptor.replayData(grid,x._gfs,projection,std::get<I>(x._tuple),transferMap);
789 template <
typename Gr
id>
790 void iteratePacks(Grid& grid)
811 template <
typename Grid,
typename X,
typename... XS>
812 void iteratePacks(Grid& grid, X& x, XS&... xs)
814 iterateTuple(grid,x,xs...);
831 template <
typename GFS,
typename... X>
834 impl::GFSWithVectors<GFS,X...> gfsWithVectors(gfs, integrationOrder, x...);
835 return gfsWithVectors;
848 template <
typename Grid,
typename... X>
855 impl::iteratePacks(grid,x...);
863 void error_fraction(
const T& x,
typename T::ElementType alpha,
typename T::ElementType beta,
864 typename T::ElementType& eta_alpha,
typename T::ElementType& eta_beta,
int verbose=0)
867 std::cout <<
"+++ error fraction: alpha=" << alpha <<
" beta=" << beta << std::endl;
869 typedef typename T::ElementType NumberType;
870 NumberType total_error = x.one_norm();
871 NumberType max_error = x.infinity_norm();
872 NumberType eta_alpha_left = 0.0;
873 NumberType eta_alpha_right = max_error;
874 NumberType eta_beta_left = 0.0;
875 NumberType eta_beta_right = max_error;
876 for (
int j=1; j<=steps; j++)
878 eta_alpha = 0.5*(eta_alpha_left+eta_alpha_right);
879 eta_beta = 0.5*(eta_beta_left+eta_beta_right);
880 NumberType sum_alpha=0.0;
881 NumberType sum_beta=0.0;
882 unsigned int alpha_count = 0;
883 unsigned int beta_count = 0;
884 for (
const auto& error : x)
886 if (error >=eta_alpha) { sum_alpha += error; alpha_count++;}
887 if (error < eta_beta) { sum_beta += error; beta_count++;}
891 std::cout <<
"+++ " << j <<
" eta_alpha=" << eta_alpha <<
" alpha_fraction=" << sum_alpha/total_error
892 <<
" elements: " << alpha_count <<
" of " << x.N() << std::endl;
893 std::cout <<
"+++ " << j <<
" eta_beta=" << eta_beta <<
" beta_fraction=" << sum_beta/total_error
894 <<
" elements: " << beta_count <<
" of " << x.N() << std::endl;
896 if (std::abs(alpha-sum_alpha/total_error) <= 0.01 && std::abs(beta-sum_beta/total_error) <= 0.01)
break;
897 if (sum_alpha>alpha*total_error)
898 eta_alpha_left = eta_alpha;
900 eta_alpha_right = eta_alpha;
901 if (sum_beta>beta*total_error)
902 eta_beta_right = eta_beta;
904 eta_beta_left = eta_beta;
908 std::cout <<
"+++ refine_threshold=" << eta_alpha
909 <<
" coarsen_threshold=" << eta_beta << std::endl;
915 void element_fraction(
const T& x,
typename T::ElementType alpha,
typename T::ElementType beta,
916 typename T::ElementType& eta_alpha,
typename T::ElementType& eta_beta,
int verbose=0)
919 typedef typename T::ElementType NumberType;
920 NumberType total_error =x.N();
921 NumberType max_error = x.infinity_norm();
922 NumberType eta_alpha_left = 0.0;
923 NumberType eta_alpha_right = max_error;
924 NumberType eta_beta_left = 0.0;
925 NumberType eta_beta_right = max_error;
926 for (
int j=1; j<=steps; j++)
928 eta_alpha = 0.5*(eta_alpha_left+eta_alpha_right);
929 eta_beta = 0.5*(eta_beta_left+eta_beta_right);
930 NumberType sum_alpha=0.0;
931 NumberType sum_beta=0.0;
932 unsigned int alpha_count = 0;
933 unsigned int beta_count = 0;
935 for (
const auto& error : x)
937 if (error>=eta_alpha) { sum_alpha += 1.0; alpha_count++;}
938 if (error< eta_beta) { sum_beta +=1.0; beta_count++;}
942 std::cout << j <<
" eta_alpha=" << eta_alpha <<
" alpha_fraction=" << sum_alpha/total_error
943 <<
" elements: " << alpha_count <<
" of " << x.N() << std::endl;
944 std::cout << j <<
" eta_beta=" << eta_beta <<
" beta_fraction=" << sum_beta/total_error
945 <<
" elements: " << beta_count <<
" of " << x.N() << std::endl;
947 if (std::abs(alpha-sum_alpha/total_error) <= 0.01 && std::abs(beta-sum_beta/total_error) <= 0.01)
break;
948 if (sum_alpha>alpha*total_error)
949 eta_alpha_left = eta_alpha;
951 eta_alpha_right = eta_alpha;
952 if (sum_beta>beta*total_error)
953 eta_beta_right = eta_beta;
955 eta_beta_left = eta_beta;
959 std::cout <<
"+++ refine_threshold=" << eta_alpha
960 <<
" coarsen_threshold=" << eta_beta << std::endl;
970 typedef typename T::ElementType NumberType;
971 NumberType total_error = x.one_norm();
972 NumberType total_elements = x.N();
973 NumberType max_error = x.infinity_norm();
974 std::vector<NumberType> left(bins,0.0);
975 std::vector<NumberType> right(bins,max_error*(1.0+1
e-8));
976 std::vector<NumberType> eta(bins);
977 std::vector<NumberType> target(bins);
978 for (
unsigned int k=0; k<bins; k++)
979 target[k]= (k+1)/((NumberType)bins);
980 for (
int j=1; j<=steps; j++)
982 for (
unsigned int k=0; k<bins; k++)
983 eta[k]= 0.5*(left[k]+right[k]);
984 std::vector<NumberType> sum(bins,0.0);
985 std::vector<int> count(bins,0);
987 for (
typename T::const_iterator it = x.begin(),
992 for (
unsigned int k=0; k<bins; k++)
1004 for (
unsigned int k=0; k<bins; k++)
1005 if (sum[k]<=target[k]*total_error)
1010 std::vector<NumberType> sum(bins,0.0);
1011 std::vector<int> count(bins,0);
1012 for (
unsigned int i=0; i<x.N(); i++)
1013 for (
unsigned int k=0; k<bins; k++)
1019 std::cout <<
"+++ error distribution" << std::endl;
1020 std::cout <<
"+++ number of elements: " << x.N() << std::endl;
1021 std::cout <<
"+++ max element error: " << max_error << std::endl;
1022 std::cout <<
"+++ total error: " << total_error << std::endl;
1023 std::cout <<
"+++ bin #elements eta sum/total " << std::endl;
1024 for (
unsigned int k=0; k<bins; k++)
1025 std::cout <<
"+++ " << k+1 <<
" " << count[k] <<
" " << eta[k] <<
" " << sum[k]/total_error << std::endl;
1028 template<
typename Gr
id,
typename X>
1029 void mark_grid (Grid &grid,
const X& x,
typename X::ElementType refine_threshold,
1030 typename X::ElementType coarsen_threshold,
int min_level = 0,
int max_level = std::numeric_limits<int>::max(),
int verbose=0)
1032 typedef typename Grid::LeafGridView GV;
1034 GV gv = grid.leafGridView();
1036 unsigned int refine_cnt=0;
1037 unsigned int coarsen_cnt=0;
1039 typedef typename X::GridFunctionSpace GFS;
1042 typedef typename X::template ConstLocalView<LFSCache> XView;
1044 LFS lfs(x.gridFunctionSpace());
1045 LFSCache lfs_cache(lfs);
1048 for(
const auto& cell : elements(gv))
1052 x_view.bind(lfs_cache);
1054 if (x_view[0]>=refine_threshold && cell.level() < max_level)
1059 if (x_view[0]<=coarsen_threshold && cell.level() > min_level)
1067 std::cout <<
"+++ mark_grid: " << refine_cnt <<
" marked for refinement, "
1068 << coarsen_cnt <<
" marked for coarsening" << std::endl;
1072 template<
typename Gr
id,
typename X>
1074 typename X::ElementType coarsen_threshold,
int verbose=0)
1076 typedef typename Grid::LeafGridView GV;
1078 GV gv = grid.leafGridView();
1080 unsigned int coarsen_cnt=0;
1082 typedef typename X::GridFunctionSpace GFS;
1085 typedef typename X::template ConstLocalView<LFSCache> XView;
1087 LFS lfs(x.gridFunctionSpace());
1088 LFSCache lfs_cache(lfs);
1091 for(
const auto& cell : elements(gv))
1095 x_view.bind(lfs_cache);
1097 if (x_view[0]>=refine_threshold)
1102 if (x_view[0]<=coarsen_threshold)
1109 std::cout <<
"+++ mark_grid_for_coarsening: "
1110 << coarsen_cnt <<
" marked for coarsening" << std::endl;
1118 double optimistic_factor;
1119 double coarsen_limit;
1120 double balance_limit;
1125 double refine_fraction_while_refinement;
1126 double coarsen_fraction_while_refinement;
1127 double coarsen_fraction_while_coarsening;
1128 double timestep_decrease_factor;
1129 double timestep_increase_factor;
1139 bool have_decreased_time_step;
1140 bool have_refined_grid;
1143 double accumulated_estimated_error_squared;
1144 double minenergy_rate;
1148 : scaling(16.0), optimistic_factor(1.0), coarsen_limit(0.5), balance_limit(0.33333),
1149 tol(tol_), T(T_), verbose(verbose_), no_adapt(false),
1150 refine_fraction_while_refinement(0.7),
1151 coarsen_fraction_while_refinement(0.2),
1152 coarsen_fraction_while_coarsening(0.2),
1153 timestep_decrease_factor(0.5), timestep_increase_factor(1.5),
1154 accept(false), adapt_dt(false),
adapt_grid(false), newdt(1.0),
1155 have_decreased_time_step(false), have_refined_grid(false),
1156 accumulated_estimated_error_squared(0.0),
1163 timestep_decrease_factor=
s;
1168 timestep_increase_factor=
s;
1173 refine_fraction_while_refinement=
s;
1178 coarsen_fraction_while_refinement=
s;
1183 coarsen_fraction_while_coarsening=
s;
1208 optimistic_factor=
s;
1258 return accumulated_estimated_error_squared;
1264 have_decreased_time_step =
false;
1265 have_refined_grid =
false;
1268 template<
typename GM,
typename X>
1269 void evaluate_estimators (GM& grid,
double time,
double dt,
const X& eta_space,
const X& eta_time,
double energy_timeslab)
1276 double spatial_error = eta_space.one_norm();
1277 double temporal_error = scaling*eta_time.one_norm();
1278 double sum = spatial_error + temporal_error;
1280 double allowed = tol*tol*(energy_timeslab+minenergy_rate*dt);
1281 q_s = spatial_error/sum;
1282 q_t = temporal_error/sum;
1289 <<
" sum/allowed=" << sum/allowed
1291 <<
" estimated error=" << sqrt(accumulated_estimated_error_squared+sum)
1292 <<
" energy_rate=" << energy_timeslab/dt
1299 accumulated_estimated_error_squared += sum;
1300 if (verbose>1) std::cout <<
"+++ no adapt mode" << std::endl;
1309 if (verbose>1) std::cout <<
"+++ accepting time step" << std::endl;
1310 accumulated_estimated_error_squared += sum;
1313 if (sum<coarsen_limit*allowed)
1316 if (q_t<balance_limit)
1319 newdt = timestep_increase_factor*dt;
1321 if (verbose>1) std::cout <<
"+++ spatial error dominates: increase time step" << std::endl;
1325 if (q_s>balance_limit)
1328 newdt = timestep_increase_factor*dt;
1330 if (verbose>1) std::cout <<
"+++ increasing time step" << std::endl;
1333 double eta_refine, eta_coarsen;
1334 if (verbose>1) std::cout <<
"+++ mark grid for coarsening" << std::endl;
1337 coarsen_fraction_while_coarsening,eta_refine,eta_coarsen);
1345 if (q_t<balance_limit)
1348 newdt = timestep_increase_factor*dt;
1350 if (verbose>1) std::cout <<
"+++ spatial error dominates: increase time step" << std::endl;
1357 if (verbose>1) std::cout <<
"+++ will redo time step" << std::endl;
1358 if (q_t>1-balance_limit)
1361 newdt = timestep_decrease_factor*dt;
1363 have_decreased_time_step =
true;
1364 if (verbose>1) std::cout <<
"+++ decreasing time step only" << std::endl;
1368 if (q_t<balance_limit)
1370 if (!have_decreased_time_step)
1373 newdt = timestep_increase_factor*dt;
1375 if (verbose>1) std::cout <<
"+++ increasing time step" << std::endl;
1381 newdt = timestep_decrease_factor*dt;
1383 have_decreased_time_step =
true;
1384 if (verbose>1) std::cout <<
"+++ decreasing time step" << std::endl;
1387 double eta_refine, eta_coarsen;
1388 if (verbose>1) std::cout <<
"+++ BINGO mark grid for refinement and coarsening" << std::endl;
1391 coarsen_fraction_while_refinement,eta_refine,eta_coarsen,0);
1404 #endif // DUNE_PDELAB_ADAPTIVITY_ADAPTIVITY_HH