Reference documentation for deal.II version Git b84270a1d4 2020-01-17 16:21:02 -0500
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 (const auto &cell : tria->active_cell_iterators())
542  // if (static_cast<unsigned int>(cell->level()) >=
543  // flags.max_refinement_level)
544  // cell->clear_refine_flag();
545  // }
546 
547  tria->execute_coarsening_and_refinement();
548  }
549 }
550 
551 
552 
553 // have a few helper functions
554 namespace
555 {
556  template <int dim>
557  void
558  mirror_refinement_flags(
559  const typename Triangulation<dim>::cell_iterator &new_cell,
560  const typename Triangulation<dim>::cell_iterator &old_cell)
561  {
562  // mirror the refinement
563  // flags from the present time level to
564  // the previous if the dual problem was
565  // used for the refinement, since the
566  // error is computed on a time-space cell
567  //
568  // we don't mirror the coarsening flags
569  // since we want stronger refinement. if
570  // this was the wrong decision, the error
571  // on the child cells of the previous
572  // time slab will indicate coarsening
573  // in the next iteration, so this is not
574  // so dangerous here.
575  //
576  // also, we only have to check whether
577  // the present cell flagged for
578  // refinement and the previous one is on
579  // the same level and also active. If it
580  // already has children, then there is
581  // no problem at all, if it is on a lower
582  // level than the present one, then it
583  // will be refined below anyway.
584  if (new_cell->is_active())
585  {
586  if (new_cell->refine_flag_set() && old_cell->is_active())
587  {
588  if (old_cell->coarsen_flag_set())
589  old_cell->clear_coarsen_flag();
590 
591  old_cell->set_refine_flag();
592  }
593 
594  return;
595  }
596 
597  if (old_cell->has_children() && new_cell->has_children())
598  {
599  Assert(old_cell->n_children() == new_cell->n_children(),
601  for (unsigned int c = 0; c < new_cell->n_children(); ++c)
602  ::mirror_refinement_flags<dim>(new_cell->child(c),
603  old_cell->child(c));
604  }
605  }
606 
607 
608 
609  template <int dim>
610  bool
611  adapt_grid_cells(const typename Triangulation<dim>::cell_iterator &cell1,
612  const typename Triangulation<dim>::cell_iterator &cell2)
613  {
614  if (cell2->has_children() && cell1->has_children())
615  {
616  bool grids_changed = false;
617 
618  Assert(cell2->n_children() == cell1->n_children(), ExcNotImplemented());
619  for (unsigned int c = 0; c < cell1->n_children(); ++c)
620  grids_changed |=
621  ::adapt_grid_cells<dim>(cell1->child(c), cell2->child(c));
622  return grids_changed;
623  }
624 
625 
626  if (!cell1->has_children() && !cell2->has_children())
627  // none of the two have children, so
628  // make sure that not one is flagged
629  // for refinement and the other for
630  // coarsening
631  {
632  if (cell1->refine_flag_set() && cell2->coarsen_flag_set())
633  {
634  cell2->clear_coarsen_flag();
635  return true;
636  }
637  else if (cell1->coarsen_flag_set() && cell2->refine_flag_set())
638  {
639  cell1->clear_coarsen_flag();
640  return true;
641  }
642 
643  return false;
644  }
645 
646 
647  if (cell1->has_children() && !cell2->has_children())
648  // cell1 has children, cell2 has not
649  // -> cell2 needs to be refined if any
650  // of cell1's children is flagged
651  // for refinement. None of them should
652  // be refined further, since then in the
653  // last round something must have gone
654  // wrong
655  //
656  // if cell2 was flagged for coarsening,
657  // we need to clear that flag in any
658  // case. The only exception would be
659  // if all children of cell1 were
660  // flagged for coarsening, but rules
661  // for coarsening are so complicated
662  // that we will not attempt to cover
663  // them. Rather accept one cell which
664  // is not coarsened...
665  {
666  bool changed_grid = false;
667  if (cell2->coarsen_flag_set())
668  {
669  cell2->clear_coarsen_flag();
670  changed_grid = true;
671  }
672 
673  if (!cell2->refine_flag_set())
674  for (unsigned int c = 0; c < cell1->n_children(); ++c)
675  if (cell1->child(c)->refine_flag_set() ||
676  cell1->child(c)->has_children())
677  {
678  cell2->set_refine_flag();
679  changed_grid = true;
680  break;
681  }
682  return changed_grid;
683  }
684 
685  if (!cell1->has_children() && cell2->has_children())
686  // same thing, other way round...
687  {
688  bool changed_grid = false;
689  if (cell1->coarsen_flag_set())
690  {
691  cell1->clear_coarsen_flag();
692  changed_grid = true;
693  }
694 
695  if (!cell1->refine_flag_set())
696  for (unsigned int c = 0; c < cell2->n_children(); ++c)
697  if (cell2->child(c)->refine_flag_set() ||
698  cell2->child(c)->has_children())
699  {
700  cell1->set_refine_flag();
701  changed_grid = true;
702  break;
703  }
704  return changed_grid;
705  }
706 
707  Assert(false, ExcInternalError());
708  return false;
709  }
710 
711 
712 
713  template <int dim>
714  bool
715  adapt_grids(Triangulation<dim> &tria1, Triangulation<dim> &tria2)
716  {
717  bool grids_changed = false;
718 
719  typename Triangulation<dim>::cell_iterator cell1 = tria1.begin(),
720  cell2 = tria2.begin();
721  typename Triangulation<dim>::cell_iterator endc;
722  endc = (tria1.n_levels() == 1 ?
723  typename Triangulation<dim>::cell_iterator(tria1.end()) :
724  tria1.begin(1));
725  for (; cell1 != endc; ++cell1, ++cell2)
726  grids_changed |= ::adapt_grid_cells<dim>(cell1, cell2);
727 
728  return grids_changed;
729  }
730 } // namespace
731 
732 
733 template <int dim>
734 void
735 TimeStepBase_Tria<dim>::refine_grid(const RefinementData refinement_data)
736 {
737  Vector<float> criteria;
739 
740  // copy the following two values since
741  // we may need modified values in the
742  // process of this function
743  double refinement_threshold = refinement_data.refinement_threshold,
744  coarsening_threshold = refinement_data.coarsening_threshold;
745 
746  // prepare an array where the criteria
747  // are stored in a sorted fashion
748  // we need this if cell number correction
749  // is switched on.
750  // the criteria are sorted in ascending
751  // order
752  // only fill it when needed
753  Vector<float> sorted_criteria;
754  // two pointers into this array denoting
755  // the position where the two thresholds
756  // are assumed
757  Vector<float>::const_iterator p_refinement_threshold = nullptr,
758  p_coarsening_threshold = nullptr;
759 
760 
761  // if we are to do some cell number
762  // correction steps, we have to find out
763  // which further cells (beyond
764  // refinement_threshold) to refine in case
765  // we need more cells, and which cells
766  // to not refine in case we need less cells
767  // (or even to coarsen, if necessary). to
768  // this end, we first define pointers into
769  // a sorted array of criteria pointing
770  // to the thresholds of refinement or
771  // coarsening; moving these pointers amounts
772  // to changing the threshold such that the
773  // number of cells flagged for refinement
774  // or coarsening would be changed by one
775  if ((timestep_no != 0) &&
776  (sweep_no >= refinement_flags.first_sweep_with_correction) &&
777  (refinement_flags.cell_number_correction_steps > 0))
778  {
779  sorted_criteria = criteria;
780  std::sort(sorted_criteria.begin(), sorted_criteria.end());
781  p_refinement_threshold =
782  Utilities::lower_bound(sorted_criteria.begin(),
783  sorted_criteria.end(),
784  static_cast<float>(refinement_threshold));
785  p_coarsening_threshold =
786  std::upper_bound(sorted_criteria.begin(),
787  sorted_criteria.end(),
788  static_cast<float>(coarsening_threshold));
789  }
790 
791 
792  // actually flag cells the first time
793  GridRefinement::refine(*tria, criteria, refinement_threshold);
794  GridRefinement::coarsen(*tria, criteria, coarsening_threshold);
795 
796  // store this number for the following
797  // since its computation is rather
798  // expensive and since it doesn't change
799  const unsigned int n_active_cells = tria->n_active_cells();
800 
801  // if not on first time level: try to
802  // adjust the number of resulting
803  // cells to those on the previous
804  // time level. Only do the cell number
805  // correction for higher sweeps and if
806  // there are sufficiently many cells
807  // already to avoid "grid stall" i.e.
808  // that the grid's evolution is hindered
809  // by the correction (this usually
810  // happens if there are very few cells,
811  // since then the number of cells touched
812  // by the correction step may exceed the
813  // number of cells which are flagged for
814  // refinement; in this case it often
815  // happens that the number of cells
816  // does not grow between sweeps, which
817  // clearly is not the wanted behaviour)
818  //
819  // however, if we do not do anything, we
820  // can get into trouble somewhen later.
821  // therefore, we also use the correction
822  // step for the first sweep or if the
823  // number of cells is between 100 and 300
824  // (unlike in the first version of the
825  // algorithm), but relax the conditions
826  // for the correction to allow deviations
827  // which are three times as high than
828  // allowed (sweep==1 || cell number<200)
829  // or twice as high (sweep==2 ||
830  // cell number<300). Also, since
831  // refinement never does any harm other
832  // than increased work, we allow for
833  // arbitrary growth of cell number if
834  // the estimated cell number is below
835  // 200.
836  //
837  // repeat this loop several times since
838  // the first estimate may not be totally
839  // correct
840  if ((timestep_no != 0) &&
841  (sweep_no >= refinement_flags.first_sweep_with_correction))
842  for (unsigned int loop = 0;
843  loop < refinement_flags.cell_number_correction_steps;
844  ++loop)
845  {
846  Triangulation<dim> *previous_tria =
847  dynamic_cast<const TimeStepBase_Tria<dim> *>(previous_timestep)->tria;
848 
849  // do one adaption step if desired
850  // (there are more coming below then
851  // also)
852  if (refinement_flags.adapt_grids)
853  ::adapt_grids<dim>(*previous_tria, *tria);
854 
855  // perform flagging of cells
856  // needed to regularize the
857  // triangulation
859  previous_tria->prepare_coarsening_and_refinement();
860 
861 
862  // now count the number of elements
863  // which will result on the previous
864  // grid after it will be refined. The
865  // number which will really result
866  // should be approximately that that we
867  // compute here, since we already
868  // performed most of the prepare*
869  // steps for the previous grid
870  //
871  // use a double value since for each
872  // four cells (in 2D) that we flagged
873  // for coarsening we result in one
874  // new. but since we loop over flagged
875  // cells, we have to subtract 3/4 of
876  // a cell for each flagged cell
878  Assert(!previous_tria->get_anisotropic_refinement_flag(),
880  double previous_cells = previous_tria->n_active_cells();
881  typename Triangulation<dim>::active_cell_iterator cell, endc;
882  cell = previous_tria->begin_active();
883  endc = previous_tria->end();
884  for (; cell != endc; ++cell)
885  if (cell->refine_flag_set())
886  previous_cells += (GeometryInfo<dim>::max_children_per_cell - 1);
887  else if (cell->coarsen_flag_set())
888  previous_cells -= static_cast<double>(
891 
892  // @p{previous_cells} now gives the
893  // number of cells which would result
894  // from the flags on the previous grid
895  // if we refined it now. However, some
896  // more flags will be set when we adapt
897  // the previous grid with this one
898  // after the flags have been set for
899  // this time level; on the other hand,
900  // we don't account for this, since the
901  // number of cells on this time level
902  // will be changed afterwards by the
903  // same way, when it is adapted to the
904  // next time level
905 
906  // now estimate the number of cells which
907  // will result on this level
908  double estimated_cells = n_active_cells;
909  cell = tria->begin_active();
910  endc = tria->end();
911  for (; cell != endc; ++cell)
912  if (cell->refine_flag_set())
913  estimated_cells += (GeometryInfo<dim>::max_children_per_cell - 1);
914  else if (cell->coarsen_flag_set())
915  estimated_cells -= static_cast<double>(
918 
919  // calculate the allowed delta in
920  // cell numbers; be more lenient
921  // if there are few cells
922  double delta_up = refinement_flags.cell_number_corridor_top,
923  delta_down = refinement_flags.cell_number_corridor_bottom;
924 
925  const std::vector<std::pair<unsigned int, double>> &relaxations =
926  (sweep_no >= refinement_flags.correction_relaxations.size() ?
927  refinement_flags.correction_relaxations.back() :
928  refinement_flags.correction_relaxations[sweep_no]);
929  for (unsigned int r = 0; r != relaxations.size(); ++r)
930  if (n_active_cells < relaxations[r].first)
931  {
932  delta_up *= relaxations[r].second;
933  delta_down *= relaxations[r].second;
934  break;
935  }
936 
937  // now, if the number of estimated
938  // cells exceeds the number of cells
939  // on the old time level by more than
940  // delta: cut the top threshold
941  //
942  // note that for each cell that
943  // we unflag we have to diminish the
944  // estimated number of cells by
945  // @p{children_per_cell}.
946  if (estimated_cells > previous_cells * (1. + delta_up))
947  {
948  // only limit the cell number
949  // if there will not be less
950  // than some number of cells
951  //
952  // also note that when using the
953  // dual estimator, the initial
954  // time level is not refined
955  // on its own, so we may not
956  // limit the number of the second
957  // time level on the basis of
958  // the initial one; since for
959  // the dual estimator, we
960  // mirror the refinement
961  // flags, the initial level
962  // will be passively refined
963  // later on.
964  if (estimated_cells > refinement_flags.min_cells_for_correction)
965  {
966  // number of cells by which the
967  // new grid is to be diminished
968  double delta_cells =
969  estimated_cells - previous_cells * (1. + delta_up);
970 
971  // if we need to reduce the
972  // number of cells, we need
973  // to raise the thresholds,
974  // i.e. move ahead in the
975  // sorted array, since this
976  // is sorted in ascending
977  // order. do so by removing
978  // cells tagged for refinement
979 
980  for (unsigned int i = 0; i < delta_cells;
982  if (p_refinement_threshold != sorted_criteria.end())
983  ++p_refinement_threshold;
984  else
985  break;
986  }
987  else
988  // too many cells, but we
989  // won't do anything about
990  // that
991  break;
992  }
993  else
994  // likewise: if the estimated number
995  // of cells is less than 90 per cent
996  // of those at the previous time level:
997  // raise threshold by refining
998  // additional cells. if we start to
999  // run into the area of cells
1000  // which are to be coarsened, we
1001  // raise the limit for these too
1002  if (estimated_cells < previous_cells * (1. - delta_down))
1003  {
1004  // number of cells by which the
1005  // new grid is to be enlarged
1006  double delta_cells =
1007  previous_cells * (1. - delta_down) - estimated_cells;
1008  // heuristics: usually, if we
1009  // add @p{delta_cells} to the
1010  // present state, we end up
1011  // with much more than only
1012  // (1-delta_down)*prev_cells
1013  // because of the effect of
1014  // regularization and because
1015  // of adaption to the
1016  // following grid. Therefore,
1017  // if we are not in the last
1018  // correction loop, we try not
1019  // to add as many cells as seem
1020  // necessary at first and hope
1021  // to get closer to the limit
1022  // this way. Only in the last
1023  // loop do we have to take the
1024  // full number to guarantee the
1025  // wanted result.
1026  //
1027  // The value 0.9 is taken from
1028  // practice, as the additional
1029  // number of cells introduced
1030  // by regularization is
1031  // approximately 10 per cent
1032  // of the flagged cells.
1033  if (loop != refinement_flags.cell_number_correction_steps - 1)
1034  delta_cells *= 0.9;
1035 
1036  // if more cells need to be
1037  // refined, we need to lower
1038  // the thresholds, i.e. to
1039  // move to the beginning
1040  // of sorted_criteria, which is
1041  // sorted in ascending order
1042  for (unsigned int i = 0; i < delta_cells;
1044  if (p_refinement_threshold != p_coarsening_threshold)
1045  --refinement_threshold;
1046  else if (p_coarsening_threshold != sorted_criteria.begin())
1047  --p_coarsening_threshold, --p_refinement_threshold;
1048  else
1049  break;
1050  }
1051  else
1052  // estimated cell number is ok,
1053  // stop correction steps
1054  break;
1055 
1056  if (p_refinement_threshold == sorted_criteria.end())
1057  {
1058  Assert(p_coarsening_threshold != p_refinement_threshold,
1059  ExcInternalError());
1060  --p_refinement_threshold;
1061  }
1062 
1063  coarsening_threshold = *p_coarsening_threshold;
1064  refinement_threshold = *p_refinement_threshold;
1065 
1066  if (coarsening_threshold >= refinement_threshold)
1067  coarsening_threshold = 0.999 * refinement_threshold;
1068 
1069  // now that we have re-adjusted
1070  // thresholds: clear all refine and
1071  // coarsening flags and do it all
1072  // over again
1073  cell = tria->begin_active();
1074  endc = tria->end();
1075  for (; cell != endc; ++cell)
1076  {
1077  cell->clear_refine_flag();
1078  cell->clear_coarsen_flag();
1079  }
1080 
1081 
1082  // flag cells finally
1083  GridRefinement::refine(*tria, criteria, refinement_threshold);
1084  GridRefinement::coarsen(*tria, criteria, coarsening_threshold);
1085  }
1086 
1087  // if step number is greater than
1088  // one: adapt this and the previous
1089  // grid to each other. Don't do so
1090  // for the initial grid because
1091  // it is always taken to be the first
1092  // grid and needs therefore no
1093  // treatment of its own.
1094  if ((timestep_no >= 1) && (refinement_flags.adapt_grids))
1095  {
1096  Triangulation<dim> *previous_tria =
1097  dynamic_cast<const TimeStepBase_Tria<dim> *>(previous_timestep)->tria;
1098  Assert(previous_tria != nullptr, ExcInternalError());
1099 
1100  // if we used the dual estimator, we
1101  // computed the error information on
1102  // a time slab, rather than on a level
1103  // of its own. we then mirror the
1104  // refinement flags we determined for
1105  // the present level to the previous
1106  // one
1107  //
1108  // do this mirroring only, if cell number
1109  // adjustment is on, since otherwise
1110  // strange things may happen
1111  if (refinement_flags.mirror_flags_to_previous_grid)
1112  {
1113  ::adapt_grids<dim>(*previous_tria, *tria);
1114 
1115  typename Triangulation<dim>::cell_iterator old_cell, new_cell, endc;
1116  old_cell = previous_tria->begin(0);
1117  new_cell = tria->begin(0);
1118  endc = tria->end(0);
1119  for (; new_cell != endc; ++new_cell, ++old_cell)
1120  ::mirror_refinement_flags<dim>(new_cell, old_cell);
1121  }
1122 
1124  previous_tria->prepare_coarsening_and_refinement();
1125 
1126  // adapt present and previous grids
1127  // to each other: flag additional
1128  // cells to avoid the previous grid
1129  // to have cells refined twice more
1130  // than the present one and vica versa.
1131  ::adapt_grids<dim>(*previous_tria, *tria);
1132  }
1133 }
1134 
1135 
1136 
1137 template <int dim>
1138 void
1140 {
1142 }
1143 
1144 
1145 
1146 template <int dim>
1147 std::size_t
1149 {
1150  return (TimeStepBase::memory_consumption() + sizeof(tria) +
1151  MemoryConsumption::memory_consumption(coarse_grid) + sizeof(flags) +
1152  sizeof(refinement_flags) +
1155 }
1156 
1157 
1158 
1159 template <int dim>
1161  : delete_and_rebuild_tria(false)
1162  , wakeup_level_to_build_grid(0)
1163  , sleep_level_to_delete_grid(0)
1164 {
1165  Assert(false, ExcInternalError());
1166 }
1167 
1168 
1169 
1170 template <int dim>
1172  const bool delete_and_rebuild_tria,
1173  const unsigned int wakeup_level_to_build_grid,
1174  const unsigned int sleep_level_to_delete_grid)
1175  : delete_and_rebuild_tria(delete_and_rebuild_tria)
1176  , wakeup_level_to_build_grid(wakeup_level_to_build_grid)
1177  , sleep_level_to_delete_grid(sleep_level_to_delete_grid)
1178 {
1179  // Assert (!delete_and_rebuild_tria || (wakeup_level_to_build_grid>=1),
1180  // ExcInvalidParameter(wakeup_level_to_build_grid));
1181  // Assert (!delete_and_rebuild_tria || (sleep_level_to_delete_grid>=1),
1182  // ExcInvalidParameter(sleep_level_to_delete_grid));
1183 }
1184 
1185 
1186 template <int dim>
1189  1, // one element, denoting the first and all subsequent sweeps
1190  std::vector<std::pair<unsigned int, double>>(1, // one element, denoting the
1191  // upper bound for the
1192  // following relaxation
1193  std::make_pair(0U, 0.)));
1194 
1195 
1196 template <int dim>
1198  const unsigned int max_refinement_level,
1199  const unsigned int first_sweep_with_correction,
1200  const unsigned int min_cells_for_correction,
1201  const double cell_number_corridor_top,
1202  const double cell_number_corridor_bottom,
1203  const CorrectionRelaxations &correction_relaxations,
1204  const unsigned int cell_number_correction_steps,
1205  const bool mirror_flags_to_previous_grid,
1206  const bool adapt_grids)
1207  : max_refinement_level(max_refinement_level)
1208  , first_sweep_with_correction(first_sweep_with_correction)
1209  , min_cells_for_correction(min_cells_for_correction)
1210  , cell_number_corridor_top(cell_number_corridor_top)
1211  , cell_number_corridor_bottom(cell_number_corridor_bottom)
1212  , correction_relaxations(correction_relaxations.size() != 0 ?
1213  correction_relaxations :
1214  default_correction_relaxations)
1215  , cell_number_correction_steps(cell_number_correction_steps)
1216  , mirror_flags_to_previous_grid(mirror_flags_to_previous_grid)
1217  , adapt_grids(adapt_grids)
1218 {
1219  Assert(cell_number_corridor_top >= 0,
1220  ExcInvalidValue(cell_number_corridor_top));
1221  Assert(cell_number_corridor_bottom >= 0,
1222  ExcInvalidValue(cell_number_corridor_bottom));
1223  Assert(cell_number_corridor_bottom <= 1,
1224  ExcInvalidValue(cell_number_corridor_bottom));
1225 }
1226 
1227 
1228 template <int dim>
1230  const double _refinement_threshold,
1231  const double _coarsening_threshold)
1232  : refinement_threshold(_refinement_threshold)
1233  ,
1234  // in some rare cases it may happen that
1235  // both thresholds are the same (e.g. if
1236  // there are many cells with the same
1237  // error indicator). That would mean that
1238  // all cells will be flagged for
1239  // refinement or coarsening, but some will
1240  // be flagged for both, namely those for
1241  // which the indicator equals the
1242  // thresholds. This is forbidden, however.
1243  //
1244  // In some rare cases with very few cells
1245  // we also could get integer round off
1246  // errors and get problems with
1247  // the top and bottom fractions.
1248  //
1249  // In these case we arbitrarily reduce the
1250  // bottom threshold by one permille below
1251  // the top threshold
1252  coarsening_threshold((_coarsening_threshold == _refinement_threshold ?
1253  _coarsening_threshold :
1254  0.999 * _coarsening_threshold))
1255 {
1258  // allow both thresholds to be zero,
1259  // since this is needed in case all indicators
1260  // are zero
1262  ((coarsening_threshold == 0) && (refinement_threshold == 0)),
1264 }
1265 
1266 
1267 
1268 /*-------------- Explicit Instantiations -------------------------------*/
1269 #include "time_dependent.inst"
1270 
1271 
1272 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:1079
virtual void sleep(const unsigned int)
virtual void start_sweep()
unsigned int n_active_cells() const
Definition: tria.cc:12602
static ::ExceptionBase & ExcPureFunctionCalled()
virtual void copy_triangulation(const Triangulation< dim, spacedim > &other_tria)
Definition: tria.cc:10479
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:1475
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:11940
const TimeSteppingData timestepping_data_dual
cell_iterator begin(const unsigned int level=0) const
Definition: tria.cc:11920
void apply_to_subranges(const RangeType &begin, const typename identity< RangeType >::type &end, const Function &f, const unsigned int grainsize)
Definition: parallel.h:444
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:12006
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:14021
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:1411
const double time
bool get_anisotropic_refinement_flag() const
Definition: tria.cc:10926
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:10880
void save_refine_flags(std::ostream &out) const
Definition: tria.cc:10812
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:1531
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()