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