Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-52.h
Go to the documentation of this file.
1
714,
715 *   const double tau,
716 *   const Vector<double> &y)
717 *   {
718 *   SparseDirectUMFPACK inverse_mass_minus_tau_Jacobian;
719 *  
720 *   mass_minus_tau_Jacobian.copy_from(mass_matrix);
721 *   mass_minus_tau_Jacobian.add(-tau, system_matrix);
722 *  
723 *   inverse_mass_minus_tau_Jacobian.initialize(mass_minus_tau_Jacobian);
724 *  
725 *   Vector<double> tmp(dof_handler.n_dofs());
726 *   mass_matrix.vmult(tmp, y);
727 *  
728 *   Vector<double> result(y);
729 *   inverse_mass_minus_tau_Jacobian.vmult(result, tmp);
730 *  
731 *   return result;
732 *   }
733 *  
734 *  
735 *  
736 * @endcode
737 *
738 *
739 * <a name="step_52-codeDiffusionoutput_resultscode"></a>
740 * <h4><code>Diffusion::output_results</code></h4>
741 *
742
743 *
744 * The following function then outputs the solution in vtu files indexed by
745 * the number of the time step and the name of the time stepping method. Of
746 * course, the (exact) result should really be the same for all time
747 * stepping method, but the output here at least allows us to compare them.
748 *
749 * @code
750 *   void Diffusion::output_results(const double time,
751 *   const unsigned int time_step,
752 *   TimeStepping::runge_kutta_method method) const
753 *   {
754 *   std::string method_name;
755 *  
756 *   switch (method)
757 *   {
759 *   {
760 *   method_name = "forward_euler";
761 *   break;
762 *   }
764 *   {
765 *   method_name = "rk3";
766 *   break;
767 *   }
769 *   {
770 *   method_name = "rk4";
771 *   break;
772 *   }
774 *   {
775 *   method_name = "backward_euler";
776 *   break;
777 *   }
779 *   {
780 *   method_name = "implicit_midpoint";
781 *   break;
782 *   }
784 *   {
785 *   method_name = "sdirk";
786 *   break;
787 *   }
788 *   case TimeStepping::HEUN_EULER:
789 *   {
790 *   method_name = "heun_euler";
791 *   break;
792 *   }
794 *   {
795 *   method_name = "bogacki_shampine";
796 *   break;
797 *   }
798 *   case TimeStepping::DOPRI:
799 *   {
800 *   method_name = "dopri";
801 *   break;
802 *   }
803 *   case TimeStepping::FEHLBERG:
804 *   {
805 *   method_name = "fehlberg";
806 *   break;
807 *   }
808 *   case TimeStepping::CASH_KARP:
809 *   {
810 *   method_name = "cash_karp";
811 *   break;
812 *   }
813 *   default:
814 *   {
815 *   break;
816 *   }
817 *   }
818 *  
819 *   DataOut<2> data_out;
820 *  
821 *   data_out.attach_dof_handler(dof_handler);
822 *   data_out.add_data_vector(solution, "solution");
823 *  
824 *   data_out.build_patches();
825 *  
826 *   data_out.set_flags(DataOutBase::VtkFlags(time, time_step));
827 *  
828 *   const std::string filename = "solution_" + method_name + "-" +
829 *   Utilities::int_to_string(time_step, 3) +
830 *   ".vtu";
831 *   std::ofstream output(filename);
832 *   data_out.write_vtu(output);
833 *  
834 *   static std::vector<std::pair<double, std::string>> times_and_names;
835 *  
836 *   static std::string method_name_prev = "";
837 *   static std::string pvd_filename;
838 *   if (method_name_prev != method_name)
839 *   {
840 *   times_and_names.clear();
841 *   method_name_prev = method_name;
842 *   pvd_filename = "solution_" + method_name + ".pvd";
843 *   }
844 *   times_and_names.emplace_back(time, filename);
845 *   std::ofstream pvd_output(pvd_filename);
846 *   DataOutBase::write_pvd_record(pvd_output, times_and_names);
847 *   }
848 *  
849 *  
850 * @endcode
851 *
852 *
853 * <a name="step_52-codeDiffusionexplicit_methodcode"></a>
854 * <h4><code>Diffusion::explicit_method</code></h4>
855 *
856
857 *
858 * This function is the driver for all the explicit methods. At the
859 * top it initializes the time stepping and the solution (by setting
860 * it to zero and then ensuring that boundary value and hanging node
861 * constraints are respected; of course, with the mesh we use here,
862 * hanging node constraints are not in fact an issue). It then calls
863 * <code>evolve_one_time_step</code> which performs one time step.
864 * Time is stored and incremented through a DiscreteTime object.
865 *
866
867 *
868 * For explicit methods, <code>evolve_one_time_step</code> needs to
869 * evaluate @f$M^{-1}(f(t,y))@f$, i.e, it needs
870 * <code>evaluate_diffusion</code>. Because
871 * <code>evaluate_diffusion</code> is a member function, it needs to
872 * be bound to <code>this</code>. After each evolution step, we
873 * again apply the correct boundary values and hanging node
874 * constraints.
875 *
876
877 *
878 * Finally, the solution is output
879 * every 10 time steps.
880 *
881 * @code
882 *   void Diffusion::explicit_method(const TimeStepping::runge_kutta_method method,
883 *   const unsigned int n_time_steps,
884 *   const double initial_time,
885 *   const double final_time)
886 *   {
887 *   const double time_step =
888 *   (final_time - initial_time) / static_cast<double>(n_time_steps);
889 *  
890 *   solution = 0.;
891 *   constraint_matrix.distribute(solution);
892 *  
893 *   TimeStepping::ExplicitRungeKutta<Vector<double>> explicit_runge_kutta(
894 *   method);
895 *   output_results(initial_time, 0, method);
896 *   DiscreteTime time(initial_time, final_time, time_step);
897 *   while (time.is_at_end() == false)
898 *   {
899 *   explicit_runge_kutta.evolve_one_time_step(
900 *   [this](const double time, const Vector<double> &y) {
901 *   return this->evaluate_diffusion(time, y);
902 *   },
903 *   time.get_current_time(),
904 *   time.get_next_step_size(),
905 *   solution);
906 *   time.advance_time();
907 *  
908 *   constraint_matrix.distribute(solution);
909 *  
910 *   if (time.get_step_number() % 10 == 0)
911 *   output_results(time.get_current_time(),
912 *   time.get_step_number(),
913 *   method);
914 *   }
915 *   }
916 *  
917 *  
918 *  
919 * @endcode
920 *
921 *
922 * <a name="step_52-codeDiffusionimplicit_methodcode"></a>
923 * <h4><code>Diffusion::implicit_method</code></h4>
924 * This function is equivalent to <code>explicit_method</code> but for
925 * implicit methods. When using implicit methods, we need to evaluate
926 * @f$M^{-1}(f(t,y))@f$ and @f$\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial
927 * y}\right)^{-1}@f$ for which we use the two member functions previously
928 * introduced.
929 *
930 * @code
931 *   void Diffusion::implicit_method(const TimeStepping::runge_kutta_method method,
932 *   const unsigned int n_time_steps,
933 *   const double initial_time,
934 *   const double final_time)
935 *   {
936 *   const double time_step =
937 *   (final_time - initial_time) / static_cast<double>(n_time_steps);
938 *  
939 *   solution = 0.;
940 *   constraint_matrix.distribute(solution);
941 *  
942 *   TimeStepping::ImplicitRungeKutta<Vector<double>> implicit_runge_kutta(
943 *   method);
944 *   output_results(initial_time, 0, method);
945 *   DiscreteTime time(initial_time, final_time, time_step);
946 *   while (time.is_at_end() == false)
947 *   {
948 *   implicit_runge_kutta.evolve_one_time_step(
949 *   [this](const double time, const Vector<double> &y) {
950 *   return this->evaluate_diffusion(time, y);
951 *   },
952 *   [this](const double time, const double tau, const Vector<double> &y) {
953 *   return this->id_minus_tau_J_inverse(time, tau, y);
954 *   },
955 *   time.get_current_time(),
956 *   time.get_next_step_size(),
957 *   solution);
958 *   time.advance_time();
959 *  
960 *   constraint_matrix.distribute(solution);
961 *  
962 *   if (time.get_step_number() % 10 == 0)
963 *   output_results(time.get_current_time(),
964 *   time.get_step_number(),
965 *   method);
966 *   }
967 *   }
968 *  
969 *  
970 *  
971 * @endcode
972 *
973 *
974 * <a name="step_52-codeDiffusionembedded_explicit_methodcode"></a>
975 * <h4><code>Diffusion::embedded_explicit_method</code></h4>
976 * This function is the driver for the embedded explicit methods. It requires
977 * more parameters:
978 * - coarsen_param: factor multiplying the current time step when the error
979 * is below the threshold.
980 * - refine_param: factor multiplying the current time step when the error
981 * is above the threshold.
982 * - min_delta: smallest time step acceptable.
983 * - max_delta: largest time step acceptable.
984 * - refine_tol: threshold above which the time step is refined.
985 * - coarsen_tol: threshold below which the time step is coarsen.
986 *
987
988 *
989 * Embedded methods use a guessed time step. If the error using this time step
990 * is too large, the time step will be reduced. If the error is below the
991 * threshold, a larger time step will be tried for the next time step.
992 * <code>delta_t_guess</code> is the guessed time step produced by the
993 * embedded method. In summary, time step size is potentially modified in
994 * three ways:
995 * - Reducing or increasing time step size within
997 * - Using the calculated <code>delta_t_guess</code>.
998 * - Automatically adjusting the step size of the last time step to ensure
999 * simulation ends precisely at <code>final_time</code>. This adjustment
1000 * is handled inside the DiscreteTime instance.
1001 *
1002 * @code
1003 *   unsigned int Diffusion::embedded_explicit_method(
1004 *   const TimeStepping::runge_kutta_method method,
1005 *   const unsigned int n_time_steps,
1006 *   const double initial_time,
1007 *   const double final_time)
1008 *   {
1009 *   const double time_step =
1010 *   (final_time - initial_time) / static_cast<double>(n_time_steps);
1011 *   const double coarsen_param = 1.2;
1012 *   const double refine_param = 0.8;
1013 *   const double min_delta = 1e-8;
1014 *   const double max_delta = 10 * time_step;
1015 *   const double refine_tol = 1e-1;
1016 *   const double coarsen_tol = 1e-5;
1017 *  
1018 *   solution = 0.;
1019 *   constraint_matrix.distribute(solution);
1020 *  
1022 *   embedded_explicit_runge_kutta(method,
1023 *   coarsen_param,
1024 *   refine_param,
1025 *   min_delta,
1026 *   max_delta,
1027 *   refine_tol,
1028 *   coarsen_tol);
1029 *   output_results(initial_time, 0, method);
1030 *   DiscreteTime time(initial_time, final_time, time_step);
1031 *   while (time.is_at_end() == false)
1032 *   {
1033 *   const double new_time =
1034 *   embedded_explicit_runge_kutta.evolve_one_time_step(
1035 *   [this](const double time, const Vector<double> &y) {
1036 *   return this->evaluate_diffusion(time, y);
1037 *   },
1038 *   time.get_current_time(),
1039 *   time.get_next_step_size(),
1040 *   solution);
1041 *   time.set_next_step_size(new_time - time.get_current_time());
1042 *   time.advance_time();
1043 *  
1044 *   constraint_matrix.distribute(solution);
1045 *  
1046 *   if (time.get_step_number() % 10 == 0)
1047 *   output_results(time.get_current_time(),
1048 *   time.get_step_number(),
1049 *   method);
1050 *  
1051 *   time.set_desired_next_step_size(
1052 *   embedded_explicit_runge_kutta.get_status().delta_t_guess);
1053 *   }
1054 *  
1055 *   return time.get_step_number();
1056 *   }
1057 *  
1058 *  
1059 *  
1060 * @endcode
1061 *
1062 *
1063 * <a name="step_52-codeDiffusionruncode"></a>
1064 * <h4><code>Diffusion::run</code></h4>
1065 *
1066
1067 *
1068 * The following is the main function of the program. At the top, we create
1069 * the grid (a @f$[0,5]\times [0,5]@f$ square) and refine it four times to get a
1070 * mesh that has 16 by 16 cells, for a total of 256. We then set the boundary
1071 * indicator to 1 for those parts of the boundary where @f$x=0@f$ and @f$x=5@f$.
1072 *
1073 * @code
1074 *   void Diffusion::run()
1075 *   {
1078 *  
1079 *   for (const auto &cell : triangulation.active_cell_iterators())
1080 *   for (const auto &face : cell->face_iterators())
1081 *   if (face->at_boundary())
1082 *   {
1083 *   if ((face->center()[0] == 0.) || (face->center()[0] == 5.))
1084 *   face->set_boundary_id(1);
1085 *   else
1086 *   face->set_boundary_id(0);
1087 *   }
1088 *  
1089 * @endcode
1090 *
1091 * Next, we set up the linear systems and fill them with content so that
1092 * they can be used throughout the time stepping process:
1093 *
1094 * @code
1095 *   setup_system();
1096 *  
1097 *   assemble_system();
1098 *  
1099 * @endcode
1100 *
1101 * Finally, we solve the diffusion problem using several of the
1102 * Runge-Kutta methods implemented in namespace TimeStepping, each time
1103 * outputting the error at the end time. (As explained in the
1104 * introduction, since the exact solution is zero at the final time, the
1105 * error equals the numerical solution and can be computed by just taking
1106 * the @f$l_2@f$ norm of the solution vector.)
1107 *
1108 * @code
1109 *   unsigned int n_steps = 0;
1110 *   const unsigned int n_time_steps = 200;
1111 *   const double initial_time = 0.;
1112 *   const double final_time = 10.;
1113 *  
1114 *   std::cout << "Explicit methods:" << std::endl;
1115 *   explicit_method(TimeStepping::FORWARD_EULER,
1116 *   n_time_steps,
1117 *   initial_time,
1118 *   final_time);
1119 *   std::cout << " Forward Euler: error=" << solution.l2_norm()
1120 *   << std::endl;
1121 *  
1122 *   explicit_method(TimeStepping::RK_THIRD_ORDER,
1123 *   n_time_steps,
1124 *   initial_time,
1125 *   final_time);
1126 *   std::cout << " Third order Runge-Kutta: error=" << solution.l2_norm()
1127 *   << std::endl;
1128 *  
1129 *   explicit_method(TimeStepping::RK_CLASSIC_FOURTH_ORDER,
1130 *   n_time_steps,
1131 *   initial_time,
1132 *   final_time);
1133 *   std::cout << " Fourth order Runge-Kutta: error=" << solution.l2_norm()
1134 *   << std::endl;
1135 *   std::cout << std::endl;
1136 *  
1137 *  
1138 *   std::cout << "Implicit methods:" << std::endl;
1139 *   implicit_method(TimeStepping::BACKWARD_EULER,
1140 *   n_time_steps,
1141 *   initial_time,
1142 *   final_time);
1143 *   std::cout << " Backward Euler: error=" << solution.l2_norm()
1144 *   << std::endl;
1145 *  
1146 *   implicit_method(TimeStepping::IMPLICIT_MIDPOINT,
1147 *   n_time_steps,
1148 *   initial_time,
1149 *   final_time);
1150 *   std::cout << " Implicit Midpoint: error=" << solution.l2_norm()
1151 *   << std::endl;
1152 *  
1153 *   implicit_method(TimeStepping::CRANK_NICOLSON,
1154 *   n_time_steps,
1155 *   initial_time,
1156 *   final_time);
1157 *   std::cout << " Crank-Nicolson: error=" << solution.l2_norm()
1158 *   << std::endl;
1159 *  
1160 *   implicit_method(TimeStepping::SDIRK_TWO_STAGES,
1161 *   n_time_steps,
1162 *   initial_time,
1163 *   final_time);
1164 *   std::cout << " SDIRK: error=" << solution.l2_norm()
1165 *   << std::endl;
1166 *   std::cout << std::endl;
1167 *  
1168 *  
1169 *   std::cout << "Embedded explicit methods:" << std::endl;
1170 *   n_steps = embedded_explicit_method(TimeStepping::HEUN_EULER,
1171 *   n_time_steps,
1172 *   initial_time,
1173 *   final_time);
1174 *   std::cout << " Heun-Euler: error=" << solution.l2_norm()
1175 *   << std::endl;
1176 *   std::cout << " steps performed=" << n_steps << std::endl;
1177 *  
1178 *   n_steps = embedded_explicit_method(TimeStepping::BOGACKI_SHAMPINE,
1179 *   n_time_steps,
1180 *   initial_time,
1181 *   final_time);
1182 *   std::cout << " Bogacki-Shampine: error=" << solution.l2_norm()
1183 *   << std::endl;
1184 *   std::cout << " steps performed=" << n_steps << std::endl;
1185 *  
1186 *   n_steps = embedded_explicit_method(TimeStepping::DOPRI,
1187 *   n_time_steps,
1188 *   initial_time,
1189 *   final_time);
1190 *   std::cout << " Dopri: error=" << solution.l2_norm()
1191 *   << std::endl;
1192 *   std::cout << " steps performed=" << n_steps << std::endl;
1193 *  
1194 *   n_steps = embedded_explicit_method(TimeStepping::FEHLBERG,
1195 *   n_time_steps,
1196 *   initial_time,
1197 *   final_time);
1198 *   std::cout << " Fehlberg: error=" << solution.l2_norm()
1199 *   << std::endl;
1200 *   std::cout << " steps performed=" << n_steps << std::endl;
1201 *  
1202 *   n_steps = embedded_explicit_method(TimeStepping::CASH_KARP,
1203 *   n_time_steps,
1204 *   initial_time,
1205 *   final_time);
1206 *   std::cout << " Cash-Karp: error=" << solution.l2_norm()
1207 *   << std::endl;
1208 *   std::cout << " steps performed=" << n_steps << std::endl;
1209 *   }
1210 *   } // namespace Step52
1211 *  
1212 *  
1213 *  
1214 * @endcode
1215 *
1216 *
1217 * <a name="step_52-Thecodemaincodefunction"></a>
1218 * <h3>The <code>main()</code> function</h3>
1219 *
1220
1221 *
1222 * The following <code>main</code> function is similar to previous examples
1223 * and need not be commented on.
1224 *
1225 * @code
1226 *   int main()
1227 *   {
1228 *   try
1229 *   {
1230 *   Step52::Diffusion diffusion;
1231 *   diffusion.run();
1232 *   }
1233 *   catch (std::exception &exc)
1234 *   {
1235 *   std::cerr << std::endl
1236 *   << std::endl
1237 *   << "----------------------------------------------------"
1238 *   << std::endl;
1239 *   std::cerr << "Exception on processing: " << std::endl
1240 *   << exc.what() << std::endl
1241 *   << "Aborting!" << std::endl
1242 *   << "----------------------------------------------------"
1243 *   << std::endl;
1244 *   return 1;
1245 *   }
1246 *   catch (...)
1247 *   {
1248 *   std::cerr << std::endl
1249 *   << std::endl
1250 *   << "----------------------------------------------------"
1251 *   << std::endl;
1252 *   std::cerr << "Unknown exception!" << std::endl
1253 *   << "Aborting!" << std::endl
1254 *   << "----------------------------------------------------"
1255 *   << std::endl;
1256 *   return 1;
1257 *   };
1258 *  
1259 *   return 0;
1260 *   }
1261 * @endcode
1262<a name="step_52-Results"></a><h1>Results</h1>
1263
1264
1265The point of this program is less to show particular results, but instead to
1266show how it is done. This we have already demonstrated simply by discussing
1267the code above. Consequently, the output the program yields is relatively
1268sparse and consists only of the console output and the solutions given in VTU
1269format for visualization.
1270
1271The console output contains both errors and, for some of the methods, the
1272number of steps they performed:
1273@code
1274Explicit methods:
1275 Forward Euler: error=1.00883
1276 Third order Runge-Kutta: error=0.000227982
1277 Fourth order Runge-Kutta: error=1.90541e-06
1278
1279Implicit methods:
1280 Backward Euler: error=1.03428
1281 Implicit Midpoint: error=0.00862702
1282 Crank-Nicolson: error=0.00862675
1283 SDIRK: error=0.0042349
1284
1285Embedded explicit methods:
1286 Heun-Euler: error=0.0073012
1287 steps performed=284
1288 Bogacki-Shampine: error=0.000408407
1289 steps performed=181
1290 Dopri: error=0.000836695
1291 steps performed=120
1292 Fehlberg: error=0.00248922
1293 steps performed=106
1294 Cash-Karp: error=0.0787735
1295 steps performed=106
1296@endcode
1297
1298As expected the higher order methods give (much) more accurate solutions. We
1299also see that the (rather inaccurate) Heun-Euler method increased the number of
1300time steps in order to satisfy the tolerance. On the other hand, the other
1301embedded methods used a lot less time steps than what was prescribed.
1302 *
1303 *
1304<a name="step_52-PlainProg"></a>
1305<h1> The plain program</h1>
1306@include "step-52.cc"
1307*/
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
void initialize(const SparsityPattern &sparsity_pattern)
double evolve_one_time_step(const std::function< VectorType(const double, const VectorType &)> &f, const std::function< VectorType(const double, const double, const VectorType &)> &id_minus_tau_J_inverse, double t, double delta_t, VectorType &y) override
void refine_global(const unsigned int times=1)
Point< 3 > center
__global__ void set(Number *val, const Number s, const size_type N)
const Event new_time
Definition event.cc:67
void write_pvd_record(std::ostream &out, const std::vector< std::pair< double, std::string > > &times_and_names)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
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)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition l2.h:57
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
int(& functions)(const void *v1, const void *v2)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation