From 6913e1f55ca6c96d0109070d30b76c64274b1676 Mon Sep 17 00:00:00 2001 From: kanschat Date: Mon, 10 Mar 2014 16:06:58 +0000 Subject: [PATCH] adjust ThetaTimestepping git-svn-id: https://svn.dealii.org/trunk@32639 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/algorithms/operator.h | 7 +- .../deal.II/algorithms/operator.templates.h | 36 ++++ .../deal.II/algorithms/theta_timestepping.h | 117 +++++------- .../algorithms/theta_timestepping.templates.h | 31 ++-- tests/algorithms/theta_01_compat.cc | 166 ++++++++++++++++++ tests/algorithms/theta_01_compat.output | 13 ++ 6 files changed, 287 insertions(+), 83 deletions(-) create mode 100644 tests/algorithms/theta_01_compat.cc create mode 100644 tests/algorithms/theta_01_compat.output diff --git a/deal.II/include/deal.II/algorithms/operator.h b/deal.II/include/deal.II/algorithms/operator.h index 0bc6575e77..6254d35019 100644 --- a/deal.II/include/deal.II/algorithms/operator.h +++ b/deal.II/include/deal.II/algorithms/operator.h @@ -196,7 +196,12 @@ namespace Algorithms OutputOperator &operator<< (unsigned int step); /** - * Output all the vectors in NamedData. + * Output all the vectors in AnyData. + */ + virtual OutputOperator &operator<< (const AnyData& vectors); + + /** + * @deprecated Output all the vectors in NamedData. */ virtual OutputOperator &operator<< (const NamedData &vectors); protected: diff --git a/deal.II/include/deal.II/algorithms/operator.templates.h b/deal.II/include/deal.II/algorithms/operator.templates.h index e2a7458569..5b16380492 100644 --- a/deal.II/include/deal.II/algorithms/operator.templates.h +++ b/deal.II/include/deal.II/algorithms/operator.templates.h @@ -119,6 +119,42 @@ namespace Algorithms os =&stream; } + template + OutputOperator & + OutputOperator::operator<< (const AnyData& vectors) + { + if (os == 0) + { + //TODO: make this possible + //deallog << ' ' << step; + //for (unsigned int i=0;iprint(deallog); + //deallog << std::endl; + } + else + { + (*os) << ' ' << step; + for (unsigned int i=0; i(i)) + { + const VECTOR& v = *vectors.entry(i); + for (unsigned int j=0; j(i)) + { + const VECTOR& v = *vectors.entry(i); + for (unsigned int j=0; j OutputOperator & OutputOperator::operator<< (const NamedData &vectors) diff --git a/deal.II/include/deal.II/algorithms/theta_timestepping.h b/deal.II/include/deal.II/algorithms/theta_timestepping.h index e1a2d0c40d..43ade372df 100644 --- a/deal.II/include/deal.II/algorithms/theta_timestepping.h +++ b/deal.II/include/deal.II/algorithms/theta_timestepping.h @@ -188,6 +188,11 @@ namespace Algorithms * should contain the initial * value in first position. out */ + virtual void operator() (AnyData &out, const AnyData &in); + + /** + * @deprecated Use the function with AnyData arguments + */ virtual void operator() (NamedData &out, const NamedData &in); /** @@ -197,11 +202,8 @@ namespace Algorithms 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. + * 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); @@ -217,34 +219,23 @@ namespace Algorithms */ const double &step_size() const; /** - * The weight between implicit - * and explicit part. + * The weight between implicit and explicit part. */ const double &theta() const; /** - * The data handed to the - * #op_explicit time stepping - * operator. + * The data handed to the #op_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. + * 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 - * #op_implicit time stepping - * operator. + * The data handed to the #op_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. + * 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; @@ -255,84 +246,61 @@ namespace Algorithms private: /** - * The object controlling the - * time step size and computing - * the new time in each step. + * 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]. + * The control parameter theta in the range [0,1]. It + * defaults to 0.5. */ double vtheta; /** - * Use adaptive #theta if - * true. + * Use adaptive #theta if true. Not yet implemented. */ bool adaptive; /** - * The data for the explicit - * part of the scheme. + * The data for the explicit part of the scheme. */ TimestepData d_explicit; /** - * The data for the implicit - * part of the scheme. + * 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(). + * 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(). + * 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, ThetaTimestepping > op_explicit; /** - * The operator solving the - * implicit part of the - * scheme. It will receive in - * its input data the vector - * "Previous time". Information on the - * timestep should be obtained - * from implicit_data(). + * The operator solving the implicit part of the scheme. It will + * receive in its input data the vector "Previous + * time". 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 + * 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, ThetaTimestepping > op_implicit; /** - * The operator writing the - * output in each time step + * The operator writing the output in each time step */ SmartPointer, ThetaTimestepping > output; }; @@ -355,6 +323,15 @@ namespace Algorithms return d_implicit; } + + template + inline + TimestepControl & + ThetaTimestepping::timestep_control () + { + return control; + } + template inline void ThetaTimestepping::set_output (OutputOperator &out) diff --git a/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h b/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h index b0f1366641..6194088cfe 100644 --- a/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h +++ b/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h @@ -64,13 +64,23 @@ namespace Algorithms template void ThetaTimestepping::operator() (NamedData &out, const NamedData &in) + { + Operator::operator() (out, in); + } + + + template + void + ThetaTimestepping::operator() (AnyData &out, const AnyData &in) { Assert(!adaptive, ExcNotImplemented()); deallog.push ("Theta"); + + VECTOR& solution = *out.entry(0); GrowingVectorMemory mem; typename VectorMemory::Pointer aux(mem); - aux->reinit(*out(0)); + aux->reinit(solution); control.restart(); @@ -79,20 +89,17 @@ namespace Algorithms // The data used to compute the // vector associated with the old // timestep - VECTOR *p = out(0); - NamedData src1; - src1.add(p, "Previous iterate"); + AnyData src1; + src1.add(&solution, "Previous iterate"); src1.merge(in); - NamedData src2; - src2.add(p, "Previous iterate"); + AnyData src2; + src2.add(&solution, "Previous iterate"); - p = aux; - NamedData out1; - out1.add(p, "Result"); - // The data provided to the inner - // solver - src2.add(p, "Previous time"); + AnyData out1; + out1.add(aux, "Solution"); + // The data provided to the inner solver + src2.add(aux, "Previous time"); src2.merge(in); if (output != 0) diff --git a/tests/algorithms/theta_01_compat.cc b/tests/algorithms/theta_01_compat.cc new file mode 100644 index 0000000000..84041252a2 --- /dev/null +++ b/tests/algorithms/theta_01_compat.cc @@ -0,0 +1,166 @@ +// --------------------------------------------------------------------- +// $Id: theta_timestepping.cc 30050 2013-07-18 20:03:31Z maier $ +// +// Copyright (C) 2005 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// See documentation of ThetaTimestepping for documentation of this example + +#include +#include +#include + +#include +#include + +#include + +using namespace dealii; +using namespace Algorithms; + + +class Explicit + : public Operator > +{ +public: + Explicit(const FullMatrix &matrix); + void operator() (NamedData*> &out, + const NamedData*> &in); + + void initialize_timestep_data(const TimestepData &); +private: + const TimestepData *timestep_data; + SmartPointer, Explicit> matrix; + FullMatrix m; +}; + + +class Implicit + : public Operator > +{ +public: + Implicit(const FullMatrix &matrix); + void operator() (NamedData*> &out, + const NamedData*> &in); + + void initialize_timestep_data(const TimestepData &); +private: + const TimestepData *timestep_data; + SmartPointer, Implicit> matrix; + FullMatrix m; +}; + +// End of declarations + +int main() +{ + std::string logname = "output"; + std::ofstream logfile(logname.c_str()); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + FullMatrix matrix(2); + matrix(0,0) = 0.; + matrix(1,1) = 0.; + matrix(0,1) = numbers::PI; + matrix(1,0) = -numbers::PI; + + OutputOperator > out; + out.initialize_stream(logfile); + + Explicit op_explicit(matrix); + Implicit op_implicit(matrix); + ThetaTimestepping > solver(op_explicit, op_implicit); + op_explicit.initialize_timestep_data(solver.explicit_data()); + op_implicit.initialize_timestep_data(solver.implicit_data()); + solver.timestep_control().start_step(0.1); + //solver.set_output(out); + + Vector value(2); + value(0) = 1.; + NamedData*> indata; + NamedData*> outdata; + Vector *p = &value; + outdata.add(p, "value"); + deallog << "Initial: " << value(0) << ' ' << value(1) << std::endl; + solver.notify(Events::initial); + solver(outdata, indata); + deallog << "Result: " << value(0) << ' ' << value(1) + << " Norm " << value.l2_norm() << std::endl; +} + + +Explicit::Explicit(const FullMatrix &M) + : + matrix(&M) +{ + m.reinit(M.m(), M.n()); +} + + +void +Explicit::initialize_timestep_data(const TimestepData &t) +{ + timestep_data = &t; +} + + +void +Explicit::operator() (NamedData*> &out, const NamedData*> &in) +{ + if (this->notifications.test(Events::initial) || this->notifications.test(Events::new_timestep_size)) + { + m.equ(-timestep_data->step, *matrix); + for (unsigned int i=0; inotifications.clear(); + unsigned int i = in.find("Previous iterate"); + m.vmult(*out(0), *in(i)); +} + + +Implicit::Implicit(const FullMatrix &M) + : + matrix(&M) +{ + m.reinit(M.m(), M.n()); +} + + +void +Implicit::initialize_timestep_data(const TimestepData &t) +{ + timestep_data = &t; +} + + +void +Implicit::operator() (NamedData*> &out, const NamedData*> &in) +{ + if (this->notifications.test(Events::initial) || this->notifications.test(Events::new_timestep_size)) + { + m.equ(timestep_data->step, *matrix); + for (unsigned int i=0; inotifications.clear(); + + unsigned int i = in.find("Previous time"); + m.vmult(*out(0), *in(i)); +} + + diff --git a/tests/algorithms/theta_01_compat.output b/tests/algorithms/theta_01_compat.output new file mode 100644 index 0000000000..7c4d9ddfc7 --- /dev/null +++ b/tests/algorithms/theta_01_compat.output @@ -0,0 +1,13 @@ + +DEAL::Initial: 1.00000 0 +DEAL:Theta::Time step:0.100000 +DEAL:Theta::Time step:0.200000 +DEAL:Theta::Time step:0.300000 +DEAL:Theta::Time step:0.400000 +DEAL:Theta::Time step:0.500000 +DEAL:Theta::Time step:0.600000 +DEAL:Theta::Time step:0.700000 +DEAL:Theta::Time step:0.800000 +DEAL:Theta::Time step:0.900000 +DEAL:Theta::Time step:1.00000 +DEAL::Result: -0.999676 0.0254599 Norm 1.00000 -- 2.39.5