deal.II version GIT relicensing-2248-g6a8f632045 2024-12-12 07:30:00+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<ObserverPointer<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 // ObserverPointer
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) ?
163 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
294
295
296
297void
300
301
302
303void
308
309
310
311void
316
317
318
319void
324
325
326
327void
332
333
334
335void
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
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 {
457 Triangulation<dim> *t = tria;
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 {
492 Triangulation<dim> *t = tria;
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 an empty 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
549 tria->execute_coarsening_and_refinement();
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
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
860 tria->prepare_coarsening_and_refinement();
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
879 Assert(!tria->get_anisotropic_refinement_flag(), ExcNotImplemented());
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
1125 tria->prepare_coarsening_and_refinement();
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{
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
std::vector< ObserverPointer< TimeStepBase, TimeDependent > > timesteps
const TimeSteppingData timestepping_data_primal
virtual void end_sweep()
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
cell_iterator begin(const unsigned int level=0) const
unsigned int n_active_cells() const
unsigned int n_levels() const
cell_iterator end() const
virtual bool prepare_coarsening_and_refinement()
active_cell_iterator begin_active(const unsigned int level=0) const
const value_type * const_iterator
Definition vector.h:141
iterator end()
iterator begin()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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:1556
#define DEAL_II_ASSERT_UNREACHABLE()
std::size_t size
Definition mpi.cc:734
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)
VectorType::value_type * begin(VectorType &V)
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
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