Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
time_dependent.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1999 - 2023 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
19
21#include <deal.II/grid/tria.h>
24
25#include <deal.II/lac/vector.h>
26
28
29#include <algorithm>
30#include <functional>
31#include <numeric>
32
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)
58 }
59 catch (...)
60 {}
61}
62
63
64void
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(),
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
126void
128{
129 insert_timestep(nullptr, new_timestep);
130}
131
132
133void
134TimeDependent::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
168void
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
179void
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
190void
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
202void
203TimeDependent::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
226void
228{
230 0U,
231 timesteps.size(),
232 [this](const unsigned int begin, const unsigned int end) {
233 this->end_sweep(begin, end);
234 },
235 1);
236}
237
238
239
240void
241TimeDependent::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
249std::size_t
251{
252 std::size_t mem =
257 for (const auto &timestep : timesteps)
259
260 return mem;
261}
262
263
264
265/* --------------------------------------------------------------------- */
266
267
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
279void
280TimeStepBase::wake_up(const unsigned int)
281{}
282
283
284
285void
286TimeStepBase::sleep(const unsigned)
287{}
288
289
290
291void
293{}
294
295
296
297void
299{}
300
301
302
303void
305{
307}
308
309
310
311void
313{
315}
316
317
318
319void
321{
323}
324
325
326
327void
329{
331}
332
333
334
335void
337{
339}
340
341
342
343double
345{
346 return time;
347}
348
349
350
351unsigned int
353{
354 return timestep_no;
355}
356
357
358
359double
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
370double
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
381void
383{
384 previous_timestep = previous;
385}
386
387
388
389void
391{
392 next_timestep = next;
393}
394
395
396
397void
398TimeStepBase::set_timestep_no(const unsigned int step_no)
399{
400 timestep_no = step_no;
401}
402
403
404
405void
406TimeStepBase::set_sweep_no(const unsigned int sweep)
407{
408 sweep_no = sweep;
409}
410
411
412
413std::size_t
415{
416 // only simple data types
417 return sizeof(*this);
418}
419
420
421
422template <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{
431}
432
433
434
435#ifndef DOXYGEN
436template <int dim>
438 const double time,
439 const Triangulation<dim> &coarse_grid,
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)
446 , refinement_flags(refinement_flags)
447{}
448#endif
449
450
451
452template <int dim>
454{
455 if (!flags.delete_and_rebuild_tria)
456 {
458 tria = nullptr;
459 delete t;
460 }
461 else
462 AssertNothrow(tria == nullptr, ExcInternalError());
463
464 coarse_grid = nullptr;
465}
466
467
468
469template <int dim>
470void
471TimeStepBase_Tria<dim>::wake_up(const unsigned int wakeup_level)
472{
473 TimeStepBase::wake_up(wakeup_level);
474
475 if (wakeup_level == flags.wakeup_level_to_build_grid)
476 if (flags.delete_and_rebuild_tria || !tria)
477 restore_grid();
478}
479
480
481
482template <int dim>
483void
484TimeStepBase_Tria<dim>::sleep(const unsigned int sleep_level)
485{
486 if (sleep_level == flags.sleep_level_to_delete_grid)
487 {
488 Assert(tria != nullptr, ExcInternalError());
489
490 if (flags.delete_and_rebuild_tria)
491 {
493 tria = nullptr;
494 delete t;
495 }
496 }
497
498 TimeStepBase::sleep(sleep_level);
499}
500
501
502
503template <int dim>
504void
506{
507 // for any of the non-initial grids
508 // store the refinement flags
509 refine_flags.emplace_back();
510 coarsen_flags.emplace_back();
511 tria->save_refine_flags(refine_flags.back());
512 tria->save_coarsen_flags(coarsen_flags.back());
513}
514
515
516
517template <int dim>
518void
520{
521 Assert(tria == nullptr, ExcGridNotDeleted());
522 Assert(refine_flags.size() == coarsen_flags.size(), ExcInternalError());
523
524 // create a virgin triangulation and
525 // set it to a copy of the coarse grid
526 tria = new Triangulation<dim>();
527 tria->copy_triangulation(*coarse_grid);
528
529 // for each of the previous refinement
530 // sweeps
531 for (unsigned int previous_sweep = 0; previous_sweep < refine_flags.size();
532 ++previous_sweep)
533 {
534 // get flags
535 tria->load_refine_flags(refine_flags[previous_sweep]);
536 tria->load_coarsen_flags(coarsen_flags[previous_sweep]);
537
538 // limit refinement depth if the user
539 // desired so
540 // if (flags.max_refinement_level != 0)
541 // {
542 // typename Triangulation<dim>::active_cell_iterator cell, endc;
543 // for (const auto &cell : tria->active_cell_iterators())
544 // if (static_cast<unsigned int>(cell->level()) >=
545 // flags.max_refinement_level)
546 // cell->clear_refine_flag();
547 // }
548
550 }
551}
552
553
554
555// have a few helper functions
556namespace
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->is_active())
587 {
588 if (new_cell->refine_flag_set() && old_cell->is_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();
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
735template <int dim>
736void
738{
739 Vector<float> criteria;
740 get_tria_refinement_criteria(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 behavior)
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();
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 (const auto &relaxation : relaxations)
932 if (n_active_cells < relaxation.first)
933 {
934 delta_up *= relaxation.second;
935 delta_down *= relaxation.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,
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
1139template <int dim>
1140void
1142{
1143 next_action = grid_refinement;
1144}
1145
1146
1147
1148template <int dim>
1149std::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
1161template <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
1172template <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
1188template <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
1198template <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{
1227}
1228
1229
1230template <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
void insert_timestep(const TimeStepBase *position, TimeStepBase *new_timestep)
virtual void start_sweep(const unsigned int sweep_no)
std::size_t memory_consumption() const
const TimeSteppingData timestepping_data_primal
virtual void end_sweep()
std::vector< SmartPointer< TimeStepBase, TimeDependent > > timesteps
void solve_dual_problem()
virtual ~TimeDependent()
void do_loop(InitFunctionObject init_function, LoopFunctionObject loop_function, const TimeSteppingData &timestepping_data, const Direction direction)
const TimeSteppingData timestepping_data_dual
void solve_primal_problem()
unsigned int sweep_no
TimeDependent(const TimeSteppingData &data_primal, const TimeSteppingData &data_dual, const TimeSteppingData &data_postprocess)
void delete_timestep(const unsigned int position)
void add_timestep(TimeStepBase *new_timestep)
const TimeSteppingData timestepping_data_postprocess
virtual void wake_up(const unsigned int wakeup_level) override
virtual void sleep(const unsigned int) override
virtual std::size_t memory_consumption() const override
void refine_grid(const RefinementData data)
virtual void init_for_refinement()
typename TimeStepBase_Tria_Flags::RefinementData< dim > RefinementData
virtual ~TimeStepBase_Tria() override
virtual std::size_t memory_consumption() const
double get_forward_timestep() const
TimeStepBase(const double time)
virtual void wake_up(const unsigned int)
void set_timestep_no(const unsigned int step_no)
void set_previous_timestep(const TimeStepBase *previous)
const TimeStepBase * previous_timestep
virtual void postprocess_timestep()
virtual void sleep(const unsigned int)
unsigned int timestep_no
void set_next_timestep(const TimeStepBase *next)
unsigned int get_timestep_no() const
unsigned int next_action
double get_backward_timestep() const
void set_sweep_no(const unsigned int sweep_no)
const double time
unsigned int sweep_no
virtual void end_sweep()
virtual void solve_primal_problem()=0
virtual void init_for_postprocessing()
virtual void start_sweep()
virtual void init_for_primal_problem()
double get_time() const
virtual void solve_dual_problem()
virtual void init_for_dual_problem()
const TimeStepBase * next_timestep
bool get_anisotropic_refinement_flag() const
virtual void copy_triangulation(const Triangulation< dim, spacedim > &other_tria)
cell_iterator begin(const unsigned int level=0) const
void load_refine_flags(std::istream &in)
unsigned int n_active_cells() const
unsigned int n_levels() const
cell_iterator end() const
void load_coarsen_flags(std::istream &out)
virtual void execute_coarsening_and_refinement()
virtual bool prepare_coarsening_and_refinement()
void save_refine_flags(std::ostream &out) const
void save_coarsen_flags(std::ostream &out) const
active_cell_iterator begin_active(const unsigned int level=0) const
const value_type * const_iterator
Definition vector.h:132
iterator end()
iterator begin()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcInvalidValue(double arg1)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcPureFunctionCalled()
#define AssertNothrow(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInvalidValue(double arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcInvalidPosition()
TriaIterator< CellAccessor< dim, spacedim > > cell_iterator
Definition tria.h:1370
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void coarsen(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold)
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
VectorType::value_type * begin(VectorType &V)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition utilities.h:1016
void apply_to_subranges(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, const Function &f, const unsigned int grainsize)
Definition parallel.h:435
TimeSteppingData(const unsigned int look_ahead, const unsigned int look_back)
RefinementData(const double refinement_threshold, const double coarsening_threshold=0)
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)
std::vector< std::vector< std::pair< unsigned int, double > > > CorrectionRelaxations
static CorrectionRelaxations default_correction_relaxations
const ::Triangulation< dim, spacedim > & tria