]> https://gitweb.dealii.org/ - dealii.git/commitdiff
ThetaTimestepping implemented and ready for testing
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 13 Jan 2010 15:21:48 +0000 (15:21 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 13 Jan 2010 15:21:48 +0000 (15:21 +0000)
git-svn-id: https://svn.dealii.org/trunk@20359 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/event.h
deal.II/deal.II/include/numerics/newton.h
deal.II/deal.II/include/numerics/operator.h
deal.II/deal.II/include/numerics/operator.templates.h
deal.II/deal.II/include/numerics/theta_timestepping.h [new file with mode: 0644]
deal.II/deal.II/include/numerics/theta_timestepping.templates.h [new file with mode: 0644]
deal.II/deal.II/include/numerics/time_step_control.h [new file with mode: 0644]
deal.II/deal.II/source/numerics/operators.all_dimensions.cc [moved from deal.II/deal.II/source/numerics/operators.cc with 79% similarity]
deal.II/deal.II/source/numerics/operators.inst.in
deal.II/deal.II/source/numerics/timestep_control.all_dimensions.cc [new file with mode: 0644]

index 22a83662161cc71f6c14f7574b233cdc7b09d00c..84f5bd51daff0ae2dcb43a8dcf7fcfde08f9a72e 100644 (file)
@@ -137,12 +137,14 @@ namespace Algorithms
 /// Events used by library operators  
   namespace Events
   {
+/// The program has just started and everything should be new
+    extern const Event initial;
 /// The current derivative leads to slow convergence of Newton's method
     extern const Event bad_derivative;
 /// The time stepping scheme starts a new time step
     extern const Event new_time;
 /// The time stepping scheme has changed the time step size
-    extern const Event new_time_step_size;
+    extern const Event new_timestep_size;
   }
   
 
index a6882c739e0b605e530379487c9e9375937a7bda..97ba03865219f807c4fc4eee114ee5fbfd6be15c 100644 (file)
@@ -15,6 +15,7 @@
 
 #include <base/smartpointer.h>
 #include <lac/solver_control.h>
+#include <numerics/operator.h>
 
 DEAL_II_NAMESPACE_OPEN
 
index 2f1c799af1f86830f061d65610d62adf9e154bce..a43b6acdbc2ad3805eabbc41d6925eff8f7b3e02 100644 (file)
@@ -70,7 +70,41 @@ namespace Algorithms
                                        */
       Event notifications;      
   };
+  
+/**
+ * An unary operator base class, intended to output the vectors in
+ * NamedData in each step of an iteration.
+ *
+ * @author Guido Kanschat, 2010
+ */
+  template <class VECTOR>
+  class OutputOperator : public Subscriptor
+  {
+    public:
+                                      /**
+                                       * Empty virtual destructor.
+                                       */
+      virtual ~OutputOperator();
+                                      /**
+                                       * Set the current step.
+                                       */
+      OutputOperator<VECTOR>& operator<< (unsigned int step);
 
+                                      /**
+                                       * Output all the vectors in NamedData.
+                                       */
+      virtual OutputOperator<VECTOR>& operator<< (const NamedData<VECTOR*> vectors);
+    protected:
+      unsigned int step;
+  };
+
+  template <class VECTOR>
+  OutputOperator<VECTOR>&
+  OutputOperator<VECTOR>::operator<< (unsigned int s)
+  {
+    step = s;
+    return *this;
+  }
 }
 
 
index 48c9ddac46933e1e2fd6dce36aebb07f2360e92a..20b7bc96dae5460b5f18ec5bfdc098ff021917af 100644 (file)
@@ -34,6 +34,26 @@ namespace Algorithms
   {
     notifications.clear();
   }
+
+
+  template <class VECTOR>
+  OutputOperator<VECTOR>::~OutputOperator()
+  {}
+  
+  
+  template <class VECTOR>
+  OutputOperator<VECTOR>&
+  OutputOperator<VECTOR>::operator<< (const NamedData<VECTOR*> vectors)
+  {
+    std::cout << ' ' << step;
+    for (unsigned int i=0;i<vectors.size();++i)
+      vectors(i)->print(std::cout);
+    std::cout << std::endl;
+    return *this;
+  }
+  
+  
+  
 }
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/deal.II/include/numerics/theta_timestepping.h b/deal.II/deal.II/include/numerics/theta_timestepping.h
new file mode 100644 (file)
index 0000000..8386033
--- /dev/null
@@ -0,0 +1,254 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#ifndef __deal2__theta_timestepping_h
+#define __deal2__theta_timestepping_h
+
+#include <base/smartpointer.h>
+#include <numerics/operator.h>
+#include <numerics/time_step_control.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+class ParameterHandler;
+
+namespace Algorithms
+{
+/**
+ * A little structure, gathering the size of a timestep and the
+ * current time. Time stepping schemes can use this to provide time
+ * step information to the classes actually performing a single step.
+ *
+ * The definition of what is considered "current time" depends on the
+ * scheme. For an explicit scheme, this is the time at the beginning
+ * of the step. For an implicit scheme, it is usually the time at the
+ * end.
+ */
+  struct TimeStepData
+  {
+/// The current time
+      double time;
+/// The current step size
+      double step;
+  };
+  
+/**
+ * Application class performing the theta timestepping scheme.
+ *
+ * The theta scheme is an abstraction of implicit and explicit Euler
+ * schemes, the Crank-Nicholson scheme and linear combinations of
+ * those. The choice of the actual scheme is controlled by the
+ * parameter #theta as follows.
+ * <ul>
+ * <li> #theta=0: explicit Euler scheme
+ * <li> #theta=1: implicit Euler scheme
+ * <li> #theta=½: Crank-Nicholson scheme
+ * </ul>
+ *
+ * For fixed #theta, the Crank-Nicholson scheme is the only second
+ * order scheme. Nevertheless, further stability may be achieved by
+ * choosing #theta larger than ½, thereby introducing a first order
+ * error term. In order to avoid a loss of convergence order, the
+ * adaptive theta scheme can be used, where <i>#theta=½+c dt</i>.
+ *
+ * Assume that we want to solve the equation <i>u' = Au</i> with a
+ * step size <i>k</i>.  A step of the theta scheme can be written as
+ *
+ * @f[
+ *   (M - \theta k A) u_{n+1} = (M + (1-\theta)k A) u_n.
+ * @f]
+ *
+ * Here, <i>M</i> is the mass matrix. We see, that the right hand side
+ * amounts to an explicit Euler step with modified step size in weak
+ * form (up to inversion of M). The left hand side corresponds to an
+ * implicit Euler step with modified step size (right hand side
+ * given). Thus, the implementation of the theta scheme will use two
+ * Operator objects, one for the explicit, one for the implicit
+ * part. Each of these will use its own TimeStepData to account for
+ * the modified step sizes (and different times if the problem is not
+ * autonomous).
+ *
+ * @author Guido Kanschat, 2010
+ */
+  template <class VECTOR>
+  class ThetaTimestepping : public Operator<VECTOR>
+  {
+    public:
+                                      /**
+                                       * Constructor, receiving the
+                                       * two operators stored in
+                                       * #op_explicit and #op_implicit. For
+                                       * their meening, see the
+                                       * description of those variables.
+                                       */
+      ThetaTimestepping (Operator<VECTOR>& op_explicit,
+                        Operator<VECTOR>& op_implicit);
+
+                                      /**
+                                       * The timestepping scheme. <tt>in</tt>
+                                       * should contain the initial
+                                       * value in first position. <tt>out</tt>
+                                       */
+      virtual void operator() (NamedData<VECTOR*>& out, const NamedData<VECTOR*>& in);
+
+                                      /**
+                                       * Register an event triggered
+                                       * by an outer iteration.
+                                       */
+      virtual void notify(const Event&);
+
+                                      /**
+                                       * Define an operator which
+                                       * will output the result in
+                                       * each step. Note that no
+                                       * output will be generated
+                                       * without this.
+                                       */
+      void set_output(OutputOperator<VECTOR>& output);
+      
+      void declare_parameters (ParameterHandler& param);
+      void initialize (ParameterHandler& param);
+                                      /**
+                                       * The current time in the
+                                       * timestepping scheme.
+                                       */
+      const double& current_time() const;
+                                      /**
+                                       * The current step size.
+                                       */
+      const double& step_size() const;
+                                      /**
+                                       * The weight between implicit
+                                       * and explicit part.
+                                       */
+      const double& theta() const;
+
+                                      /**
+                                       * The data handed to the
+                                       * #explicit time stepping
+                                       * operator.
+                                       *
+                                       * The time in here is the time
+                                       * at the beginning of the
+                                       * current step, the time step
+                                       * is (1-#theta) times the
+                                       * actual time step.
+                                       */
+      const TimeStepData& explicit_data() const;
+      
+                                      /**
+                                       * The data handed to the
+                                       * #implicit time stepping
+                                       * operator.
+                                       *
+                                       * The time in here is the time
+                                       * at the beginning of the
+                                       * current step, the time step
+                                       * is #theta times the
+                                       * actual time step.
+                                       */
+      const TimeStepData& implicit_data() const;
+
+                                      /**
+                                       * Allow access to the control object.
+                                       */
+      TimestepControl& timestep_control();
+      
+    private:
+                                      /**
+                                       * The object controlling the
+                                       * time step size and computing
+                                       * the new time in each step.
+                                       */
+      TimestepControl control;
+      
+                                      /**
+                                       * The control parameter theta in the
+                                       * range <tt>[0,1]</tt>.
+                                       */
+      double vtheta;
+                                      /**
+                                       * Use adaptive #theta if
+                                       * <tt>true</tt>.
+                                       */
+      bool adaptive;
+
+                                      /**
+                                       * The data for the explicit
+                                       * part of the scheme.
+                                       */
+      TimeStepData d_explicit;
+      
+                                      /**
+                                       * The data for the implicit
+                                       * part of the scheme.
+                                       */
+      TimeStepData d_implicit;
+      
+      
+                                      /**
+                                       * The operator computing the
+                                       * explicit part of the
+                                       * scheme. This will receive in
+                                       * its input data the value at
+                                       * the current time with name
+                                       * "Current time solution". It
+                                       * should obtain the current
+                                       * time and time step size from
+                                       * explicit_data().
+                                       *
+                                       * Its return value is
+                                       * <i>Mu+cAu</i>, where
+                                       * <i>u</i> is the current
+                                       * state vector, <i>M</i> the
+                                       * mass matrix, <i>A</i> the
+                                       * operator in space and
+                                       * <i>c</i> is the time step
+                                       * size in explicit_data().
+                                       */
+      SmartPointer<Operator<VECTOR> > op_explicit;
+      
+                                      /**
+                                       * The operator solving the
+                                       * implicit part of the
+                                       * scheme. It will receive in
+                                       * its input data the vector
+                                       * "Previous time
+                                       * data". Information on the
+                                       * timestep should be obtained
+                                       * from implicit_data().
+                                       *
+                                       * Its return value is the
+                                       * solution <i>u</i> of
+                                       * <i>Mu-cAu=f</i>, where
+                                       * <i>f</i> is the dual space
+                                       * vector found in the
+                                       * "Previous time" entry of the
+                                       * input data, <i>M</i> the
+                                       * mass matrix, <i>A</i> the
+                                       * operator in space and
+                                       * <i>c</i> is the time step
+                                       * size in explicit_data().
+                                       */
+      SmartPointer<Operator<VECTOR> > op_implicit;
+
+                                      /**
+                                       * The operator writing the
+                                       * output in each time step
+                                       */
+      SmartPointer<OutputOperator<VECTOR> > output;
+  };
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/deal.II/include/numerics/theta_timestepping.templates.h b/deal.II/deal.II/include/numerics/theta_timestepping.templates.h
new file mode 100644 (file)
index 0000000..2d324b8
--- /dev/null
@@ -0,0 +1,123 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2006, 2007, 2008, 2009, 2010 by Guido Kanschat
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <numerics/theta_timestepping.h>
+#include <lac/vector_memory.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Algorithms
+{
+  template <class VECTOR>
+  ThetaTimestepping<VECTOR>::ThetaTimestepping (Operator<VECTOR>& e, Operator<VECTOR>& i)
+                 : op_explicit(&e), op_implicit(&i)
+  {}
+
+  
+  template <class VECTOR>
+  void
+  ThetaTimestepping<VECTOR>::notify(const Event& e)
+  {
+    op_explicit->notify(e);
+    op_implicit->notify(e);
+  }
+
+  template <class VECTOR>
+  void
+  ThetaTimestepping<VECTOR>::declare_parameters(ParameterHandler& param)
+  {
+    param.enter_subsection("ThetaTimestepping");
+    TimestepControl::declare_parameters (param);
+    param.declare_entry("Theta", ".5", Patterns::Double());
+    param.declare_entry("Adaptive", "false", Patterns::Bool());
+    param.leave_subsection();
+  }
+
+  template <class VECTOR>
+  void
+  ThetaTimestepping<VECTOR>::initialize (ParameterHandler& param)
+  {
+    param.enter_subsection("ThetaTimestepping");
+    control.parse_parameters (param);
+    vtheta = param.get_double("Theta");
+    adaptive = param.get_bool("Adaptive");
+    param.leave_subsection ();
+  }
+
+
+  template <class VECTOR>
+  void
+  ThetaTimestepping<VECTOR>::operator() (NamedData<VECTOR*>& out, const NamedData<VECTOR*>& in)
+  {
+    Assert(!adaptive, ExcNotImplemented());
+  
+    deallog.push ("Theta");
+    GrowingVectorMemory<VECTOR> mem;
+    typename VectorMemory<VECTOR>::Pointer aux(mem);
+    aux->reinit(*out(0));
+    
+    control.restart();
+  
+    bool first = true;
+    d_explicit.time = control.now();
+
+                                    // The data used to compute the
+                                    // vector associated with the old
+                                    // timestep
+    VECTOR* p = out(0);
+    NamedData<VECTOR*> src1;
+    src1.add(p, "Previous time");
+    src1.merge(in);
+    p = aux;
+    NamedData<VECTOR*> out1;
+    out1.add(p, "Result");
+                                    // The data provided to the inner
+                                    // solver
+    NamedData<VECTOR*> src2;
+    src2.add(p, "Previous time data");
+    src2.merge(in);
+    
+    if (output != 0)
+      (*output) << 0U << out;
+    
+    for (unsigned int count = 1; d_explicit.time < control.final(); ++count)
+      {      
+       const bool step_change = control.advance();
+       d_implicit.time = control.now();
+       d_explicit.step = (1.-vtheta)*control.step();
+       d_implicit.step = vtheta*control.step();
+       deallog << "Time step:" << d_implicit.time << std::endl;
+
+       op_explicit->notify(Events::new_time);
+       op_implicit->notify(Events::new_time);
+       if (step_change)
+         {
+           op_explicit->notify(Events::new_timestep_size);
+           op_implicit->notify(Events::new_timestep_size);
+         }
+      
+                                        // Compute
+                                        // (I + (1-theta)dt A) u
+       (*op_explicit)(out1, src1);
+       (*op_implicit)(out, src2);
+      
+       if (output != 0 && control.print())
+         (*output) << count << out;
+       
+       first = false;
+       d_explicit.time = control.now();
+      }
+    deallog.pop();
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/deal.II/include/numerics/time_step_control.h b/deal.II/deal.II/include/numerics/time_step_control.h
new file mode 100644 (file)
index 0000000..4646703
--- /dev/null
@@ -0,0 +1,323 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2010 by Guido Kanschat
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#ifndef __deal2__time_step_control_h
+#define __deal2__time_step_control_h
+
+#include <base/subscriptor.h>
+#include <base/smartpointer.h>
+#include <lac/vector_memory.h>
+#include <cstdio>
+
+DEAL_II_NAMESPACE_OPEN
+
+class ParameterHandler;
+
+namespace Algorithms
+{
+/**
+ * Control class for timestepping schemes. Its main task is
+ * determining the size of the next time step and the according point
+ * in the time interval. Additionally, it controls writing the
+ * solution to a file.
+ *
+ * The size of the next time step is determined as follows:
+ * <ol>
+ * <li> According to the strategy, the step size is tentatively added to the
+ * current time.
+ * <li> If the resulting time exceeds the final time of the interval,
+ * the step size is reduced in order to meet this time.
+ * <li> If the resulting time is below the final time by just a
+ * fraction of the step size, the step size is increased in order to
+ * meet this time.
+ * <li> The resulting step size is used from the current time.
+ * </ol>
+ *
+ * The variable #print_step can be used to control the amount of
+ * output generated by the timestepping scheme.
+ */
+  class TimestepControl : public Subscriptor
+  {
+    public:
+                                      /**
+                                       * The time stepping
+                                       * strategies. These are
+                                       * controlled by the value of
+                                       * tolerance() and start_step().
+                                       */
+      enum Strategy
+      {
+                                            /**
+                                             * Choose a uniform time
+                                             * step size. The step size
+                                             * is determined by
+                                             * start_step(),
+                                             * tolerance() is ignored.
+                                             */
+           uniform,
+                                            /**
+                                             * Start with the time step
+                                             * size given by
+                                             * start_step() and double it
+                                             * in every
+                                             * step. tolerance() is
+                                             * ignored.
+                                             *
+                                             * This strategy is
+                                             * intended for
+                                             * pseudo-timestepping
+                                             * schemes computing a
+                                             * stationary limit.
+                                             */
+           doubling
+      };
+
+                                      /**
+                                       * Constructor setting default
+                                       * values
+                                       */
+      TimestepControl (double start = 0.,
+                      double final = 1.,
+                      double tolerance = 1.e-2,
+                      double start_step = 1.e-2,
+                      double print_step = -1.,
+                      double max_step = 1.);
+
+                                      /**
+                                       * Declare the control
+                                       * parameters for parameter
+                                       * handler.
+                                       */
+      static void declare_parameters (ParameterHandler& param);
+                                      /**
+                                       * Read the control parameters
+                                       * from a parameter handler.
+                                       */
+      void parse_parameters (ParameterHandler& param);
+
+                                      /**
+                                       * The left end of the time interval.
+                                       */
+      double start () const;
+                                      /**
+                                       * The right end of the time
+                                       * interval. The control
+                                       * mechanism ensures that the
+                                       * final time step ends at this
+                                       * point.
+                                       */
+      double final () const;
+                                      /**
+                                       * The tolerance value
+                                       * controlling the time steps.
+                                       */
+      double tolerance () const;
+                                      /**
+                                       * The size of the current time
+                                       * step.
+                                       */
+      double step () const;
+      
+                                      /**
+                                       * The current time.
+                                       */
+      double now () const;
+
+                                      /**
+                                       * Compute the size of the next
+                                       * step size and return true if
+                                       * it differs from the current
+                                       * step size. Advance the
+                                       * current time by the new step
+                                       * size.
+                                       */
+      bool advance ();
+
+                                      /**
+                                       * Set start value.
+                                       */
+      void start (double);
+                                      /**
+                                       * Set final time value.
+                                       */
+      void final (double);
+                                      /**
+                                       * Set tolerance
+                                       */
+      void tolerance (double);
+                                      /**
+                                       * Set strategy.
+                                       */
+      void strategy (Strategy);
+    
+                                      /**
+                                       * Set size of the first
+                                       * step. This may be overwritten
+                                       * by the time stepping strategy.
+                                       */
+      void start_step (double);
+    
+                                      /**
+                                       * Set size of the maximum
+                                       * step size.
+                                       */
+      void max_step (double);
+
+                                      /**
+                                       * Set now() equal to
+                                       * start(). Initialize step() and
+                                       * print() to their initial
+                                       * values.
+                                       */
+      void restart ();
+                                      /**
+                                       * Return true if this timestep
+                                       * should be written to disk.
+                                       */
+      bool print ();
+                                      /**
+                                       * Set the output name template.
+                                       */
+      void file_name_format (const char*);
+      const char* file_name_format ();
+    private:
+    
+      double start_val;
+      double final_val;
+      double tolerance_val;
+      Strategy strategy_val;
+      double start_step_val;
+      double max_step_val;
+      double min_step_val;
+                                      /**
+                                       * The size of the current time
+                                       * step. This may differ from
+                                       * #step_val, if we aimed at #final_val.
+                                       */
+      double current_step_val;
+      double step_val;
+    
+      double now_val;
+      double print_step;
+      double next_print_val;
+    
+      char format[30];
+  };
+
+
+  inline double
+  TimestepControl::start () const
+  {
+    return start_val;
+  }
+
+
+  inline double
+  TimestepControl::final () const
+  {
+    return final_val;
+  }
+
+
+  inline double
+  TimestepControl::step () const
+  {
+    return current_step_val;
+  }
+
+
+  inline double
+  TimestepControl::tolerance () const
+  {
+    return tolerance_val;
+  }
+
+
+  inline double
+  TimestepControl::now () const
+  {
+    return now_val;
+  }
+
+
+  inline void
+  TimestepControl::start (double t)
+  {
+    start_val = t;
+  }
+
+
+  inline void
+  TimestepControl::final (double t)
+  {
+    final_val = t;
+  }
+
+
+  inline void
+  TimestepControl::tolerance (double t)
+  {
+    tolerance_val = t;
+  }
+
+
+  inline void
+  TimestepControl::strategy (Strategy t)
+  {
+    strategy_val = t;
+  }
+
+
+  inline void
+  TimestepControl::start_step (double t)
+  {
+    start_step_val = t;
+  }
+
+
+  inline void
+  TimestepControl::max_step (double t)
+  {
+    max_step_val = t;
+  }
+
+
+  inline void
+  TimestepControl::restart ()
+  {
+    now_val = start_val;
+    step_val = start_step_val;
+    current_step_val = step_val;
+    if (print_step > 0.)
+      next_print_val = now_val + print_step;
+    else
+      next_print_val = now_val - 1.;
+  }
+
+
+  inline void
+  TimestepControl::file_name_format (const char* fmt)
+  {
+    strcpy(format, fmt);
+  }
+
+
+  inline const char*
+  TimestepControl::file_name_format ()
+  {
+    return format;
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
similarity index 79%
rename from deal.II/deal.II/source/numerics/operators.cc
rename to deal.II/deal.II/source/numerics/operators.all_dimensions.cc
index 9246625d715e85a4eb4355ae4c52a81ee2cfbd52..c42d9573a974ffc2afb51e20c73a63a473c725b9 100644 (file)
 
 #include <numerics/operator.templates.h>
 #include <numerics/newton.templates.h>
+#include <numerics/theta_timestepping.templates.h>
 
 #include <lac/vector.h>
 #include <lac/block_vector.h>
 
-#ifdef DEAL_II_USE_PETSC
-#  include <lac/petsc_vector.h>
-#  include <lac/petsc_block_vector.h>
-#endif
-
-#ifdef DEAL_II_USE_TRILINOS
-#  include <lac/trilinos_vector.h>
-#  include <lac/trilinos_block_vector.h>
-#endif
-
 DEAL_II_NAMESPACE_OPEN
 
 using namespace Algorithms;
index 75979cd3cbbb3870dc51c5f4f46203d9b5e4294a..7b5086980fd9d74323c5fbe6e94931f44b11c79f 100644 (file)
@@ -15,4 +15,5 @@ for (VEC : SERIAL_VECTORS)
 {  
   template class Operator<VEC>;
   template class Newton<VEC>;
+  template class ThetaTimestepping<VEC>;
 }
diff --git a/deal.II/deal.II/source/numerics/timestep_control.all_dimensions.cc b/deal.II/deal.II/source/numerics/timestep_control.all_dimensions.cc
new file mode 100644 (file)
index 0000000..fa05e6d
--- /dev/null
@@ -0,0 +1,137 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2006, 2007, 2010 by Guido Kanschat
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <numerics/time_step_control.h>
+#include <base/parameter_handler.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+using namespace Algorithms;
+
+TimestepControl::TimestepControl (double start,
+                                 double final,
+                                 double tolerance,
+                                 double start_step,
+                                 double print_step,
+                                 double max_step)
+  : start_val(start),
+    final_val(final),
+    tolerance_val(tolerance),
+    strategy_val(uniform),
+    start_step_val(start_step),
+    max_step_val(max_step),
+    min_step_val(0),
+    current_step_val(start_step),
+    step_val(start_step),
+    print_step(print_step)
+{
+  now_val = start_val;
+  strcpy(format, "T.%06.3f");
+}
+
+
+
+void TimestepControl::declare_parameters (ParameterHandler& param)
+{
+  param.declare_entry ("Start", "0.", Patterns::Double());
+  param.declare_entry ("Final", "1.", Patterns::Double());
+  param.declare_entry ("First step", "1.e-2", Patterns::Double());
+  param.declare_entry ("Max step", "1.", Patterns::Double());
+  param.declare_entry ("Tolerance", "1.e-2", Patterns::Double());
+  param.declare_entry ("Print step", "-1.", Patterns::Double());
+  param.declare_entry ("Strategy", "uniform",
+                      Patterns::Selection("uniform|doubling"));
+}
+
+
+
+
+void TimestepControl::parse_parameters (ParameterHandler& param)
+{
+  start (param.get_double ("Start"));
+  start_step (param.get_double ("First step"));
+  max_step (param.get_double ("Max step"));
+  final (param.get_double ("Final"));
+  tolerance (param.get_double ("Tolerance"));
+  print_step = param.get_double ("Print step");
+  const std::string strat = param.get("Strategy");
+  if (strat == std::string("uniform"))
+    strategy_val = uniform;
+  else if (strat == std::string("doubling"))
+    strategy_val = doubling;
+}
+
+
+
+
+bool
+TimestepControl::advance ()
+{
+  bool changed = false;
+  double s = step_val;
+
+                                  // Do time step control, but not in
+                                  // first step.
+  if (now_val != start())
+    {
+      if (strategy_val == doubling && 2*s <= tolerance_val)
+       s *= 2;
+      if (s > max_step_val)
+       s = max_step_val;
+    }
+
+                                  // Try incrementing time by s
+  double h = now_val + s;
+  changed = s != step_val;
+  
+  step_val = s;
+  current_step_val = s;
+                                  // If we just missed the final
+                                  // time, increase the step size a
+                                  // bit. This way, we avoid a very
+                                  // small final step. If the step
+                                  // shot over the final time, adjust
+                                  // it so we hit the final time
+                                  // exactly.
+  double s1 = .01*s;
+  if (h > final_val-s1)
+    {
+      current_step_val = final_val - now_val;
+      h = final_val;
+      changed = true;
+    }
+  
+  now_val = h;
+  return changed;
+}
+
+
+bool TimestepControl::print ()
+{
+  if (print_step == 0.)
+    return false;
+  if (print_step < 0.)
+    return true;
+  
+  bool result = (now_val >= next_print_val);
+
+  if (result)
+    {
+      next_print_val += print_step;
+      if (next_print_val > final_val)
+       next_print_val = final_val;
+    }
+  return result;
+}
+
+DEAL_II_NAMESPACE_CLOSE
+

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.