/// 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;
}
#include <base/smartpointer.h>
#include <lac/solver_control.h>
+#include <numerics/operator.h>
DEAL_II_NAMESPACE_OPEN
*/
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;
+ }
}
{
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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
#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;
{
template class Operator<VEC>;
template class Newton<VEC>;
+ template class ThetaTimestepping<VEC>;
}
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
+