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