From 757e2eec8038cd91d01227e7fc0528ac3a1e2b95 Mon Sep 17 00:00:00 2001 From: Reza Rastak Date: Wed, 25 Mar 2020 21:54:23 -0700 Subject: [PATCH] Step-21 modified to use DiscreteTime --- doc/doxygen/tutorial/tutorial.h.in | 2 +- examples/step-21/doc/intro.dox | 4 ++ examples/step-21/doc/results.dox | 6 +++ examples/step-21/step-21.cc | 68 ++++++++++++++++------------ include/deal.II/base/discrete_time.h | 4 +- 5 files changed, 53 insertions(+), 31 deletions(-) diff --git a/doc/doxygen/tutorial/tutorial.h.in b/doc/doxygen/tutorial/tutorial.h.in index fb3fb8e87f..ca7a1dcfdd 100644 --- a/doc/doxygen/tutorial/tutorial.h.in +++ b/doc/doxygen/tutorial/tutorial.h.in @@ -271,7 +271,7 @@ * complicated block solvers. Simple explicit (forward Euler) time * stepping. *
Keywords: TensorFunction, FE_RaviartThomas, - * VectorTools::project() + * VectorTools::project(), DiscreteTime * * * diff --git a/examples/step-21/doc/intro.dox b/examples/step-21/doc/intro.dox index 2848a6bf44..5a62eecbc4 100644 --- a/examples/step-21/doc/intro.dox +++ b/examples/step-21/doc/intro.dox @@ -250,6 +250,10 @@ cell term to get an equation as follows: (S^n,\sigma)_\Omega + \triangle t \sum_K \left(F(S^n) q^{n+1}, \sigma\right)_K. @f} +We introduce an object of type DiscreteTime in order to keep track of the +current value of time and time step in the code. This class encapsulates many +complexities regarding adjusting time step size and stopping at a specified +final time. diff --git a/examples/step-21/doc/results.dox b/examples/step-21/doc/results.dox index 66345340db..1e353e26b2 100644 --- a/examples/step-21/doc/results.dox +++ b/examples/step-21/doc/results.dox @@ -1,5 +1,11 @@

Results

+The code as presented here does not actually compute the results +found on the web page. The reason is, that even on a decent +computer it runs more than a day. If you want to reproduce these +results, modify the end time of the DiscreteTime object to `250` within the +constructor of TwoPhaseFlowProblem. + If we run the program, we get the following kind of output: @code Number of active cells: 1024 diff --git a/examples/step-21/step-21.cc b/examples/step-21/step-21.cc index 6c22304499..9a45cf21fb 100644 --- a/examples/step-21/step-21.cc +++ b/examples/step-21/step-21.cc @@ -66,6 +66,10 @@ // such functionality: #include +// Additionally, we use the class DiscreteTime to perform +// operations related to time incrementation. +#include + // The last step is as in all previous programs: namespace Step21 { @@ -91,7 +95,9 @@ namespace Step21 // // The rest of the class should be pretty much obvious. The // viscosity variable stores the viscosity $\mu$ that enters - // several of the formulas in the nonlinear equations. + // several of the formulas in the nonlinear equations. The variable + // time keeps track of the time information within the + // simulation. template class TwoPhaseFlowProblem { @@ -119,8 +125,7 @@ namespace Step21 const unsigned int n_refinement_steps; - double time_step; - unsigned int timestep_number; + DiscreteTime time; double viscosity; BlockVector solution; @@ -500,7 +505,9 @@ namespace Step21 // First for the constructor. We use $RT_k \times DQ_k \times DQ_k$ // spaces. The time step is set to zero initially, but will be computed // before it is needed first, as described in a subsection of the - // introduction. + // introduction. The time object internally prevents itself from being + // incremented with $dt = 0$, forcing us to set a non-zero desired size for + // $dt$ before advancing time. template TwoPhaseFlowProblem::TwoPhaseFlowProblem(const unsigned int degree) : degree(degree) @@ -512,8 +519,7 @@ namespace Step21 1) , dof_handler(triangulation) , n_refinement_steps(5) - , time_step(0) - , timestep_number(1) + , time(/*start time*/ 0., /*end time*/ 1., /*time step*/ 0.) , viscosity(0.2) {} @@ -827,10 +833,11 @@ namespace Step21 const Tensor<1, dim> grad_phi_i_s = fe_values[saturation].gradient(i, q); - local_rhs(i) += (time_step * fractional_flow(old_s, viscosity) * - present_u * grad_phi_i_s + - old_s * phi_i_s) * - fe_values.JxW(q); + local_rhs(i) += + (time.get_next_step_size() * fractional_flow(old_s, viscosity) * + present_u * grad_phi_i_s + + old_s * phi_i_s) * + fe_values.JxW(q); } // Secondly, we have to deal with the flux parts on the face @@ -884,7 +891,7 @@ namespace Step21 for (unsigned int i = 0; i < dofs_per_cell; ++i) local_rhs(i) -= - time_step * normal_flux * + time.get_next_step_size() * normal_flux * fractional_flow((is_outflow_q_point == true ? old_solution_values_face[q](dim + 1) : neighbor_saturation[q]), @@ -963,9 +970,11 @@ namespace Step21 // // The maximal velocity we compute using a helper function to compute the // maximal velocity defined below, and with all this we can evaluate our - // new time step length: - time_step = - std::pow(0.5, double(n_refinement_steps)) / get_maximal_velocity(); + // new time step length. The method DiscreteTime::set_next_time_step() is + // used to assign the new value of the time step to the DiscreteTime + // object. + time.set_next_step_size(std::pow(0.5, double(n_refinement_steps)) / + get_maximal_velocity()); // The next step is to assemble the right hand side, and then to pass // everything on for solution. At the end, we project back saturations @@ -1005,7 +1014,7 @@ namespace Step21 template void TwoPhaseFlowProblem::output_results() const { - if (timestep_number % 5 != 0) + if (time.get_step_number() % 5 != 0) return; std::vector solution_names; @@ -1031,7 +1040,8 @@ namespace Step21 data_out.build_patches(degree + 1); std::ofstream output("solution-" + - Utilities::int_to_string(timestep_number, 4) + ".vtk"); + Utilities::int_to_string(time.get_step_number(), 4) + + ".vtk"); data_out.write_vtk(output); } @@ -1118,13 +1128,15 @@ namespace Step21 // // The second point worth mentioning is that we only compute the length of // the present time step in the middle of solving the linear system - // corresponding to each time step. We can therefore output the present end + // corresponding to each time step. We can therefore output the present // time of a time step only at the end of the time step. - // - // The function as it is here does actually not compute the results - // found on the web page. The reason is, that even on a decent - // computer it runs more than a day. If you want to reproduce these - // results, set the final time at the end of the do loop to 250. + // We increment time by calling the method DiscreteTime::advance_time() + // inside the loop. Since we are reporting the time and dt after we + // increment it, we have to call the method + // DiscreteTime::get_previous_step_size() instead of + // DiscreteTime::get_next_step_size(). After many steps, when the simulation + // reaches the end time, the last dt is chosen by the DiscreteTime class in + // such a way that the last step finishes exactly at the end time. template void TwoPhaseFlowProblem::run() { @@ -1141,11 +1153,9 @@ namespace Step21 old_solution); } - double time = 0; - do { - std::cout << "Timestep " << timestep_number << std::endl; + std::cout << "Timestep " << time.get_step_number() + 1 << std::endl; assemble_system(); @@ -1153,13 +1163,13 @@ namespace Step21 output_results(); - time += time_step; - ++timestep_number; - std::cout << " Now at t=" << time << ", dt=" << time_step << '.' + time.advance_time(); + std::cout << " Now at t=" << time.get_current_time() + << ", dt=" << time.get_previous_step_size() << '.' << std::endl << std::endl; } - while (time <= 1.); + while (time.is_at_end() == false); } } // namespace Step21 diff --git a/include/deal.II/base/discrete_time.h b/include/deal.II/base/discrete_time.h index de7fa92bc2..c5ad7dc3d7 100644 --- a/include/deal.II/base/discrete_time.h +++ b/include/deal.II/base/discrete_time.h @@ -24,7 +24,9 @@ DEAL_II_NAMESPACE_OPEN * This class provides a means to keep track of the simulation time in a * time-dependent simulation. It manages stepping forward from a start time * $T_{\text{start}}$ to an end time $T_{\text{end}}$. It also allows adjusting - * the time step size during the simulation. + * the time step size during the simulation. This class provides the necessary + * interface to be incorporated in any time-dependent simulation. As an + * example, the usage of this class is demonstrated in step-21. * * This class provides a number of invariants that are guaranteed to be * true at all times. -- 2.39.5