Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
time_dependent.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/memory_consumption.h>
17 #include <deal.II/base/parallel.h>
18 #include <deal.II/base/utilities.h>
19 
20 #include <deal.II/grid/grid_refinement.h>
21 #include <deal.II/grid/tria.h>
22 #include <deal.II/grid/tria_accessor.h>
23 #include <deal.II/grid/tria_iterator.h>
24 
25 #include <deal.II/lac/vector.h>
26 
27 #include <deal.II/numerics/time_dependent.h>
28 
29 #include <algorithm>
30 #include <functional>
31 #include <numeric>
32 
33 DEAL_II_NAMESPACE_OPEN
34 
36  const unsigned int look_back)
37  : look_ahead(look_ahead)
38  , look_back(look_back)
39 {}
40 
41 
43  const TimeSteppingData &data_dual,
44  const TimeSteppingData &data_postprocess)
45  : sweep_no(numbers::invalid_unsigned_int)
46  , timestepping_data_primal(data_primal)
47  , timestepping_data_dual(data_dual)
48  , timestepping_data_postprocess(data_postprocess)
49 {}
50 
51 
53 {
54  try
55  {
56  while (timesteps.size() != 0)
57  delete_timestep(0);
58  }
59  catch (...)
60  {}
61 }
62 
63 
64 void
66  TimeStepBase * new_timestep)
67 {
68  Assert((std::find(timesteps.begin(), timesteps.end(), position) !=
69  timesteps.end()) ||
70  (position == nullptr),
72  // first insert the new time step
73  // into the doubly linked list
74  // of timesteps
75  if (position == nullptr)
76  {
77  // at the end
78  new_timestep->set_next_timestep(nullptr);
79  if (timesteps.size() > 0)
80  {
81  timesteps.back()->set_next_timestep(new_timestep);
82  new_timestep->set_previous_timestep(timesteps.back());
83  }
84  else
85  new_timestep->set_previous_timestep(nullptr);
86  }
87  else if (position == timesteps[0])
88  {
89  // at the beginning
90  new_timestep->set_previous_timestep(nullptr);
91  if (timesteps.size() > 0)
92  {
93  timesteps[0]->set_previous_timestep(new_timestep);
94  new_timestep->set_next_timestep(timesteps[0]);
95  }
96  else
97  new_timestep->set_next_timestep(nullptr);
98  }
99  else
100  {
101  // inner time step
102  const std::vector<SmartPointer<TimeStepBase, TimeDependent>>::iterator
103  insert_position =
104  std::find(timesteps.begin(), timesteps.end(), position);
105  // check iterators again to satisfy coverity: both insert_position and
106  // insert_position - 1 must be valid iterators
107  Assert(insert_position != timesteps.begin() &&
108  insert_position != timesteps.end(),
109  ExcInternalError());
110 
111  (*(insert_position - 1))->set_next_timestep(new_timestep);
112  new_timestep->set_previous_timestep(*(insert_position - 1));
113  new_timestep->set_next_timestep(*insert_position);
114  (*insert_position)->set_previous_timestep(new_timestep);
115  }
116 
117  // finally enter it into the
118  // array
119  timesteps.insert((position == nullptr ?
120  timesteps.end() :
121  std::find(timesteps.begin(), timesteps.end(), position)),
122  new_timestep);
123 }
124 
125 
126 void
128 {
129  insert_timestep(nullptr, new_timestep);
130 }
131 
132 
133 void
134 TimeDependent::delete_timestep(const unsigned int position)
135 {
136  Assert(position < timesteps.size(), ExcInvalidPosition());
137 
138  // Remember time step object for
139  // later deletion and unlock
140  // SmartPointer
141  TimeStepBase *t = timesteps[position];
142  timesteps[position] = nullptr;
143  // Now delete unsubscribed object
144  delete t;
145 
146  timesteps.erase(timesteps.begin() + position);
147 
148  // reset "next" pointer of previous
149  // time step if possible
150  //
151  // note that if now position==size,
152  // then we deleted the last time step
153  if (position != 0)
154  timesteps[position - 1]->set_next_timestep(
155  (position < timesteps.size()) ?
156  timesteps[position] :
158 
159  // same for "previous" pointer of next
160  // time step
161  if (position < timesteps.size())
162  timesteps[position]->set_previous_timestep(
163  (position != 0) ? timesteps[position - 1] :
165 }
166 
167 
168 void
170 {
171  do_loop(
172  [](TimeStepBase *const time_step) { time_step->init_for_primal_problem(); },
173  [](TimeStepBase *const time_step) { time_step->solve_primal_problem(); },
175  forward);
176 }
177 
178 
179 void
181 {
182  do_loop(
183  [](TimeStepBase *const time_step) { time_step->init_for_dual_problem(); },
184  [](TimeStepBase *const time_step) { time_step->init_for_dual_problem(); },
186  backward);
187 }
188 
189 
190 void
192 {
193  do_loop(
194  [](TimeStepBase *const time_step) { time_step->init_for_postprocessing(); },
195  [](TimeStepBase *const time_step) { time_step->postprocess_timestep(); },
197  forward);
198 }
199 
200 
201 
202 void
203 TimeDependent::start_sweep(const unsigned int s)
204 {
205  sweep_no = s;
206 
207  // reset the number each
208  // time step has, since some time
209  // steps might have been added since
210  // the last time we visited them
211  //
212  // also set the sweep we will
213  // process in the sequel
214  for (unsigned int step = 0; step < timesteps.size(); ++step)
215  {
216  timesteps[step]->set_timestep_no(step);
217  timesteps[step]->set_sweep_no(sweep_no);
218  }
219 
220  for (const auto &timestep : timesteps)
221  timestep->start_sweep();
222 }
223 
224 
225 
226 void
228 {
230  timesteps.size(),
231  [this](const unsigned int begin,
232  const unsigned int end) {
233  this->end_sweep(begin, end);
234  },
235  1);
236 }
237 
238 
239 
240 void
241 TimeDependent::end_sweep(const unsigned int begin, const unsigned int end)
242 {
243  for (unsigned int step = begin; step < end; ++step)
244  timesteps[step]->end_sweep();
245 }
246 
247 
248 
249 std::size_t
251 {
252  std::size_t mem =
257  for (const auto &timestep : timesteps)
258  mem += MemoryConsumption::memory_consumption(*timestep);
259 
260  return mem;
261 }
262 
263 
264 
265 /* --------------------------------------------------------------------- */
266 
267 
268 TimeStepBase::TimeStepBase(const double time)
269  : previous_timestep(nullptr)
270  , next_timestep(nullptr)
271  , sweep_no(numbers::invalid_unsigned_int)
272  , timestep_no(numbers::invalid_unsigned_int)
273  , time(time)
274  , next_action(numbers::invalid_unsigned_int)
275 {}
276 
277 
278 
279 void
280 TimeStepBase::wake_up(const unsigned int)
281 {}
282 
283 
284 
285 void
286 TimeStepBase::sleep(const unsigned)
287 {}
288 
289 
290 
291 void
293 {}
294 
295 
296 
297 void
299 {}
300 
301 
302 
303 void
305 {
307 }
308 
309 
310 
311 void
313 {
315 }
316 
317 
318 
319 void
321 {
323 }
324 
325 
326 
327 void
329 {
330  Assert(false, ExcPureFunctionCalled());
331 }
332 
333 
334 
335 void
337 {
338  Assert(false, ExcPureFunctionCalled());
339 }
340 
341 
342 
343 double
345 {
346  return time;
347 }
348 
349 
350 
351 unsigned int
353 {
354  return timestep_no;
355 }
356 
357 
358 
359 double
361 {
362  Assert(previous_timestep != nullptr,
363  ExcMessage("The backward time step cannot be computed because "
364  "there is no previous time step."));
365  return time - previous_timestep->time;
366 }
367 
368 
369 
370 double
372 {
373  Assert(next_timestep != nullptr,
374  ExcMessage("The forward time step cannot be computed because "
375  "there is no next time step."));
376  return next_timestep->time - time;
377 }
378 
379 
380 
381 void
383 {
384  previous_timestep = previous;
385 }
386 
387 
388 
389 void
391 {
392  next_timestep = next;
393 }
394 
395 
396 
397 void
398 TimeStepBase::set_timestep_no(const unsigned int step_no)
399 {
400  timestep_no = step_no;
401 }
402 
403 
404 
405 void
406 TimeStepBase::set_sweep_no(const unsigned int sweep)
407 {
408  sweep_no = sweep;
409 }
410 
411 
412 
413 std::size_t
415 {
416  // only simple data types
417  return sizeof(*this);
418 }
419 
420 
421 
422 template <int dim>
424  : TimeStepBase(0)
425  , tria(nullptr, typeid(*this).name())
426  , coarse_grid(nullptr, typeid(*this).name())
427  , flags()
428  , refinement_flags(0)
429 {
430  Assert(false, ExcPureFunctionCalled());
431 }
432 
433 
434 
435 template <int dim>
437  const double time,
439  const Flags & flags,
440  const RefinementFlags & refinement_flags)
441  : TimeStepBase(time)
442  , tria(nullptr, typeid(*this).name())
443  , coarse_grid(&coarse_grid, typeid(*this).name())
444  , flags(flags)
446 {}
447 
448 
449 
450 template <int dim>
452 {
453  if (!flags.delete_and_rebuild_tria)
454  {
456  tria = nullptr;
457  delete t;
458  }
459  else
460  AssertNothrow(tria == nullptr, ExcInternalError());
461 
462  coarse_grid = nullptr;
463 }
464 
465 
466 
467 template <int dim>
468 void
469 TimeStepBase_Tria<dim>::wake_up(const unsigned int wakeup_level)
470 {
471  TimeStepBase::wake_up(wakeup_level);
472 
473  if (wakeup_level == flags.wakeup_level_to_build_grid)
474  if (flags.delete_and_rebuild_tria || !tria)
475  restore_grid();
476 }
477 
478 
479 
480 template <int dim>
481 void
482 TimeStepBase_Tria<dim>::sleep(const unsigned int sleep_level)
483 {
484  if (sleep_level == flags.sleep_level_to_delete_grid)
485  {
486  Assert(tria != nullptr, ExcInternalError());
487 
488  if (flags.delete_and_rebuild_tria)
489  {
491  tria = nullptr;
492  delete t;
493  }
494  }
495 
496  TimeStepBase::sleep(sleep_level);
497 }
498 
499 
500 
501 template <int dim>
502 void
504 {
505  // for any of the non-initial grids
506  // store the refinement flags
507  refine_flags.emplace_back();
508  coarsen_flags.emplace_back();
511 }
512 
513 
514 
515 template <int dim>
516 void
518 {
519  Assert(tria == nullptr, ExcGridNotDeleted());
520  Assert(refine_flags.size() == coarsen_flags.size(), ExcInternalError());
521 
522  // create a virgin triangulation and
523  // set it to a copy of the coarse grid
524  tria = new Triangulation<dim>();
525  tria->copy_triangulation(*coarse_grid);
526 
527  // for each of the previous refinement
528  // sweeps
529  for (unsigned int previous_sweep = 0; previous_sweep < refine_flags.size();
530  ++previous_sweep)
531  {
532  // get flags
533  tria->load_refine_flags(refine_flags[previous_sweep]);
534  tria->load_coarsen_flags(coarsen_flags[previous_sweep]);
535 
536  // limit refinement depth if the user
537  // desired so
538  // if (flags.max_refinement_level != 0)
539  // {
540  // typename Triangulation<dim>::active_cell_iterator cell, endc;
541  // for (cell = tria->begin_active(),
542  // endc = tria->end();
543  // cell!=endc; ++cell)
544  // if (static_cast<unsigned int>(cell->level()) >=
545  // flags.max_refinement_level)
546  // cell->clear_refine_flag();
547  // }
548 
549  tria->execute_coarsening_and_refinement();
550  }
551 }
552 
553 
554 
555 // have a few helper functions
556 namespace
557 {
558  template <int dim>
559  void
560  mirror_refinement_flags(
561  const typename Triangulation<dim>::cell_iterator &new_cell,
562  const typename Triangulation<dim>::cell_iterator &old_cell)
563  {
564  // mirror the refinement
565  // flags from the present time level to
566  // the previous if the dual problem was
567  // used for the refinement, since the
568  // error is computed on a time-space cell
569  //
570  // we don't mirror the coarsening flags
571  // since we want stronger refinement. if
572  // this was the wrong decision, the error
573  // on the child cells of the previous
574  // time slab will indicate coarsening
575  // in the next iteration, so this is not
576  // so dangerous here.
577  //
578  // also, we only have to check whether
579  // the present cell flagged for
580  // refinement and the previous one is on
581  // the same level and also active. If it
582  // already has children, then there is
583  // no problem at all, if it is on a lower
584  // level than the present one, then it
585  // will be refined below anyway.
586  if (new_cell->active())
587  {
588  if (new_cell->refine_flag_set() && old_cell->active())
589  {
590  if (old_cell->coarsen_flag_set())
591  old_cell->clear_coarsen_flag();
592 
593  old_cell->set_refine_flag();
594  }
595 
596  return;
597  }
598 
599  if (old_cell->has_children() && new_cell->has_children())
600  {
601  Assert(old_cell->n_children() == new_cell->n_children(),
603  for (unsigned int c = 0; c < new_cell->n_children(); ++c)
604  ::mirror_refinement_flags<dim>(new_cell->child(c),
605  old_cell->child(c));
606  }
607  }
608 
609 
610 
611  template <int dim>
612  bool
613  adapt_grid_cells(const typename Triangulation<dim>::cell_iterator &cell1,
614  const typename Triangulation<dim>::cell_iterator &cell2)
615  {
616  if (cell2->has_children() && cell1->has_children())
617  {
618  bool grids_changed = false;
619 
620  Assert(cell2->n_children() == cell1->n_children(), ExcNotImplemented());
621  for (unsigned int c = 0; c < cell1->n_children(); ++c)
622  grids_changed |=
623  ::adapt_grid_cells<dim>(cell1->child(c), cell2->child(c));
624  return grids_changed;
625  }
626 
627 
628  if (!cell1->has_children() && !cell2->has_children())
629  // none of the two have children, so
630  // make sure that not one is flagged
631  // for refinement and the other for
632  // coarsening
633  {
634  if (cell1->refine_flag_set() && cell2->coarsen_flag_set())
635  {
636  cell2->clear_coarsen_flag();
637  return true;
638  }
639  else if (cell1->coarsen_flag_set() && cell2->refine_flag_set())
640  {
641  cell1->clear_coarsen_flag();
642  return true;
643  }
644 
645  return false;
646  }
647 
648 
649  if (cell1->has_children() && !cell2->has_children())
650  // cell1 has children, cell2 has not
651  // -> cell2 needs to be refined if any
652  // of cell1's children is flagged
653  // for refinement. None of them should
654  // be refined further, since then in the
655  // last round something must have gone
656  // wrong
657  //
658  // if cell2 was flagged for coarsening,
659  // we need to clear that flag in any
660  // case. The only exception would be
661  // if all children of cell1 were
662  // flagged for coarsening, but rules
663  // for coarsening are so complicated
664  // that we will not attempt to cover
665  // them. Rather accept one cell which
666  // is not coarsened...
667  {
668  bool changed_grid = false;
669  if (cell2->coarsen_flag_set())
670  {
671  cell2->clear_coarsen_flag();
672  changed_grid = true;
673  }
674 
675  if (!cell2->refine_flag_set())
676  for (unsigned int c = 0; c < cell1->n_children(); ++c)
677  if (cell1->child(c)->refine_flag_set() ||
678  cell1->child(c)->has_children())
679  {
680  cell2->set_refine_flag();
681  changed_grid = true;
682  break;
683  }
684  return changed_grid;
685  }
686 
687  if (!cell1->has_children() && cell2->has_children())
688  // same thing, other way round...
689  {
690  bool changed_grid = false;
691  if (cell1->coarsen_flag_set())
692  {
693  cell1->clear_coarsen_flag();
694  changed_grid = true;
695  }
696 
697  if (!cell1->refine_flag_set())
698  for (unsigned int c = 0; c < cell2->n_children(); ++c)
699  if (cell2->child(c)->refine_flag_set() ||
700  cell2->child(c)->has_children())
701  {
702  cell1->set_refine_flag();
703  changed_grid = true;
704  break;
705  }
706  return changed_grid;
707  }
708 
709  Assert(false, ExcInternalError());
710  return false;
711  }
712 
713 
714 
715  template <int dim>
716  bool
717  adapt_grids(Triangulation<dim> &tria1, Triangulation<dim> &tria2)
718  {
719  bool grids_changed = false;
720 
721  typename Triangulation<dim>::cell_iterator cell1 = tria1.begin(),
722  cell2 = tria2.begin();
723  typename Triangulation<dim>::cell_iterator endc;
724  endc = (tria1.n_levels() == 1 ?
725  typename Triangulation<dim>::cell_iterator(tria1.end()) :
726  tria1.begin(1));
727  for (; cell1 != endc; ++cell1, ++cell2)
728  grids_changed |= ::adapt_grid_cells<dim>(cell1, cell2);
729 
730  return grids_changed;
731  }
732 } // namespace
733 
734 
735 template <int dim>
736 void
737 TimeStepBase_Tria<dim>::refine_grid(const RefinementData refinement_data)
738 {
739  Vector<float> criteria;
741 
742  // copy the following two values since
743  // we may need modified values in the
744  // process of this function
745  double refinement_threshold = refinement_data.refinement_threshold,
746  coarsening_threshold = refinement_data.coarsening_threshold;
747 
748  // prepare an array where the criteria
749  // are stored in a sorted fashion
750  // we need this if cell number correction
751  // is switched on.
752  // the criteria are sorted in ascending
753  // order
754  // only fill it when needed
755  Vector<float> sorted_criteria;
756  // two pointers into this array denoting
757  // the position where the two thresholds
758  // are assumed
759  Vector<float>::const_iterator p_refinement_threshold = nullptr,
760  p_coarsening_threshold = nullptr;
761 
762 
763  // if we are to do some cell number
764  // correction steps, we have to find out
765  // which further cells (beyond
766  // refinement_threshold) to refine in case
767  // we need more cells, and which cells
768  // to not refine in case we need less cells
769  // (or even to coarsen, if necessary). to
770  // this end, we first define pointers into
771  // a sorted array of criteria pointing
772  // to the thresholds of refinement or
773  // coarsening; moving these pointers amounts
774  // to changing the threshold such that the
775  // number of cells flagged for refinement
776  // or coarsening would be changed by one
777  if ((timestep_no != 0) &&
778  (sweep_no >= refinement_flags.first_sweep_with_correction) &&
779  (refinement_flags.cell_number_correction_steps > 0))
780  {
781  sorted_criteria = criteria;
782  std::sort(sorted_criteria.begin(), sorted_criteria.end());
783  p_refinement_threshold =
784  Utilities::lower_bound(sorted_criteria.begin(),
785  sorted_criteria.end(),
786  static_cast<float>(refinement_threshold));
787  p_coarsening_threshold =
788  std::upper_bound(sorted_criteria.begin(),
789  sorted_criteria.end(),
790  static_cast<float>(coarsening_threshold));
791  }
792 
793 
794  // actually flag cells the first time
795  GridRefinement::refine(*tria, criteria, refinement_threshold);
796  GridRefinement::coarsen(*tria, criteria, coarsening_threshold);
797 
798  // store this number for the following
799  // since its computation is rather
800  // expensive and since it doesn't change
801  const unsigned int n_active_cells = tria->n_active_cells();
802 
803  // if not on first time level: try to
804  // adjust the number of resulting
805  // cells to those on the previous
806  // time level. Only do the cell number
807  // correction for higher sweeps and if
808  // there are sufficiently many cells
809  // already to avoid "grid stall" i.e.
810  // that the grid's evolution is hindered
811  // by the correction (this usually
812  // happens if there are very few cells,
813  // since then the number of cells touched
814  // by the correction step may exceed the
815  // number of cells which are flagged for
816  // refinement; in this case it often
817  // happens that the number of cells
818  // does not grow between sweeps, which
819  // clearly is not the wanted behaviour)
820  //
821  // however, if we do not do anything, we
822  // can get into trouble somewhen later.
823  // therefore, we also use the correction
824  // step for the first sweep or if the
825  // number of cells is between 100 and 300
826  // (unlike in the first version of the
827  // algorithm), but relax the conditions
828  // for the correction to allow deviations
829  // which are three times as high than
830  // allowed (sweep==1 || cell number<200)
831  // or twice as high (sweep==2 ||
832  // cell number<300). Also, since
833  // refinement never does any harm other
834  // than increased work, we allow for
835  // arbitrary growth of cell number if
836  // the estimated cell number is below
837  // 200.
838  //
839  // repeat this loop several times since
840  // the first estimate may not be totally
841  // correct
842  if ((timestep_no != 0) &&
843  (sweep_no >= refinement_flags.first_sweep_with_correction))
844  for (unsigned int loop = 0;
845  loop < refinement_flags.cell_number_correction_steps;
846  ++loop)
847  {
848  Triangulation<dim> *previous_tria =
849  dynamic_cast<const TimeStepBase_Tria<dim> *>(previous_timestep)->tria;
850 
851  // do one adaption step if desired
852  // (there are more coming below then
853  // also)
854  if (refinement_flags.adapt_grids)
855  ::adapt_grids<dim>(*previous_tria, *tria);
856 
857  // perform flagging of cells
858  // needed to regularize the
859  // triangulation
861  previous_tria->prepare_coarsening_and_refinement();
862 
863 
864  // now count the number of elements
865  // which will result on the previous
866  // grid after it will be refined. The
867  // number which will really result
868  // should be approximately that that we
869  // compute here, since we already
870  // performed most of the prepare*
871  // steps for the previous grid
872  //
873  // use a double value since for each
874  // four cells (in 2D) that we flagged
875  // for coarsening we result in one
876  // new. but since we loop over flagged
877  // cells, we have to subtract 3/4 of
878  // a cell for each flagged cell
880  Assert(!previous_tria->get_anisotropic_refinement_flag(),
882  double previous_cells = previous_tria->n_active_cells();
883  typename Triangulation<dim>::active_cell_iterator cell, endc;
884  cell = previous_tria->begin_active();
885  endc = previous_tria->end();
886  for (; cell != endc; ++cell)
887  if (cell->refine_flag_set())
888  previous_cells += (GeometryInfo<dim>::max_children_per_cell - 1);
889  else if (cell->coarsen_flag_set())
890  previous_cells -= static_cast<double>(
893 
894  // @p{previous_cells} now gives the
895  // number of cells which would result
896  // from the flags on the previous grid
897  // if we refined it now. However, some
898  // more flags will be set when we adapt
899  // the previous grid with this one
900  // after the flags have been set for
901  // this time level; on the other hand,
902  // we don't account for this, since the
903  // number of cells on this time level
904  // will be changed afterwards by the
905  // same way, when it is adapted to the
906  // next time level
907 
908  // now estimate the number of cells which
909  // will result on this level
910  double estimated_cells = n_active_cells;
911  cell = tria->begin_active();
912  endc = tria->end();
913  for (; cell != endc; ++cell)
914  if (cell->refine_flag_set())
915  estimated_cells += (GeometryInfo<dim>::max_children_per_cell - 1);
916  else if (cell->coarsen_flag_set())
917  estimated_cells -= static_cast<double>(
920 
921  // calculate the allowed delta in
922  // cell numbers; be more lenient
923  // if there are few cells
924  double delta_up = refinement_flags.cell_number_corridor_top,
925  delta_down = refinement_flags.cell_number_corridor_bottom;
926 
927  const std::vector<std::pair<unsigned int, double>> &relaxations =
928  (sweep_no >= refinement_flags.correction_relaxations.size() ?
929  refinement_flags.correction_relaxations.back() :
930  refinement_flags.correction_relaxations[sweep_no]);
931  for (unsigned int r = 0; r != relaxations.size(); ++r)
932  if (n_active_cells < relaxations[r].first)
933  {
934  delta_up *= relaxations[r].second;
935  delta_down *= relaxations[r].second;
936  break;
937  }
938 
939  // now, if the number of estimated
940  // cells exceeds the number of cells
941  // on the old time level by more than
942  // delta: cut the top threshold
943  //
944  // note that for each cell that
945  // we unflag we have to diminish the
946  // estimated number of cells by
947  // @p{children_per_cell}.
948  if (estimated_cells > previous_cells * (1. + delta_up))
949  {
950  // only limit the cell number
951  // if there will not be less
952  // than some number of cells
953  //
954  // also note that when using the
955  // dual estimator, the initial
956  // time level is not refined
957  // on its own, so we may not
958  // limit the number of the second
959  // time level on the basis of
960  // the initial one; since for
961  // the dual estimator, we
962  // mirror the refinement
963  // flags, the initial level
964  // will be passively refined
965  // later on.
966  if (estimated_cells > refinement_flags.min_cells_for_correction)
967  {
968  // number of cells by which the
969  // new grid is to be diminished
970  double delta_cells =
971  estimated_cells - previous_cells * (1. + delta_up);
972 
973  // if we need to reduce the
974  // number of cells, we need
975  // to raise the thresholds,
976  // i.e. move ahead in the
977  // sorted array, since this
978  // is sorted in ascending
979  // order. do so by removing
980  // cells tagged for refinement
981 
982  for (unsigned int i = 0; i < delta_cells;
984  if (p_refinement_threshold != sorted_criteria.end())
985  ++p_refinement_threshold;
986  else
987  break;
988  }
989  else
990  // too many cells, but we
991  // won't do anything about
992  // that
993  break;
994  }
995  else
996  // likewise: if the estimated number
997  // of cells is less than 90 per cent
998  // of those at the previous time level:
999  // raise threshold by refining
1000  // additional cells. if we start to
1001  // run into the area of cells
1002  // which are to be coarsened, we
1003  // raise the limit for these too
1004  if (estimated_cells < previous_cells * (1. - delta_down))
1005  {
1006  // number of cells by which the
1007  // new grid is to be enlarged
1008  double delta_cells =
1009  previous_cells * (1. - delta_down) - estimated_cells;
1010  // heuristics: usually, if we
1011  // add @p{delta_cells} to the
1012  // present state, we end up
1013  // with much more than only
1014  // (1-delta_down)*prev_cells
1015  // because of the effect of
1016  // regularization and because
1017  // of adaption to the
1018  // following grid. Therefore,
1019  // if we are not in the last
1020  // correction loop, we try not
1021  // to add as many cells as seem
1022  // necessary at first and hope
1023  // to get closer to the limit
1024  // this way. Only in the last
1025  // loop do we have to take the
1026  // full number to guarantee the
1027  // wanted result.
1028  //
1029  // The value 0.9 is taken from
1030  // practice, as the additional
1031  // number of cells introduced
1032  // by regularization is
1033  // approximately 10 per cent
1034  // of the flagged cells.
1035  if (loop != refinement_flags.cell_number_correction_steps - 1)
1036  delta_cells *= 0.9;
1037 
1038  // if more cells need to be
1039  // refined, we need to lower
1040  // the thresholds, i.e. to
1041  // move to the beginning
1042  // of sorted_criteria, which is
1043  // sorted in ascending order
1044  for (unsigned int i = 0; i < delta_cells;
1046  if (p_refinement_threshold != p_coarsening_threshold)
1047  --refinement_threshold;
1048  else if (p_coarsening_threshold != sorted_criteria.begin())
1049  --p_coarsening_threshold, --p_refinement_threshold;
1050  else
1051  break;
1052  }
1053  else
1054  // estimated cell number is ok,
1055  // stop correction steps
1056  break;
1057 
1058  if (p_refinement_threshold == sorted_criteria.end())
1059  {
1060  Assert(p_coarsening_threshold != p_refinement_threshold,
1061  ExcInternalError());
1062  --p_refinement_threshold;
1063  }
1064 
1065  coarsening_threshold = *p_coarsening_threshold;
1066  refinement_threshold = *p_refinement_threshold;
1067 
1068  if (coarsening_threshold >= refinement_threshold)
1069  coarsening_threshold = 0.999 * refinement_threshold;
1070 
1071  // now that we have re-adjusted
1072  // thresholds: clear all refine and
1073  // coarsening flags and do it all
1074  // over again
1075  cell = tria->begin_active();
1076  endc = tria->end();
1077  for (; cell != endc; ++cell)
1078  {
1079  cell->clear_refine_flag();
1080  cell->clear_coarsen_flag();
1081  }
1082 
1083 
1084  // flag cells finally
1085  GridRefinement::refine(*tria, criteria, refinement_threshold);
1086  GridRefinement::coarsen(*tria, criteria, coarsening_threshold);
1087  }
1088 
1089  // if step number is greater than
1090  // one: adapt this and the previous
1091  // grid to each other. Don't do so
1092  // for the initial grid because
1093  // it is always taken to be the first
1094  // grid and needs therefore no
1095  // treatment of its own.
1096  if ((timestep_no >= 1) && (refinement_flags.adapt_grids))
1097  {
1098  Triangulation<dim> *previous_tria =
1099  dynamic_cast<const TimeStepBase_Tria<dim> *>(previous_timestep)->tria;
1100  Assert(previous_tria != nullptr, ExcInternalError());
1101 
1102  // if we used the dual estimator, we
1103  // computed the error information on
1104  // a time slab, rather than on a level
1105  // of its own. we then mirror the
1106  // refinement flags we determined for
1107  // the present level to the previous
1108  // one
1109  //
1110  // do this mirroring only, if cell number
1111  // adjustment is on, since otherwise
1112  // strange things may happen
1113  if (refinement_flags.mirror_flags_to_previous_grid)
1114  {
1115  ::adapt_grids<dim>(*previous_tria, *tria);
1116 
1117  typename Triangulation<dim>::cell_iterator old_cell, new_cell, endc;
1118  old_cell = previous_tria->begin(0);
1119  new_cell = tria->begin(0);
1120  endc = tria->end(0);
1121  for (; new_cell != endc; ++new_cell, ++old_cell)
1122  ::mirror_refinement_flags<dim>(new_cell, old_cell);
1123  }
1124 
1126  previous_tria->prepare_coarsening_and_refinement();
1127 
1128  // adapt present and previous grids
1129  // to each other: flag additional
1130  // cells to avoid the previous grid
1131  // to have cells refined twice more
1132  // than the present one and vica versa.
1133  ::adapt_grids<dim>(*previous_tria, *tria);
1134  }
1135 }
1136 
1137 
1138 
1139 template <int dim>
1140 void
1142 {
1144 }
1145 
1146 
1147 
1148 template <int dim>
1149 std::size_t
1151 {
1152  return (TimeStepBase::memory_consumption() + sizeof(tria) +
1153  MemoryConsumption::memory_consumption(coarse_grid) + sizeof(flags) +
1154  sizeof(refinement_flags) +
1157 }
1158 
1159 
1160 
1161 template <int dim>
1163  : delete_and_rebuild_tria(false)
1164  , wakeup_level_to_build_grid(0)
1165  , sleep_level_to_delete_grid(0)
1166 {
1167  Assert(false, ExcInternalError());
1168 }
1169 
1170 
1171 
1172 template <int dim>
1174  const bool delete_and_rebuild_tria,
1175  const unsigned int wakeup_level_to_build_grid,
1176  const unsigned int sleep_level_to_delete_grid)
1177  : delete_and_rebuild_tria(delete_and_rebuild_tria)
1178  , wakeup_level_to_build_grid(wakeup_level_to_build_grid)
1179  , sleep_level_to_delete_grid(sleep_level_to_delete_grid)
1180 {
1181  // Assert (!delete_and_rebuild_tria || (wakeup_level_to_build_grid>=1),
1182  // ExcInvalidParameter(wakeup_level_to_build_grid));
1183  // Assert (!delete_and_rebuild_tria || (sleep_level_to_delete_grid>=1),
1184  // ExcInvalidParameter(sleep_level_to_delete_grid));
1185 }
1186 
1187 
1188 template <int dim>
1191  1, // one element, denoting the first and all subsequent sweeps
1192  std::vector<std::pair<unsigned int, double>>(1, // one element, denoting the
1193  // upper bound for the
1194  // following relaxation
1195  std::make_pair(0U, 0.)));
1196 
1197 
1198 template <int dim>
1200  const unsigned int max_refinement_level,
1201  const unsigned int first_sweep_with_correction,
1202  const unsigned int min_cells_for_correction,
1203  const double cell_number_corridor_top,
1204  const double cell_number_corridor_bottom,
1205  const CorrectionRelaxations &correction_relaxations,
1206  const unsigned int cell_number_correction_steps,
1207  const bool mirror_flags_to_previous_grid,
1208  const bool adapt_grids)
1209  : max_refinement_level(max_refinement_level)
1210  , first_sweep_with_correction(first_sweep_with_correction)
1211  , min_cells_for_correction(min_cells_for_correction)
1212  , cell_number_corridor_top(cell_number_corridor_top)
1213  , cell_number_corridor_bottom(cell_number_corridor_bottom)
1214  , correction_relaxations(correction_relaxations.size() != 0 ?
1215  correction_relaxations :
1216  default_correction_relaxations)
1217  , cell_number_correction_steps(cell_number_correction_steps)
1218  , mirror_flags_to_previous_grid(mirror_flags_to_previous_grid)
1219  , adapt_grids(adapt_grids)
1220 {
1221  Assert(cell_number_corridor_top >= 0,
1222  ExcInvalidValue(cell_number_corridor_top));
1223  Assert(cell_number_corridor_bottom >= 0,
1224  ExcInvalidValue(cell_number_corridor_bottom));
1225  Assert(cell_number_corridor_bottom <= 1,
1226  ExcInvalidValue(cell_number_corridor_bottom));
1227 }
1228 
1229 
1230 template <int dim>
1232  const double _refinement_threshold,
1233  const double _coarsening_threshold)
1234  : refinement_threshold(_refinement_threshold)
1235  ,
1236  // in some rare cases it may happen that
1237  // both thresholds are the same (e.g. if
1238  // there are many cells with the same
1239  // error indicator). That would mean that
1240  // all cells will be flagged for
1241  // refinement or coarsening, but some will
1242  // be flagged for both, namely those for
1243  // which the indicator equals the
1244  // thresholds. This is forbidden, however.
1245  //
1246  // In some rare cases with very few cells
1247  // we also could get integer round off
1248  // errors and get problems with
1249  // the top and bottom fractions.
1250  //
1251  // In these case we arbitrarily reduce the
1252  // bottom threshold by one permille below
1253  // the top threshold
1254  coarsening_threshold((_coarsening_threshold == _refinement_threshold ?
1255  _coarsening_threshold :
1256  0.999 * _coarsening_threshold))
1257 {
1260  // allow both thresholds to be zero,
1261  // since this is needed in case all indicators
1262  // are zero
1264  ((coarsening_threshold == 0) && (refinement_threshold == 0)),
1266 }
1267 
1268 
1269 
1270 /*-------------- Explicit Instantiations -------------------------------*/
1271 #include "time_dependent.inst"
1272 
1273 
1274 DEAL_II_NAMESPACE_CLOSE
void set_sweep_no(const unsigned int sweep_no)
virtual std::size_t memory_consumption() const
virtual void solve_primal_problem()=0
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1063
virtual void sleep(const unsigned int)
virtual void start_sweep()
unsigned int n_active_cells() const
Definition: tria.cc:12601
static ::ExceptionBase & ExcPureFunctionCalled()
virtual void copy_triangulation(const Triangulation< dim, spacedim > &other_tria)
Definition: tria.cc:10477
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:446
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1471
void refine_grid(const RefinementData data)
unsigned int next_action
static ::ExceptionBase & ExcInvalidValue(double arg1)
TimeDependent(const TimeSteppingData &data_primal, const TimeSteppingData &data_dual, const TimeSteppingData &data_postprocess)
virtual void solve_dual_problem()
void set_timestep_no(const unsigned int step_no)
void solve_primal_problem()
void add_timestep(TimeStepBase *new_timestep)
void solve_dual_problem()
const TimeStepBase * next_timestep
const TimeSteppingData timestepping_data_postprocess
const unsigned int wakeup_level_to_build_grid
virtual void wake_up(const unsigned int wakeup_level) override
double get_forward_timestep() const
void set_next_timestep(const TimeStepBase *next)
double get_time() const
const TimeStepBase * previous_timestep
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:11939
const TimeSteppingData timestepping_data_dual
cell_iterator begin(const unsigned int level=0) const
Definition: tria.cc:11919
void apply_to_subranges(const RangeType &begin, const typename identity< RangeType >::type &end, const Function &f, const unsigned int grainsize)
Definition: parallel.h:447
unsigned int n_levels() const
SmartPointer< const Triangulation< dim, dim >, TimeStepBase_Tria< dim > > coarse_grid
virtual std::size_t memory_consumption() const override
void set_previous_timestep(const TimeStepBase *previous)
cell_iterator end() const
Definition: tria.cc:12005
unsigned int timestep_no
unsigned int get_timestep_no() const
void do_loop(InitFunctionObject init_function, LoopFunctionObject loop_function, const TimeSteppingData &timestepping_data, const Direction direction)
const RefinementFlags refinement_flags
virtual bool prepare_coarsening_and_refinement()
Definition: tria.cc:14025
static ::ExceptionBase & ExcMessage(std::string arg1)
std::vector< std::vector< std::pair< unsigned int, double > >> CorrectionRelaxations
std::size_t memory_consumption() const
virtual void init_for_refinement()
RefinementData(const double refinement_threshold, const double coarsening_threshold=0)
#define Assert(cond, exc)
Definition: exceptions.h:1407
const double time
bool get_anisotropic_refinement_flag() const
Definition: tria.cc:10925
std::vector< SmartPointer< TimeStepBase, TimeDependent > > timesteps
virtual void get_tria_refinement_criteria(Vector< float > &criteria) const =0
virtual ~TimeDependent()
const TimeSteppingData timestepping_data_primal
static ::ExceptionBase & ExcInvalidPosition()
virtual void start_sweep(const unsigned int sweep_no)
void save_coarsen_flags(std::ostream &out) const
Definition: tria.cc:10879
void save_refine_flags(std::ostream &out) const
Definition: tria.cc:10811
virtual void init_for_postprocessing()
static ::ExceptionBase & ExcGridNotDeleted()
void coarsen(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold)
static CorrectionRelaxations default_correction_relaxations
static ::ExceptionBase & ExcInvalidValue(double arg1)
double get_backward_timestep() const
const unsigned int sleep_level_to_delete_grid
void delete_timestep(const unsigned int position)
virtual void wake_up(const unsigned int)
RefinementFlags(const unsigned int max_refinement_level=0, const unsigned int first_sweep_with_correction=0, const unsigned int min_cells_for_correction=0, const double cell_number_corridor_top=(1<< dim), const double cell_number_corridor_bottom=1, const CorrectionRelaxations &correction_relaxations=CorrectionRelaxations(), const unsigned int cell_number_correction_steps=0, const bool mirror_flags_to_previous_grid=false, const bool adapt_grids=false)
typename TimeStepBase_Tria_Flags::Flags< dim > Flags
virtual void init_for_dual_problem()
virtual void init_for_primal_problem()
unsigned int sweep_no
static ::ExceptionBase & ExcNotImplemented()
void insert_timestep(const TimeStepBase *position, TimeStepBase *new_timestep)
std::vector< std::vector< bool > > refine_flags
virtual void end_sweep()
SmartPointer< Triangulation< dim, dim >, TimeStepBase_Tria< dim > > tria
TriaIterator< CellAccessor< dim, spacedim > > cell_iterator
Definition: tria.h:1529
std::vector< std::vector< bool > > coarsen_flags
virtual void postprocess_timestep()
virtual ~TimeStepBase_Tria() override
virtual void end_sweep()
unsigned int sweep_no
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
virtual void sleep(const unsigned int) override
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
TimeStepBase(const double time)
TimeSteppingData(const unsigned int look_ahead, const unsigned int look_back)
static ::ExceptionBase & ExcInternalError()