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