Reference documentation for deal.II version GIT 01a9543f1b 2023-12-05 20:40: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\}}\)
mg_smoother.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_mg_smoother_h
17 #define dealii_mg_smoother_h
18 
19 
20 #include <deal.II/base/config.h>
21 
24 
27 
29 
30 #include <vector>
31 
33 
34 /*
35  * MGSmootherBase is defined in mg_base.h
36  */
37 
49 template <typename VectorType>
50 class MGSmoother : public MGSmootherBase<VectorType>
51 {
52 public:
56  MGSmoother(const unsigned int steps = 1,
57  const bool variable = false,
58  const bool symmetric = false,
59  const bool transpose = false);
60 
64  void
65  set_steps(const unsigned int);
66 
70  void
71  set_variable(const bool);
72 
76  void
77  set_symmetric(const bool);
78 
83  void
84  set_transpose(const bool);
85 
90  void
91  set_debug(const unsigned int level);
92 
93 protected:
101 
106  unsigned int steps;
107 
112  bool variable;
113 
118  bool symmetric;
119 
124  bool transpose;
125 
129  unsigned int debug;
130 };
131 
132 
139 template <typename VectorType>
140 class MGSmootherIdentity : public MGSmootherBase<VectorType>
141 {
142 public:
148  virtual void
149  smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;
150 
151  virtual void
152  clear();
153 };
154 
155 
156 namespace mg
157 {
184  template <typename RelaxationType, typename VectorType>
185  class SmootherRelaxation : public MGLevelObject<RelaxationType>,
186  public MGSmoother<VectorType>
187  {
188  public:
192  SmootherRelaxation(const unsigned int steps = 1,
193  const bool variable = false,
194  const bool symmetric = false,
195  const bool transpose = false);
196 
205  template <typename MatrixType2>
206  void
208  const typename RelaxationType::AdditionalData &additional_data =
209  typename RelaxationType::AdditionalData());
210 
218  template <typename MatrixType2, typename DataType>
219  void
221  const MGLevelObject<DataType> &additional_data);
222 
226  void
227  clear() override;
228 
232  virtual void
233  smooth(const unsigned int level,
234  VectorType &u,
235  const VectorType &rhs) const override;
236 
253  virtual void
254  apply(const unsigned int level,
255  VectorType &u,
256  const VectorType &rhs) const override;
257 
261  std::size_t
263  };
264 } // namespace mg
265 
295 template <typename MatrixType, typename RelaxationType, typename VectorType>
296 class MGSmootherRelaxation : public MGSmoother<VectorType>
297 {
298 public:
302  MGSmootherRelaxation(const unsigned int steps = 1,
303  const bool variable = false,
304  const bool symmetric = false,
305  const bool transpose = false);
306 
315  template <typename MatrixType2>
316  void
318  const typename RelaxationType::AdditionalData &additional_data =
319  typename RelaxationType::AdditionalData());
320 
329  template <typename MatrixType2, typename DataType>
330  void
332  const MGLevelObject<DataType> &additional_data);
333 
343  template <typename MatrixType2, typename DataType>
344  void
346  const DataType &additional_data,
347  const unsigned int block_row,
348  const unsigned int block_col);
349 
359  template <typename MatrixType2, typename DataType>
360  void
362  const MGLevelObject<DataType> &additional_data,
363  const unsigned int block_row,
364  const unsigned int block_col);
365 
369  void
370  clear();
371 
375  virtual void
376  smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;
377 
393  virtual void
394  apply(const unsigned int level, VectorType &u, const VectorType &rhs) const;
395 
400 
404  std::size_t
406 
407 
408 private:
413 };
414 
415 
416 
445 template <typename MatrixType, typename PreconditionerType, typename VectorType>
446 class MGSmootherPrecondition : public MGSmoother<VectorType>
447 {
448 public:
452  MGSmootherPrecondition(const unsigned int steps = 1,
453  const bool variable = false,
454  const bool symmetric = false,
455  const bool transpose = false);
456 
466  template <typename MatrixType2>
467  void
469  const typename PreconditionerType::AdditionalData &
470  additional_data = typename PreconditionerType::AdditionalData());
471 
479  template <typename MatrixType2>
480  void
482 
492  template <typename MatrixType2, typename DataType>
493  void
495  const MGLevelObject<DataType> &additional_data);
496 
507  template <typename MatrixType2, typename DataType>
508  void
510  const DataType &additional_data,
511  const unsigned int block_row,
512  const unsigned int block_col);
513 
524  template <typename MatrixType2, typename DataType>
525  void
527  const MGLevelObject<DataType> &additional_data,
528  const unsigned int block_row,
529  const unsigned int block_col);
530 
534  void
535  clear() override;
536 
540  virtual void
541  smooth(const unsigned int level,
542  VectorType &u,
543  const VectorType &rhs) const override;
544 
560  virtual void
561  apply(const unsigned int level,
562  VectorType &u,
563  const VectorType &rhs) const override;
564 
569 
573  std::size_t
575 
576 
577 private:
582 };
583 
586 /* ------------------------------- Inline functions --------------------------
587  */
588 
589 #ifndef DOXYGEN
590 
591 template <typename VectorType>
592 inline void
593 MGSmootherIdentity<VectorType>::smooth(const unsigned int,
594  VectorType &,
595  const VectorType &) const
596 {}
597 
598 template <typename VectorType>
599 inline void
601 {}
602 
603 //---------------------------------------------------------------------------
604 
605 template <typename VectorType>
606 inline MGSmoother<VectorType>::MGSmoother(const unsigned int steps,
607  const bool variable,
608  const bool symmetric,
609  const bool transpose)
610  : steps(steps)
611  , variable(variable)
614  , debug(0)
615 {}
616 
617 
618 template <typename VectorType>
619 inline void
620 MGSmoother<VectorType>::set_steps(const unsigned int s)
621 {
622  steps = s;
623 }
624 
625 
626 template <typename VectorType>
627 inline void
628 MGSmoother<VectorType>::set_debug(const unsigned int s)
629 {
630  debug = s;
631 }
632 
633 
634 template <typename VectorType>
635 inline void
637 {
638  variable = flag;
639 }
640 
641 
642 template <typename VectorType>
643 inline void
645 {
646  symmetric = flag;
647 }
648 
649 
650 template <typename VectorType>
651 inline void
653 {
654  transpose = flag;
655 }
656 
657 //----------------------------------------------------------------------//
658 
659 namespace mg
660 {
661  template <typename RelaxationType, typename VectorType>
663  const unsigned int steps,
664  const bool variable,
665  const bool symmetric,
666  const bool transpose)
667  : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
668  {}
669 
670 
671  template <typename RelaxationType, typename VectorType>
672  inline void
673  SmootherRelaxation<RelaxationType, VectorType>::clear()
674  {
676  }
677 
678 
679  template <typename RelaxationType, typename VectorType>
680  template <typename MatrixType2>
681  inline void
682  SmootherRelaxation<RelaxationType, VectorType>::initialize(
684  const typename RelaxationType::AdditionalData &data)
685  {
686  const unsigned int min = m.min_level();
687  const unsigned int max = m.max_level();
688 
689  this->resize(min, max);
690 
691  for (unsigned int i = min; i <= max; ++i)
692  (*this)[i].initialize(Utilities::get_underlying_value(m[i]), data);
693  }
694 
695 
696  template <typename RelaxationType, typename VectorType>
697  template <typename MatrixType2, typename DataType>
698  inline void
699  SmootherRelaxation<RelaxationType, VectorType>::initialize(
701  const MGLevelObject<DataType> &data)
702  {
703  const unsigned int min = std::max(m.min_level(), data.min_level());
704  const unsigned int max = std::min(m.max_level(), data.max_level());
705 
706  this->resize(min, max);
707 
708  for (unsigned int i = min; i <= max; ++i)
709  (*this)[i].initialize(Utilities::get_underlying_value(m[i]), data[i]);
710  }
711 
712 
713  template <typename RelaxationType, typename VectorType>
714  inline void
715  SmootherRelaxation<RelaxationType, VectorType>::smooth(
716  const unsigned int level,
717  VectorType &u,
718  const VectorType &rhs) const
719  {
720  unsigned int maxlevel = this->max_level();
721  unsigned int steps2 = this->steps;
722 
723  if (this->variable)
724  steps2 *= (1 << (maxlevel - level));
725 
726  bool T = this->transpose;
727  if (this->symmetric && (steps2 % 2 == 0))
728  T = false;
729  if (this->debug > 0)
730  deallog << 'S' << level << ' ';
731 
732  for (unsigned int i = 0; i < steps2; ++i)
733  {
734  if (T)
735  (*this)[level].Tstep(u, rhs);
736  else
737  (*this)[level].step(u, rhs);
738  if (this->symmetric)
739  T = !T;
740  }
741  }
742 
743 
744  template <typename RelaxationType, typename VectorType>
745  inline void
747  const unsigned int level,
748  VectorType &u,
749  const VectorType &rhs) const
750  {
751  unsigned int maxlevel = this->max_level();
752  unsigned int steps2 = this->steps;
753 
754  if (this->variable)
755  steps2 *= (1 << (maxlevel - level));
756 
757  bool T = this->transpose;
758  if (this->symmetric && (steps2 % 2 == 0))
759  T = false;
760  if (this->debug > 0)
761  deallog << 'S' << level << ' ';
762 
763  if (T)
764  (*this)[level].Tvmult(u, rhs);
765  else
766  (*this)[level].vmult(u, rhs);
767  if (this->symmetric)
768  T = !T;
769  for (unsigned int i = 1; i < steps2; ++i)
770  {
771  if (T)
772  (*this)[level].Tstep(u, rhs);
773  else
774  (*this)[level].step(u, rhs);
775  if (this->symmetric)
776  T = !T;
777  }
778  }
779 
780 
781  template <typename RelaxationType, typename VectorType>
782  inline std::size_t
784  {
785  return sizeof(*this) - sizeof(MGLevelObject<RelaxationType>) +
787  this->vector_memory.memory_consumption();
788  }
789 } // namespace mg
790 
791 
792 //----------------------------------------------------------------------//
793 
794 template <typename MatrixType, typename RelaxationType, typename VectorType>
796  MGSmootherRelaxation(const unsigned int steps,
797  const bool variable,
798  const bool symmetric,
799  const bool transpose)
800  : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
801 {}
802 
803 
804 
805 template <typename MatrixType, typename RelaxationType, typename VectorType>
806 inline void
808 {
809  smoothers.clear_elements();
810 
811  unsigned int i = matrices.min_level(), max_level = matrices.max_level();
812  for (; i <= max_level; ++i)
813  matrices[i] = LinearOperator<VectorType>();
814 }
815 
816 
817 template <typename MatrixType, typename RelaxationType, typename VectorType>
818 template <typename MatrixType2>
819 inline void
822  const typename RelaxationType::AdditionalData &data)
823 {
824  const unsigned int min = m.min_level();
825  const unsigned int max = m.max_level();
826 
827  matrices.resize(min, max);
828  smoothers.resize(min, max);
829 
830  for (unsigned int i = min; i <= max; ++i)
831  {
832  // Workaround: Unfortunately, not every "m[i]" object has a rich
833  // enough interface to populate reinit_(domain|range)_vector. Thus,
834  // apply an empty LinearOperator exemplar.
835  matrices[i] =
836  linear_operator<VectorType>(LinearOperator<VectorType>(),
838  smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data);
839  }
840 }
841 
842 template <typename MatrixType, typename RelaxationType, typename VectorType>
843 template <typename MatrixType2, typename DataType>
844 inline void
847  const MGLevelObject<DataType> &data)
848 {
849  const unsigned int min = m.min_level();
850  const unsigned int max = m.max_level();
851 
852  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
853  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
854 
855  matrices.resize(min, max);
856  smoothers.resize(min, max);
857 
858  for (unsigned int i = min; i <= max; ++i)
859  {
860  // Workaround: Unfortunately, not every "m[i]" object has a rich
861  // enough interface to populate reinit_(domain|range)_vector. Thus,
862  // apply an empty LinearOperator exemplar.
863  matrices[i] =
864  linear_operator<VectorType>(LinearOperator<VectorType>(),
866  smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data[i]);
867  }
868 }
869 
870 template <typename MatrixType, typename RelaxationType, typename VectorType>
871 template <typename MatrixType2, typename DataType>
872 inline void
875  const DataType &data,
876  const unsigned int row,
877  const unsigned int col)
878 {
879  const unsigned int min = m.min_level();
880  const unsigned int max = m.max_level();
881 
882  matrices.resize(min, max);
883  smoothers.resize(min, max);
884 
885  for (unsigned int i = min; i <= max; ++i)
886  {
887  // Workaround: Unfortunately, not every "m[i]" object has a rich
888  // enough interface to populate reinit_(domain|range)_vector. Thus,
889  // apply an empty LinearOperator exemplar.
890  matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
891  m[i].block(row, col));
892  smoothers[i].initialize(m[i].block(row, col), data);
893  }
894 }
895 
896 template <typename MatrixType, typename RelaxationType, typename VectorType>
897 template <typename MatrixType2, typename DataType>
898 inline void
901  const MGLevelObject<DataType> &data,
902  const unsigned int row,
903  const unsigned int col)
904 {
905  const unsigned int min = m.min_level();
906  const unsigned int max = m.max_level();
907 
908  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
909  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
910 
911  matrices.resize(min, max);
912  smoothers.resize(min, max);
913 
914  for (unsigned int i = min; i <= max; ++i)
915  {
916  // Workaround: Unfortunately, not every "m[i]" object has a rich
917  // enough interface to populate reinit_(domain|range)_vector. Thus,
918  // apply an empty LinearOperator exemplar.
919  matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
920  m[i].block(row, col));
921  smoothers[i].initialize(m[i].block(row, col), data[i]);
922  }
923 }
924 
925 
926 template <typename MatrixType, typename RelaxationType, typename VectorType>
927 inline void
929  const unsigned int level,
930  VectorType &u,
931  const VectorType &rhs) const
932 {
933  unsigned int maxlevel = smoothers.max_level();
934  unsigned int steps2 = this->steps;
935 
936  if (this->variable)
937  steps2 *= (1 << (maxlevel - level));
938 
939  bool T = this->transpose;
940  if (this->symmetric && (steps2 % 2 == 0))
941  T = false;
942  if (this->debug > 0)
943  deallog << 'S' << level << ' ';
944 
945  for (unsigned int i = 0; i < steps2; ++i)
946  {
947  if (T)
948  smoothers[level].Tstep(u, rhs);
949  else
950  smoothers[level].step(u, rhs);
951  if (this->symmetric)
952  T = !T;
953  }
954 }
955 
956 
957 template <typename MatrixType, typename RelaxationType, typename VectorType>
958 inline void
960  const unsigned int level,
961  VectorType &u,
962  const VectorType &rhs) const
963 {
964  unsigned int maxlevel = smoothers.max_level();
965  unsigned int steps2 = this->steps;
966 
967  if (this->variable)
968  steps2 *= (1 << (maxlevel - level));
969 
970  bool T = this->transpose;
971  if (this->symmetric && (steps2 % 2 == 0))
972  T = false;
973  if (this->debug > 0)
974  deallog << 'S' << level << ' ';
975 
976  if (T)
977  smoothers[level].Tvmult(u, rhs);
978  else
979  smoothers[level].vmult(u, rhs);
980  if (this->symmetric)
981  T = !T;
982  for (unsigned int i = 1; i < steps2; ++i)
983  {
984  if (T)
985  smoothers[level].Tstep(u, rhs);
986  else
987  smoothers[level].step(u, rhs);
988  if (this->symmetric)
989  T = !T;
990  }
991 }
992 
993 
994 
995 template <typename MatrixType, typename RelaxationType, typename VectorType>
996 inline std::size_t
998  memory_consumption() const
999 {
1000  return sizeof(*this) + matrices.memory_consumption() +
1001  smoothers.memory_consumption() +
1002  this->vector_memory.memory_consumption();
1003 }
1004 
1005 
1006 //----------------------------------------------------------------------//
1007 
1008 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1010  MGSmootherPrecondition(const unsigned int steps,
1011  const bool variable,
1012  const bool symmetric,
1013  const bool transpose)
1014  : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
1015 {}
1016 
1017 
1018 
1019 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1020 inline void
1022 {
1023  smoothers.clear_elements();
1024 
1025  unsigned int i = matrices.min_level(), max_level = matrices.max_level();
1026  for (; i <= max_level; ++i)
1027  matrices[i] = LinearOperator<VectorType>();
1028 }
1029 
1030 
1031 
1032 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1033 template <typename MatrixType2>
1034 inline void
1036  const MGLevelObject<MatrixType2> &m,
1037  const typename PreconditionerType::AdditionalData &data)
1038 {
1039  const unsigned int min = m.min_level();
1040  const unsigned int max = m.max_level();
1041 
1042  matrices.resize(min, max);
1043  smoothers.resize(min, max);
1044 
1045  for (unsigned int i = min; i <= max; ++i)
1046  {
1047  // Workaround: Unfortunately, not every "m[i]" object has a rich
1048  // enough interface to populate reinit_(domain|range)_vector. Thus,
1049  // apply an empty LinearOperator exemplar.
1050  matrices[i] =
1051  linear_operator<VectorType>(LinearOperator<VectorType>(),
1053  smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data);
1054  }
1055 }
1056 
1057 
1058 
1059 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1060 template <typename MatrixType2>
1061 inline void
1064 {
1065  const unsigned int min = m.min_level();
1066  const unsigned int max = m.max_level();
1067 
1068  matrices.resize(min, max);
1069  smoothers.resize(min, max);
1070 
1071  for (unsigned int i = min; i <= max; ++i)
1072  {
1073  // Workaround: Unfortunately, not every "m[i]" object has a rich
1074  // enough interface to populate reinit_(domain|range)_vector. Thus,
1075  // apply an empty LinearOperator exemplar.
1076  matrices[i] =
1077  linear_operator<VectorType>(LinearOperator<VectorType>(),
1079  }
1080 }
1081 
1082 
1083 
1084 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1085 template <typename MatrixType2, typename DataType>
1086 inline void
1088  const MGLevelObject<MatrixType2> &m,
1089  const MGLevelObject<DataType> &data)
1090 {
1091  const unsigned int min = m.min_level();
1092  const unsigned int max = m.max_level();
1093 
1094  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1095  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1096 
1097  matrices.resize(min, max);
1098  smoothers.resize(min, max);
1099 
1100  for (unsigned int i = min; i <= max; ++i)
1101  {
1102  // Workaround: Unfortunately, not every "m[i]" object has a rich
1103  // enough interface to populate reinit_(domain|range)_vector. Thus,
1104  // apply an empty LinearOperator exemplar.
1105  matrices[i] =
1106  linear_operator<VectorType>(LinearOperator<VectorType>(),
1108  smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data[i]);
1109  }
1110 }
1111 
1112 
1113 
1114 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1115 template <typename MatrixType2, typename DataType>
1116 inline void
1118  const MGLevelObject<MatrixType2> &m,
1119  const DataType &data,
1120  const unsigned int row,
1121  const unsigned int col)
1122 {
1123  const unsigned int min = m.min_level();
1124  const unsigned int max = m.max_level();
1125 
1126  matrices.resize(min, max);
1127  smoothers.resize(min, max);
1128 
1129  for (unsigned int i = min; i <= max; ++i)
1130  {
1131  matrices[i] = &(m[i].block(row, col));
1132  smoothers[i].initialize(m[i].block(row, col), data);
1133  }
1134 }
1135 
1136 
1137 
1138 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1139 template <typename MatrixType2, typename DataType>
1140 inline void
1142  const MGLevelObject<MatrixType2> &m,
1143  const MGLevelObject<DataType> &data,
1144  const unsigned int row,
1145  const unsigned int col)
1146 {
1147  const unsigned int min = m.min_level();
1148  const unsigned int max = m.max_level();
1149 
1150  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1151  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1152 
1153  matrices.resize(min, max);
1154  smoothers.resize(min, max);
1155 
1156  for (unsigned int i = min; i <= max; ++i)
1157  {
1158  matrices[i] = &(m[i].block(row, col));
1159  smoothers[i].initialize(m[i].block(row, col), data[i]);
1160  }
1161 }
1162 
1163 
1164 
1165 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1166 inline void
1168  const unsigned int level,
1169  VectorType &u,
1170  const VectorType &rhs) const
1171 {
1172  unsigned int maxlevel = matrices.max_level();
1173  unsigned int steps2 = this->steps;
1174 
1175  if (this->variable)
1176  steps2 *= (1 << (maxlevel - level));
1177 
1178  typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1179  typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1180 
1181  r->reinit(u, true);
1182  d->reinit(u, true);
1183 
1184  bool T = this->transpose;
1185  if (this->symmetric && (steps2 % 2 == 0))
1186  T = false;
1187  if (this->debug > 0)
1188  deallog << 'S' << level << ' ';
1189 
1190  for (unsigned int i = 0; i < steps2; ++i)
1191  {
1192  if (T)
1193  {
1194  if (this->debug > 0)
1195  deallog << 'T';
1196  matrices[level].Tvmult(*r, u);
1197  r->sadd(-1., 1., rhs);
1198  if (this->debug > 2)
1199  deallog << ' ' << r->l2_norm() << ' ';
1200  smoothers[level].Tvmult(*d, *r);
1201  if (this->debug > 1)
1202  deallog << ' ' << d->l2_norm() << ' ';
1203  }
1204  else
1205  {
1206  if (this->debug > 0)
1207  deallog << 'N';
1208  matrices[level].vmult(*r, u);
1209  r->sadd(-1., rhs);
1210  if (this->debug > 2)
1211  deallog << ' ' << r->l2_norm() << ' ';
1212  smoothers[level].vmult(*d, *r);
1213  if (this->debug > 1)
1214  deallog << ' ' << d->l2_norm() << ' ';
1215  }
1216  u += *d;
1217  if (this->symmetric)
1218  T = !T;
1219  }
1220  if (this->debug > 0)
1221  deallog << std::endl;
1222 }
1223 
1224 
1225 
1226 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1227 inline void
1229  const unsigned int level,
1230  VectorType &u,
1231  const VectorType &rhs) const
1232 {
1233  unsigned int maxlevel = matrices.max_level();
1234  unsigned int steps2 = this->steps;
1235 
1236  if (this->variable)
1237  steps2 *= (1 << (maxlevel - level));
1238 
1239  bool T = this->transpose;
1240  if (this->symmetric && (steps2 % 2 == 0))
1241  T = false;
1242  if (this->debug > 0)
1243  deallog << 'S' << level << ' ';
1244 
1245  // first step where we overwrite the result
1246  if (this->debug > 2)
1247  deallog << ' ' << rhs.l2_norm() << ' ';
1248  if (this->debug > 0)
1249  deallog << (T ? 'T' : 'N');
1250  if (T)
1251  smoothers[level].Tvmult(u, rhs);
1252  else
1253  smoothers[level].vmult(u, rhs);
1254  if (this->debug > 1)
1255  deallog << ' ' << u.l2_norm() << ' ';
1256  if (this->symmetric)
1257  T = !T;
1258 
1259  typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1260  typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1261 
1262  if (steps2 > 1)
1263  {
1264  r->reinit(u, true);
1265  d->reinit(u, true);
1266  }
1267 
1268  for (unsigned int i = 1; i < steps2; ++i)
1269  {
1270  if (T)
1271  {
1272  if (this->debug > 0)
1273  deallog << 'T';
1274  matrices[level].Tvmult(*r, u);
1275  r->sadd(-1., 1., rhs);
1276  if (this->debug > 2)
1277  deallog << ' ' << r->l2_norm() << ' ';
1278  smoothers[level].Tvmult(*d, *r);
1279  if (this->debug > 1)
1280  deallog << ' ' << d->l2_norm() << ' ';
1281  }
1282  else
1283  {
1284  if (this->debug > 0)
1285  deallog << 'N';
1286  matrices[level].vmult(*r, u);
1287  r->sadd(-1., rhs);
1288  if (this->debug > 2)
1289  deallog << ' ' << r->l2_norm() << ' ';
1290  smoothers[level].vmult(*d, *r);
1291  if (this->debug > 1)
1292  deallog << ' ' << d->l2_norm() << ' ';
1293  }
1294  u += *d;
1295  if (this->symmetric)
1296  T = !T;
1297  }
1298  if (this->debug > 0)
1299  deallog << std::endl;
1300 }
1301 
1302 
1303 
1304 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1305 inline std::size_t
1307  memory_consumption() const
1308 {
1309  return sizeof(*this) + matrices.memory_consumption() +
1310  smoothers.memory_consumption() +
1311  this->vector_memory.memory_consumption();
1312 }
1313 
1314 
1315 #endif // DOXYGEN
1316 
1318 
1319 #endif
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
std::size_t memory_consumption() const
unsigned int max_level() const
unsigned int min_level() const
virtual void clear()
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const
std::size_t memory_consumption() const
MGLevelObject< PreconditionerType > smoothers
Definition: mg_smoother.h:568
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DataType > &additional_data)
void initialize_matrices(const MGLevelObject< MatrixType2 > &matrices)
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const override
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename PreconditionerType::AdditionalData &additional_data=typename PreconditionerType::AdditionalData())
virtual void apply(const unsigned int level, VectorType &u, const VectorType &rhs) const override
void clear() override
void initialize(const MGLevelObject< MatrixType2 > &matrices, const DataType &additional_data, const unsigned int block_row, const unsigned int block_col)
MGLevelObject< LinearOperator< VectorType > > matrices
Definition: mg_smoother.h:581
MGSmootherPrecondition(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DataType > &additional_data, const unsigned int block_row, const unsigned int block_col)
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const
MGLevelObject< LinearOperator< VectorType > > matrices
Definition: mg_smoother.h:412
void initialize(const MGLevelObject< MatrixType2 > &matrices, const DataType &additional_data, const unsigned int block_row, const unsigned int block_col)
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DataType > &additional_data)
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DataType > &additional_data, const unsigned int block_row, const unsigned int block_col)
std::size_t memory_consumption() const
virtual void apply(const unsigned int level, VectorType &u, const VectorType &rhs) const
MGLevelObject< RelaxationType > smoothers
Definition: mg_smoother.h:399
MGSmootherRelaxation(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename RelaxationType::AdditionalData &additional_data=typename RelaxationType::AdditionalData())
bool transpose
Definition: mg_smoother.h:124
GrowingVectorMemory< VectorType > vector_memory
Definition: mg_smoother.h:100
MGSmoother(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
bool variable
Definition: mg_smoother.h:112
void set_debug(const unsigned int level)
void set_steps(const unsigned int)
unsigned int debug
Definition: mg_smoother.h:129
bool symmetric
Definition: mg_smoother.h:118
void set_symmetric(const bool)
void set_transpose(const bool)
unsigned int steps
Definition: mg_smoother.h:106
void set_variable(const bool)
VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
std::size_t memory_consumption() const
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DataType > &additional_data)
SmootherRelaxation(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const override
virtual void apply(const unsigned int level, VectorType &u, const VectorType &rhs) const override
void clear() override
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename RelaxationType::AdditionalData &additional_data=typename RelaxationType::AdditionalData())
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
unsigned int level
Definition: grid_out.cc:4617
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1631
LogStream deallog
Definition: logstream.cc:37
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
@ symmetric
Matrix is symmetric.
static const char T
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T & get_underlying_value(T &p)
Definition: utilities.h:1601
Definition: mg.h:82