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