From 2b6109b61b8210d94d9e6072c730f4ff4e70a837 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Wed, 13 Jan 2010 15:21:48 +0000 Subject: [PATCH] ThetaTimestepping implemented and ready for testing git-svn-id: https://svn.dealii.org/trunk@20359 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/event.h | 4 +- deal.II/deal.II/include/numerics/newton.h | 1 + deal.II/deal.II/include/numerics/operator.h | 34 ++ .../include/numerics/operator.templates.h | 20 ++ .../include/numerics/theta_timestepping.h | 254 ++++++++++++++ .../numerics/theta_timestepping.templates.h | 123 +++++++ .../include/numerics/time_step_control.h | 323 ++++++++++++++++++ ...erators.cc => operators.all_dimensions.cc} | 11 +- .../deal.II/source/numerics/operators.inst.in | 1 + .../timestep_control.all_dimensions.cc | 137 ++++++++ 10 files changed, 897 insertions(+), 11 deletions(-) create mode 100644 deal.II/deal.II/include/numerics/theta_timestepping.h create mode 100644 deal.II/deal.II/include/numerics/theta_timestepping.templates.h create mode 100644 deal.II/deal.II/include/numerics/time_step_control.h rename deal.II/deal.II/source/numerics/{operators.cc => operators.all_dimensions.cc} (79%) create mode 100644 deal.II/deal.II/source/numerics/timestep_control.all_dimensions.cc diff --git a/deal.II/deal.II/include/numerics/event.h b/deal.II/deal.II/include/numerics/event.h index 22a8366216..84f5bd51da 100644 --- a/deal.II/deal.II/include/numerics/event.h +++ b/deal.II/deal.II/include/numerics/event.h @@ -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; } diff --git a/deal.II/deal.II/include/numerics/newton.h b/deal.II/deal.II/include/numerics/newton.h index a6882c739e..97ba038652 100644 --- a/deal.II/deal.II/include/numerics/newton.h +++ b/deal.II/deal.II/include/numerics/newton.h @@ -15,6 +15,7 @@ #include #include +#include DEAL_II_NAMESPACE_OPEN diff --git a/deal.II/deal.II/include/numerics/operator.h b/deal.II/deal.II/include/numerics/operator.h index 2f1c799af1..a43b6acdbc 100644 --- a/deal.II/deal.II/include/numerics/operator.h +++ b/deal.II/deal.II/include/numerics/operator.h @@ -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 OutputOperator : public Subscriptor + { + public: + /** + * Empty virtual destructor. + */ + virtual ~OutputOperator(); + /** + * Set the current step. + */ + OutputOperator& operator<< (unsigned int step); + /** + * Output all the vectors in NamedData. + */ + virtual OutputOperator& operator<< (const NamedData vectors); + protected: + unsigned int step; + }; + + template + OutputOperator& + OutputOperator::operator<< (unsigned int s) + { + step = s; + return *this; + } } diff --git a/deal.II/deal.II/include/numerics/operator.templates.h b/deal.II/deal.II/include/numerics/operator.templates.h index 48c9ddac46..20b7bc96da 100644 --- a/deal.II/deal.II/include/numerics/operator.templates.h +++ b/deal.II/deal.II/include/numerics/operator.templates.h @@ -34,6 +34,26 @@ namespace Algorithms { notifications.clear(); } + + + template + OutputOperator::~OutputOperator() + {} + + + template + OutputOperator& + OutputOperator::operator<< (const NamedData vectors) + { + std::cout << ' ' << step; + for (unsigned int i=0;iprint(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 index 0000000000..8386033983 --- /dev/null +++ b/deal.II/deal.II/include/numerics/theta_timestepping.h @@ -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 +#include +#include + +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. + *
    + *
  • #theta=0: explicit Euler scheme + *
  • #theta=1: implicit Euler scheme + *
  • #theta=½: Crank-Nicholson scheme + *
+ * + * 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 #theta=½+c dt. + * + * Assume that we want to solve the equation u' = Au with a + * step size k. 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, M 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 ThetaTimestepping : public Operator + { + public: + /** + * Constructor, receiving the + * two operators stored in + * #op_explicit and #op_implicit. For + * their meening, see the + * description of those variables. + */ + ThetaTimestepping (Operator& op_explicit, + Operator& op_implicit); + + /** + * The timestepping scheme. in + * should contain the initial + * value in first position. out + */ + virtual void operator() (NamedData& out, const NamedData& 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& 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 [0,1]. + */ + double vtheta; + /** + * Use adaptive #theta if + * true. + */ + 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 + * Mu+cAu, where + * u is the current + * state vector, M the + * mass matrix, A the + * operator in space and + * c is the time step + * size in explicit_data(). + */ + SmartPointer > 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 u of + * Mu-cAu=f, where + * f is the dual space + * vector found in the + * "Previous time" entry of the + * input data, M the + * mass matrix, A the + * operator in space and + * c is the time step + * size in explicit_data(). + */ + SmartPointer > op_implicit; + + /** + * The operator writing the + * output in each time step + */ + SmartPointer > 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 index 0000000000..2d324b80b8 --- /dev/null +++ b/deal.II/deal.II/include/numerics/theta_timestepping.templates.h @@ -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 +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Algorithms +{ + template + ThetaTimestepping::ThetaTimestepping (Operator& e, Operator& i) + : op_explicit(&e), op_implicit(&i) + {} + + + template + void + ThetaTimestepping::notify(const Event& e) + { + op_explicit->notify(e); + op_implicit->notify(e); + } + + template + void + ThetaTimestepping::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 + void + ThetaTimestepping::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 + void + ThetaTimestepping::operator() (NamedData& out, const NamedData& in) + { + Assert(!adaptive, ExcNotImplemented()); + + deallog.push ("Theta"); + GrowingVectorMemory mem; + typename VectorMemory::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 src1; + src1.add(p, "Previous time"); + src1.merge(in); + p = aux; + NamedData out1; + out1.add(p, "Result"); + // The data provided to the inner + // solver + NamedData 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 index 0000000000..4646703144 --- /dev/null +++ b/deal.II/deal.II/include/numerics/time_step_control.h @@ -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 +#include +#include +#include + +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: + *
    + *
  1. According to the strategy, the step size is tentatively added to the + * current time. + *
  2. If the resulting time exceeds the final time of the interval, + * the step size is reduced in order to meet this time. + *
  3. 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. + *
  4. The resulting step size is used from the current time. + *
+ * + * 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 diff --git a/deal.II/deal.II/source/numerics/operators.cc b/deal.II/deal.II/source/numerics/operators.all_dimensions.cc 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 9246625d71..c42d9573a9 100644 --- a/deal.II/deal.II/source/numerics/operators.cc +++ b/deal.II/deal.II/source/numerics/operators.all_dimensions.cc @@ -17,20 +17,11 @@ #include #include +#include #include #include -#ifdef DEAL_II_USE_PETSC -# include -# include -#endif - -#ifdef DEAL_II_USE_TRILINOS -# include -# include -#endif - DEAL_II_NAMESPACE_OPEN using namespace Algorithms; diff --git a/deal.II/deal.II/source/numerics/operators.inst.in b/deal.II/deal.II/source/numerics/operators.inst.in index 75979cd3cb..7b5086980f 100644 --- a/deal.II/deal.II/source/numerics/operators.inst.in +++ b/deal.II/deal.II/source/numerics/operators.inst.in @@ -15,4 +15,5 @@ for (VEC : SERIAL_VECTORS) { template class Operator; template class Newton; + template class ThetaTimestepping; } 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 index 0000000000..fa05e6dbe7 --- /dev/null +++ b/deal.II/deal.II/source/numerics/timestep_control.all_dimensions.cc @@ -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 +#include + +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 + -- 2.39.5