]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Step-21 modified to use DiscreteTime
authorReza Rastak <rastak@stanford.edu>
Thu, 26 Mar 2020 04:54:23 +0000 (21:54 -0700)
committerReza Rastak <rastak@stanford.edu>
Thu, 23 Apr 2020 01:35:01 +0000 (18:35 -0700)
doc/doxygen/tutorial/tutorial.h.in
examples/step-21/doc/intro.dox
examples/step-21/doc/results.dox
examples/step-21/step-21.cc
include/deal.II/base/discrete_time.h

index fb3fb8e87ff7ea84f1e57e7f4d13e4d3abd937a6..ca7a1dcfddf7ea90622168b7bdf37cca37ae2b77 100644 (file)
  *       complicated block solvers. Simple explicit (forward Euler) time 
  *       stepping.
  *       <br/> Keywords: TensorFunction, FE_RaviartThomas,
- *       VectorTools::project()
+ *       VectorTools::project(), DiscreteTime
  *       </td></tr>
  *
  *   <tr valign="top">
index 2848a6bf44e811db2b613aaeb9f36bfd38771d7d..5a62eecbc4d9fac8c9643fe63ed4364f6ef3af5b 100644 (file)
@@ -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.
 
 
 
index 66345340db832c7e77514956275f5bf40f4f4251..1e353e26b2b0d4837e64ed5c47db41001fb8eb20 100644 (file)
@@ -1,5 +1,11 @@
 <h1>Results</h1>
 
+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
index 6c22304499dd46f4beb5dc79234030a1a37179c9..9a45cf21fbf07a843ddfdcd89822646f0565850f 100644 (file)
 // such functionality:
 #include <deal.II/base/tensor_function.h>
 
+// Additionally, we use the class <code>DiscreteTime</code> to perform
+// operations related to time incrementation.
+#include <deal.II/base/discrete_time.h>
+
 // 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
   // <code>viscosity</code> 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
+  // <code>time</code> keeps track of the time information within the
+  // simulation.
   template <int dim>
   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<double> 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 <int dim>
   TwoPhaseFlowProblem<dim>::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 <int dim>
   void TwoPhaseFlowProblem<dim>::output_results() const
   {
-    if (timestep_number % 5 != 0)
+    if (time.get_step_number() % 5 != 0)
       return;
 
     std::vector<std::string> 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 <int dim>
   void TwoPhaseFlowProblem<dim>::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
 
index de7fa92bc2ef708e3aa9a8bac07913a85dc8f38c..c5ad7dc3d7c1b4519ea245a0bd6766d3d5bfdfde 100644 (file)
@@ -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.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.