716 * mass_minus_tau_Jacobian.add(-tau, system_matrix);
718 * inverse_mass_minus_tau_Jacobian.
initialize(mass_minus_tau_Jacobian);
724 * inverse_mass_minus_tau_Jacobian.
vmult(result, tmp);
734 * <a name=
"codeDiffusionoutput_resultscode"></a>
735 * <h4><code>Diffusion::output_results</code></h4>
739 * The following
function then outputs the solution in
vtu files indexed by
740 * the number of the time step and the name of the time stepping method. Of
741 * course, the (exact) result should really be the same
for all time
742 * stepping method, but the output here at least allows us to compare them.
745 *
void Diffusion::output_results(
const double time,
746 *
const unsigned int time_step,
749 * std::string method_name;
755 * method_name =
"forward_euler";
760 * method_name =
"rk3";
765 * method_name =
"rk4";
770 * method_name =
"backward_euler";
775 * method_name =
"implicit_midpoint";
780 * method_name =
"sdirk";
785 * method_name =
"heun_euler";
790 * method_name =
"bocacki_shampine";
795 * method_name =
"dopri";
800 * method_name =
"fehlberg";
805 * method_name =
"cash_karp";
823 *
const std::string filename =
"solution_" + method_name +
"-" +
826 * std::ofstream output(filename);
829 *
static std::vector<std::pair<double, std::string>> times_and_names;
831 *
static std::string method_name_prev =
"";
832 *
static std::string pvd_filename;
833 *
if (method_name_prev != method_name)
835 * times_and_names.clear();
836 * method_name_prev = method_name;
837 * pvd_filename =
"solution_" + method_name +
".pvd";
839 * times_and_names.emplace_back(time, filename);
840 * std::ofstream pvd_output(pvd_filename);
848 * <a name=
"codeDiffusionexplicit_methodcode"></a>
849 * <h4><code>Diffusion::explicit_method</code></h4>
853 * This
function is the driver
for all the
explicit methods. At the
854 * top it initializes the time stepping and the solution (by setting
855 * it to
zero and then ensuring that boundary
value and hanging node
856 * constraints are respected; of course, with the mesh we use here,
857 * hanging node constraints are not in fact an issue). It then calls
858 * <code>evolve_one_time_step</code> which performs
one time step.
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
872 * Finally, the solution is output
873 * every 10 time steps.
877 *
const unsigned int n_time_steps,
878 *
const double initial_time,
879 *
const double final_time)
881 *
const double time_step =
882 * (final_time - initial_time) /
static_cast<double>(n_time_steps);
883 *
double time = initial_time;
886 * constraint_matrix.distribute(solution);
890 * output_results(time, 0, method);
891 *
for (
unsigned int i = 0; i < n_time_steps; ++i)
893 * time = explicit_runge_kutta.evolve_one_time_step(
895 *
return this->evaluate_diffusion(time, y);
901 * constraint_matrix.distribute(solution);
903 *
if ((i + 1) % 10 == 0)
904 * output_results(time, i + 1, method);
913 * <a name=
"codeDiffusionimplicit_methodcode"></a>
914 * <h4><code>Diffusion::implicit_method</code></h4>
915 * This
function is equivalent to <code>explicit_method</code> but
for
916 * implicit methods. When
using implicit methods, we need to evaluate
917 * @f$M^{-1}(f(t,y))@f$ and @f$\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial
918 * y}\right)^{-1}@f$
for which we use the two member
functions previously
923 *
const unsigned int n_time_steps,
924 *
const double initial_time,
925 *
const double final_time)
927 *
const double time_step =
928 * (final_time - initial_time) /
static_cast<double>(n_time_steps);
929 *
double time = initial_time;
932 * constraint_matrix.distribute(solution);
936 * output_results(time, 0, method);
937 *
for (
unsigned int i = 0; i < n_time_steps; ++i)
939 * time = implicit_runge_kutta.evolve_one_time_step(
941 *
return this->evaluate_diffusion(time, y);
943 * [
this](
const double time,
const double tau,
const Vector<double> &y) {
944 *
return this->id_minus_tau_J_inverse(time, tau, y);
950 * constraint_matrix.distribute(solution);
952 *
if ((i + 1) % 10 == 0)
953 * output_results(time, i + 1, method);
962 * <a name=
"codeDiffusionembedded_explicit_methodcode"></a>
963 * <h4><code>Diffusion::embedded_explicit_method</code></h4>
964 * This
function is the driver
for the embedded
explicit methods. It requires
966 * - coarsen_param: factor multiplying the current time step when the error
967 * is below the threshold.
968 * - refine_param: factor multiplying the current time step when the error
969 * is above the threshold.
970 * - min_delta: smallest time step acceptable.
971 * - max_delta: largest time step acceptable.
972 * - refine_tol: threshold above which the time step is refined.
973 * - coarsen_tol: threshold below which the time step is
coarsen.
974 * Embedded methods use a guessed time step. If the error
using this time step
975 * is too large, the time step will be reduced. If the error is below the
976 * threshold, a larger time step will be tried
for the next time step.
977 * <code>delta_t_guess</code> is the guessed time step produced by the
981 *
unsigned int Diffusion::embedded_explicit_method(
983 *
const unsigned int n_time_steps,
984 *
const double initial_time,
985 *
const double final_time)
988 * (final_time - initial_time) /
static_cast<double>(n_time_steps);
989 *
double time = initial_time;
990 *
const double coarsen_param = 1.2;
991 *
const double refine_param = 0.8;
992 *
const double min_delta = 1
e-8;
993 *
const double max_delta = 10 * time_step;
994 *
const double refine_tol = 1
e-1;
995 *
const double coarsen_tol = 1
e-5;
998 * constraint_matrix.distribute(solution);
1001 * embedded_explicit_runge_kutta(method,
1008 * output_results(time, 0, method);
1012 * Now
for the time
loop. The last time step is chosen such that the
final
1013 * time is exactly reached.
1016 *
unsigned int n_steps = 0;
1017 *
while (time < final_time)
1019 *
if (time + time_step > final_time)
1020 * time_step = final_time - time;
1022 * time = embedded_explicit_runge_kutta.evolve_one_time_step(
1024 *
return this->evaluate_diffusion(time, y);
1030 * constraint_matrix.distribute(solution);
1032 *
if ((n_steps + 1) % 10 == 0)
1033 * output_results(time, n_steps + 1, method);
1035 * time_step = embedded_explicit_runge_kutta.get_status().delta_t_guess;
1047 * <a name=
"codeDiffusionruncode"></a>
1052 * The following is the main
function of the program. At the top, we create
1053 * the grid (a [0,5]x[0,5] square) and
refine it four times to get a mesh
1054 * that has 16 by 16 cells,
for a total of 256. We then
set the boundary
1055 * indicator to 1
for those parts of the boundary where @f$x=0@f$ and @f$x=5@f$.
1063 *
for (
const auto &cell :
triangulation.active_cell_iterators())
1064 *
for (
const auto &face : cell->face_iterators())
1065 *
if (face->at_boundary())
1067 *
if ((face->center()[0] == 0.) || (face->center()[0] == 5.))
1068 * face->set_boundary_id(1);
1070 * face->set_boundary_id(0);
1075 * Next, we
set up the linear systems and fill them with content so that
1076 * they can be used throughout the time stepping process:
1081 * assemble_system();
1085 * Finally, we solve the diffusion problem
using several of the
1086 * Runge-Kutta methods implemented in
namespace TimeStepping, each time
1087 * outputting the error at the
end time. (As explained in the
1088 * introduction, since the exact solution is
zero at the
final time, the
1089 * error equals the numerical solution and can be computed by just taking
1090 * the @f$l_2@f$
norm of the solution vector.)
1093 *
unsigned int n_steps = 0;
1094 *
const unsigned int n_time_steps = 200;
1095 *
const double initial_time = 0.;
1096 *
const double final_time = 10.;
1098 * std::cout <<
"Explicit methods:" << std::endl;
1103 * std::cout <<
" Forward Euler: error=" << solution.l2_norm()
1110 * std::cout <<
" Third order Runge-Kutta: error=" << solution.l2_norm()
1117 * std::cout <<
" Fourth order Runge-Kutta: error=" << solution.l2_norm()
1119 * std::cout << std::endl;
1122 * std::cout <<
"Implicit methods:" << std::endl;
1127 * std::cout <<
" Backward Euler: error=" << solution.l2_norm()
1134 * std::cout <<
" Implicit Midpoint: error=" << solution.l2_norm()
1141 * std::cout <<
" Crank-Nicolson: error=" << solution.l2_norm()
1148 * std::cout <<
" SDIRK: error=" << solution.l2_norm()
1150 * std::cout << std::endl;
1153 * std::cout <<
"Embedded explicit methods:" << std::endl;
1158 * std::cout <<
" Heun-Euler: error=" << solution.l2_norm()
1160 * std::cout <<
" steps performed=" << n_steps << std::endl;
1166 * std::cout <<
" Bogacki-Shampine: error=" << solution.l2_norm()
1168 * std::cout <<
" steps performed=" << n_steps << std::endl;
1174 * std::cout <<
" Dopri: error=" << solution.l2_norm()
1176 * std::cout <<
" steps performed=" << n_steps << std::endl;
1182 * std::cout <<
" Fehlberg: error=" << solution.l2_norm()
1184 * std::cout <<
" steps performed=" << n_steps << std::endl;
1190 * std::cout <<
" Cash-Karp: error=" << solution.l2_norm()
1192 * std::cout <<
" steps performed=" << n_steps << std::endl;
1201 * <a name=
"Thecodemaincodefunction"></a>
1202 * <h3>The <code>main()</code>
function</h3>
1206 * The following <code>main</code>
function is similar to previous examples
1207 * and need not be commented on.
1214 * Step52::Diffusion diffusion;
1217 *
catch (std::exception &exc)
1219 * std::cerr << std::endl
1221 * <<
"----------------------------------------------------"
1223 * std::cerr <<
"Exception on processing: " << std::endl
1224 * << exc.what() << std::endl
1225 * <<
"Aborting!" << std::endl
1226 * <<
"----------------------------------------------------"
1232 * std::cerr << std::endl
1234 * <<
"----------------------------------------------------"
1236 * std::cerr <<
"Unknown exception!" << std::endl
1237 * <<
"Aborting!" << std::endl
1238 * <<
"----------------------------------------------------"
1246 <a name=
"Results"></a><h1>Results</h1>
1249 The
point of
this program is less to show particular results, but instead to
1250 show how it is done. This we have already demonstrated simply by discussing
1251 the code above. Consequently, the output the program yields is relatively
1252 sparse and consists only of the console output and the solutions given in VTU
1253 format
for visualization.
1255 The console output contains both errors and,
for some of the methods, the
1256 number of steps they performed:
1259 Forward Euler: error=1.00883
1260 Third order Runge-Kutta: error=0.000227982
1261 Fourth order Runge-Kutta: error=1.90541e-06
1264 Backward Euler: error=1.03428
1265 Implicit Midpoint: error=0.00862702
1266 Crank-Nicolson: error=0.00862675
1267 SDIRK: error=0.0042349
1269 Embedded %
explicit methods:
1270 Heun-Euler: error=0.0073012
1272 Bogacki-Shampine: error=0.000403281
1274 Dopri: error=0.0165485
1276 Fehlberg: error=0.00104926
1278 Cash-Karp: error=8.59366e-07
1283 As expected the higher order methods give (much) more accurate solutions. We
1284 also see that the (rather inaccurate) Heun-Euler method increased the number of
1285 time steps in order to satisfy the tolerance. On the other hand, the other
1286 embedded methods used a lot less time steps than what was prescribed.
1289 <a name=
"PlainProg"></a>
1290 <h1> The plain program</h1>
1291 @include
"step-52.cc"