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 for (
size_type i = 0; i < coarse_phi.size(); ++i)
265 for (
size_type j = 0; j < inverse_mass_matrix.M(); ++j)
266 x.axpy(inverse_mass_matrix[i][j],coarse_phi[j]);
267 (*_u_coarse)[coarse_offset + i] += factor * (x * val);
288 size_type max_level =
_lfs.gridFunctionSpace().gridView().grid().maxLevel();
306 for (
const auto& child : descendantElements(
_ancestor,max_level))
323 TypeTree::applyToTree(
_lfs,*
this);
333 TransferMap& transfer_map,
334 std::size_t int_order = 2)
337 ,
_id_set(gfs.gridView().grid().localIdSet())
354 typename DOFVector::template ConstLocalView<LFSCache>
_u_view;
366 template<
typename GFS,
typename DOFVector,
typename CountVector>
368 :
public TypeTree::TreeVisitor
369 ,
public TypeTree::DynamicTraversal
377 using IDSet =
typename EntitySet::Traits::GridView::Grid::LocalIdSet;
380 typedef typename DOFVector::ElementType
RF;
385 using DF =
typename EntitySet::Traits::CoordinateField;
387 template<
typename FiniteElement>
390 using Range =
typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
392 template<
typename X,
typename Y>
414 mutable std::vector<Range>
_phi;
420 template<
typename LeafLFS,
typename TreePath>
421 void leaf(
const LeafLFS& leaf_lfs, TreePath treePath)
423 using FiniteElement =
typename LeafLFS::Traits::FiniteElementType;
425 auto& fem = leaf_lfs.gridFunctionSpace().finiteElementMap();
432 _u_tmp.resize(fe.localBasis().size());
434 fe.localInterpolation().interpolate(f,
_u_tmp);
461 TypeTree::applyToTree(
_lfs,*
this);
475 ,
_id_set(gfs.entitySet().gridView().grid().localIdSet())
487 typename DOFVector::template LocalView<LFSCache>
_u_view;
488 typename CountVector::template LocalView<LFSCache>
_uc_view;
512 template<
class Gr
id,
class GFSU,
class U,
class Projection>
515 typedef typename Grid::LeafGridView LeafGridView;
516 typedef typename LeafGridView::template Codim<0>
517 ::template Partition<Dune::Interior_Partition>::Iterator LeafIterator;
518 typedef typename Grid::template Codim<0>::Entity Element;
519 typedef typename Grid::LocalIdSet IDSet;
520 typedef typename IDSet::IdType ID;
523 typedef std::unordered_map<ID,std::vector<typename U::ElementType> >
MapType;
531 : _leaf_offset_cache(gfs)
543 Visitor visitor(gfsu,projection,u,_leaf_offset_cache,transfer_map);
546 for(
const auto& cell : elements(gfsu.entitySet(),Partitions::interior))
555 void replayData(Grid& grid, GFSU& gfsu, Projection& projection, U& u,
const MapType& transfer_map)
557 const IDSet& id_set = grid.localIdSet();
560 CountVector uc(gfsu,0);
563 Visitor visitor(gfsu,u,uc,_leaf_offset_cache);
566 for (
const auto& cell : elements(gfsu.entitySet(),Partitions::interior))
568 Element ancestor = cell;
570 typename MapType::const_iterator map_it;
571 while ((map_it = transfer_map.find(id_set.id(ancestor))) == transfer_map.end())
573 if (!ancestor.hasFather())
575 "transferMap of GridAdaptor didn't contain ancestor of element with id " << id_set.id(ancestor));
576 ancestor = ancestor.father();
579 visitor(cell,ancestor,map_it->second);
583 DOFHandle addHandle1(gfsu,u);
584 gfsu.entitySet().gridView().communicate(addHandle1,
585 Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
587 CountHandle addHandle2(gfsu,uc);
588 gfsu.entitySet().gridView().communicate(addHandle2,
589 Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
592 typename CountVector::iterator ucit = uc.begin();
593 for (
typename U::iterator uit = u.begin(), uend = u.end(); uit != uend; ++uit, ++ucit)
594 (*uit) /= ((*ucit) > 0 ? (*ucit) : 1.0);
621 template<
class Gr
id,
class GFS,
class X>
625 Projection projection(gfs,int_order);
634 grid_adaptor.
backupData(grid,gfs,projection,x1,transferMap1);
644 grid_adaptor.
replayData(grid,gfs,projection,x1,transferMap1);
661 template<
class Gr
id,
class GFS,
class X>
662 void adapt_grid (Grid& grid, GFS& gfs, X& x1, X& x2,
int int_order)
665 Projection projection(gfs,int_order);
674 grid_adaptor.
backupData(grid,gfs,projection,x1,transferMap1);
676 grid_adaptor.
backupData(grid,gfs,projection,x2,transferMap2);
686 grid_adaptor.
replayData(grid,gfs,projection,x1,transferMap1);
688 grid_adaptor.
replayData(grid,gfs,projection,x2,transferMap2);
698 template <
typename G,
typename... X>
699 struct GFSWithVectors
703 using Tuple = std::tuple<X&...>;
705 GFSWithVectors (GFS& gfs,
int integrationOrder, X&... x) :
707 _integrationOrder(integrationOrder),
712 int _integrationOrder;
717 template <
typename Gr
id>
718 void iteratePacks(Grid& grid);
719 template <
typename Grid,
typename X,
typename... XS>
720 void iteratePacks(Grid& grid, X& x, XS&... xs);
725 template<std::size_t I = 0,
typename Grid,
typename X,
typename... XS>
727 iterateTuple(Grid& grid, X& x, XS&... xs)
730 iteratePacks(grid,xs...);
745 template<std::size_t I = 0,
typename Grid,
typename X,
typename... XS>
747 iterateTuple(Grid& grid, X& x, XS&... xs)
750 using GFS =
typename X::GFS;
751 using Tuple =
typename X::Tuple;
752 using V =
typename std::decay<typename std::tuple_element<I,Tuple>::type>::type;
759 Projection projection(x._gfs,x._integrationOrder);
760 GridAdaptor<Grid,GFS,V,Projection> gridAdaptor(x._gfs);
764 gridAdaptor.backupData(grid,x._gfs,projection,std::get<I>(x._tuple),transferMap);
768 iterateTuple<I + 1, Grid, X, XS...>(grid,x,xs...);
772 std::get<I>(x._tuple) = V(x._gfs,0.0);
773 gridAdaptor.replayData(grid,x._gfs,projection,std::get<I>(x._tuple),transferMap);
779 template <
typename Gr
id>
780 void iteratePacks(Grid& grid)
801 template <
typename Grid,
typename X,
typename... XS>
802 void iteratePacks(Grid& grid, X& x, XS&... xs)
804 iterateTuple(grid,x,xs...);
821 template <
typename GFS,
typename... X>
824 impl::GFSWithVectors<GFS,X...> gfsWithVectors(gfs, integrationOrder, x...);
825 return gfsWithVectors;
838 template <
typename Grid,
typename... X>
845 impl::iteratePacks(grid,x...);
853 void error_fraction(
const T& x,
typename T::ElementType alpha,
typename T::ElementType beta,
854 typename T::ElementType& eta_alpha,
typename T::ElementType& eta_beta,
int verbose=0)
857 std::cout <<
"+++ error fraction: alpha=" << alpha <<
" beta=" << beta << std::endl;
859 typedef typename T::ElementType NumberType;
860 NumberType total_error = x.one_norm();
861 NumberType max_error = x.infinity_norm();
862 NumberType eta_alpha_left = 0.0;
863 NumberType eta_alpha_right = max_error;
864 NumberType eta_beta_left = 0.0;
865 NumberType eta_beta_right = max_error;
866 for (
int j=1; j<=steps; j++)
868 eta_alpha = 0.5*(eta_alpha_left+eta_alpha_right);
869 eta_beta = 0.5*(eta_beta_left+eta_beta_right);
870 NumberType sum_alpha=0.0;
871 NumberType sum_beta=0.0;
872 unsigned int alpha_count = 0;
873 unsigned int beta_count = 0;
874 for (
const auto& error : x)
876 if (error >=eta_alpha) { sum_alpha += error; alpha_count++;}
877 if (error < eta_beta) { sum_beta += error; beta_count++;}
881 std::cout <<
"+++ " << j <<
" eta_alpha=" << eta_alpha <<
" alpha_fraction=" << sum_alpha/total_error
882 <<
" elements: " << alpha_count <<
" of " << x.N() << std::endl;
883 std::cout <<
"+++ " << j <<
" eta_beta=" << eta_beta <<
" beta_fraction=" << sum_beta/total_error
884 <<
" elements: " << beta_count <<
" of " << x.N() << std::endl;
886 if (std::abs(alpha-sum_alpha/total_error) <= 0.01 && std::abs(beta-sum_beta/total_error) <= 0.01)
break;
887 if (sum_alpha>alpha*total_error)
888 eta_alpha_left = eta_alpha;
890 eta_alpha_right = eta_alpha;
891 if (sum_beta>beta*total_error)
892 eta_beta_right = eta_beta;
894 eta_beta_left = eta_beta;
898 std::cout <<
"+++ refine_threshold=" << eta_alpha
899 <<
" coarsen_threshold=" << eta_beta << std::endl;
905 void element_fraction(
const T& x,
typename T::ElementType alpha,
typename T::ElementType beta,
906 typename T::ElementType& eta_alpha,
typename T::ElementType& eta_beta,
int verbose=0)
909 typedef typename T::ElementType NumberType;
910 NumberType total_error =x.N();
911 NumberType max_error = x.infinity_norm();
912 NumberType eta_alpha_left = 0.0;
913 NumberType eta_alpha_right = max_error;
914 NumberType eta_beta_left = 0.0;
915 NumberType eta_beta_right = max_error;
916 for (
int j=1; j<=steps; j++)
918 eta_alpha = 0.5*(eta_alpha_left+eta_alpha_right);
919 eta_beta = 0.5*(eta_beta_left+eta_beta_right);
920 NumberType sum_alpha=0.0;
921 NumberType sum_beta=0.0;
922 unsigned int alpha_count = 0;
923 unsigned int beta_count = 0;
925 for (
const auto& error : x)
927 if (error>=eta_alpha) { sum_alpha += 1.0; alpha_count++;}
928 if (error< eta_beta) { sum_beta +=1.0; beta_count++;}
932 std::cout << j <<
" eta_alpha=" << eta_alpha <<
" alpha_fraction=" << sum_alpha/total_error
933 <<
" elements: " << alpha_count <<
" of " << x.N() << std::endl;
934 std::cout << j <<
" eta_beta=" << eta_beta <<
" beta_fraction=" << sum_beta/total_error
935 <<
" elements: " << beta_count <<
" of " << x.N() << std::endl;
937 if (std::abs(alpha-sum_alpha/total_error) <= 0.01 && std::abs(beta-sum_beta/total_error) <= 0.01)
break;
938 if (sum_alpha>alpha*total_error)
939 eta_alpha_left = eta_alpha;
941 eta_alpha_right = eta_alpha;
942 if (sum_beta>beta*total_error)
943 eta_beta_right = eta_beta;
945 eta_beta_left = eta_beta;
949 std::cout <<
"+++ refine_threshold=" << eta_alpha
950 <<
" coarsen_threshold=" << eta_beta << std::endl;
960 typedef typename T::ElementType NumberType;
961 NumberType total_error = x.one_norm();
962 NumberType total_elements = x.N();
963 NumberType max_error = x.infinity_norm();
964 std::vector<NumberType> left(bins,0.0);
965 std::vector<NumberType> right(bins,max_error*(1.0+1
e-8));
966 std::vector<NumberType> eta(bins);
967 std::vector<NumberType> target(bins);
968 for (
unsigned int k=0; k<bins; k++)
969 target[k]= (k+1)/((NumberType)bins);
970 for (
int j=1; j<=steps; j++)
972 for (
unsigned int k=0; k<bins; k++)
973 eta[k]= 0.5*(left[k]+right[k]);
974 std::vector<NumberType> sum(bins,0.0);
975 std::vector<int> count(bins,0);
977 for (
typename T::const_iterator it = x.begin(),
982 for (
unsigned int k=0; k<bins; k++)
994 for (
unsigned int k=0; k<bins; k++)
995 if (sum[k]<=target[k]*total_error)
1000 std::vector<NumberType> sum(bins,0.0);
1001 std::vector<int> count(bins,0);
1002 for (
unsigned int i=0; i<x.N(); i++)
1003 for (
unsigned int k=0; k<bins; k++)
1009 std::cout <<
"+++ error distribution" << std::endl;
1010 std::cout <<
"+++ number of elements: " << x.N() << std::endl;
1011 std::cout <<
"+++ max element error: " << max_error << std::endl;
1012 std::cout <<
"+++ total error: " << total_error << std::endl;
1013 std::cout <<
"+++ bin #elements eta sum/total " << std::endl;
1014 for (
unsigned int k=0; k<bins; k++)
1015 std::cout <<
"+++ " << k+1 <<
" " << count[k] <<
" " << eta[k] <<
" " << sum[k]/total_error << std::endl;
1018 template<
typename Gr
id,
typename X>
1019 void mark_grid (Grid &grid,
const X& x,
typename X::ElementType refine_threshold,
1020 typename X::ElementType coarsen_threshold,
int min_level = 0,
int max_level = std::numeric_limits<int>::max(),
int verbose=0)
1022 typedef typename Grid::LeafGridView GV;
1024 GV gv = grid.leafGridView();
1026 unsigned int refine_cnt=0;
1027 unsigned int coarsen_cnt=0;
1029 typedef typename X::GridFunctionSpace GFS;
1032 typedef typename X::template ConstLocalView<LFSCache> XView;
1034 LFS lfs(x.gridFunctionSpace());
1035 LFSCache lfs_cache(lfs);
1038 for(
const auto& cell : elements(gv))
1042 x_view.bind(lfs_cache);
1044 if (x_view[0]>=refine_threshold && cell.level() < max_level)
1049 if (x_view[0]<=coarsen_threshold && cell.level() > min_level)
1057 std::cout <<
"+++ mark_grid: " << refine_cnt <<
" marked for refinement, "
1058 << coarsen_cnt <<
" marked for coarsening" << std::endl;
1062 template<
typename Gr
id,
typename X>
1064 typename X::ElementType coarsen_threshold,
int verbose=0)
1066 typedef typename Grid::LeafGridView GV;
1068 GV gv = grid.leafGridView();
1070 unsigned int coarsen_cnt=0;
1072 typedef typename X::GridFunctionSpace GFS;
1075 typedef typename X::template ConstLocalView<LFSCache> XView;
1077 LFS lfs(x.gridFunctionSpace());
1078 LFSCache lfs_cache(lfs);
1081 for(
const auto& cell : elements(gv))
1085 x_view.bind(lfs_cache);
1087 if (x_view[0]>=refine_threshold)
1092 if (x_view[0]<=coarsen_threshold)
1099 std::cout <<
"+++ mark_grid_for_coarsening: "
1100 << coarsen_cnt <<
" marked for coarsening" << std::endl;
1108 double optimistic_factor;
1109 double coarsen_limit;
1110 double balance_limit;
1115 double refine_fraction_while_refinement;
1116 double coarsen_fraction_while_refinement;
1117 double coarsen_fraction_while_coarsening;
1118 double timestep_decrease_factor;
1119 double timestep_increase_factor;
1129 bool have_decreased_time_step;
1130 bool have_refined_grid;
1133 double accumulated_estimated_error_squared;
1134 double minenergy_rate;
1138 : scaling(16.0), optimistic_factor(1.0), coarsen_limit(0.5), balance_limit(0.33333),
1139 tol(tol_), T(T_), verbose(verbose_), no_adapt(false),
1140 refine_fraction_while_refinement(0.7),
1141 coarsen_fraction_while_refinement(0.2),
1142 coarsen_fraction_while_coarsening(0.2),
1143 timestep_decrease_factor(0.5), timestep_increase_factor(1.5),
1144 accept(false), adapt_dt(false),
adapt_grid(false), newdt(1.0),
1145 have_decreased_time_step(false), have_refined_grid(false),
1146 accumulated_estimated_error_squared(0.0),
1153 timestep_decrease_factor=
s;
1158 timestep_increase_factor=
s;
1163 refine_fraction_while_refinement=
s;
1168 coarsen_fraction_while_refinement=
s;
1173 coarsen_fraction_while_coarsening=
s;
1198 optimistic_factor=
s;
1248 return accumulated_estimated_error_squared;
1254 have_decreased_time_step =
false;
1255 have_refined_grid =
false;
1258 template<
typename GM,
typename X>
1259 void evaluate_estimators (GM& grid,
double time,
double dt,
const X& eta_space,
const X& eta_time,
double energy_timeslab)
1266 double spatial_error = eta_space.one_norm();
1267 double temporal_error = scaling*eta_time.one_norm();
1268 double sum = spatial_error + temporal_error;
1270 double allowed = tol*tol*(energy_timeslab+minenergy_rate*dt);
1271 q_s = spatial_error/sum;
1272 q_t = temporal_error/sum;
1279 <<
" sum/allowed=" << sum/allowed
1281 <<
" estimated error=" << sqrt(accumulated_estimated_error_squared+sum)
1282 <<
" energy_rate=" << energy_timeslab/dt
1289 accumulated_estimated_error_squared += sum;
1290 if (verbose>1) std::cout <<
"+++ no adapt mode" << std::endl;
1299 if (verbose>1) std::cout <<
"+++ accepting time step" << std::endl;
1300 accumulated_estimated_error_squared += sum;
1303 if (sum<coarsen_limit*allowed)
1306 if (q_t<balance_limit)
1309 newdt = timestep_increase_factor*dt;
1311 if (verbose>1) std::cout <<
"+++ spatial error dominates: increase time step" << std::endl;
1315 if (q_s>balance_limit)
1318 newdt = timestep_increase_factor*dt;
1320 if (verbose>1) std::cout <<
"+++ increasing time step" << std::endl;
1323 double eta_refine, eta_coarsen;
1324 if (verbose>1) std::cout <<
"+++ mark grid for coarsening" << std::endl;
1327 coarsen_fraction_while_coarsening,eta_refine,eta_coarsen);
1335 if (q_t<balance_limit)
1338 newdt = timestep_increase_factor*dt;
1340 if (verbose>1) std::cout <<
"+++ spatial error dominates: increase time step" << std::endl;
1347 if (verbose>1) std::cout <<
"+++ will redo time step" << std::endl;
1348 if (q_t>1-balance_limit)
1351 newdt = timestep_decrease_factor*dt;
1353 have_decreased_time_step =
true;
1354 if (verbose>1) std::cout <<
"+++ decreasing time step only" << std::endl;
1358 if (q_t<balance_limit)
1360 if (!have_decreased_time_step)
1363 newdt = timestep_increase_factor*dt;
1365 if (verbose>1) std::cout <<
"+++ increasing time step" << std::endl;
1371 newdt = timestep_decrease_factor*dt;
1373 have_decreased_time_step =
true;
1374 if (verbose>1) std::cout <<
"+++ decreasing time step" << std::endl;
1377 double eta_refine, eta_coarsen;
1378 if (verbose>1) std::cout <<
"+++ BINGO mark grid for refinement and coarsening" << std::endl;
1381 coarsen_fraction_while_refinement,eta_refine,eta_coarsen,0);
1394 #endif // DUNE_PDELAB_ADAPTIVITY_ADAPTIVITY_HH