From: Guido Kanschat Date: Tue, 12 Jan 2010 16:56:38 +0000 (+0000) Subject: start algorithms library X-Git-Tag: v8.0.0~6642 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e6eff760005d946e95c6ae1b09abc63d2f86238e;p=dealii.git start algorithms library git-svn-id: https://svn.dealii.org/trunk@20350 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/event.h b/deal.II/deal.II/include/numerics/event.h new file mode 100644 index 0000000000..7074e0065b --- /dev/null +++ b/deal.II/deal.II/include/numerics/event.h @@ -0,0 +1,287 @@ +//--------------------------------------------------------------------------- +// $Id: timestepping.cc 822 2010-01-11 08:57:30Z kanschat $ +// +// 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__event_h +#define __deal2__event_h + +#include + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Algorithms +{ +/** + * Objects of this kind are used to notify interior applications of + * changes provoked by an outer loop. They are handed to the + * application through Operator::notify() and it is up to the + * actual application how to handle them. + * + * Event is organized as an extensible binary enumerator. Every class + * can add its own events using assign(). A typical code example is + * + *
+ * class A
+ * {
+ *   static Event event;
+ * };
+ *
+ * Event A::event = Event::assign("Event for A");
+ * 
+ */ + class Event + { + public: + /** + * This function registers a + * new event type and assigns a + * unique identifier to it. The + * result of this function + * should be stored for later + * use. + */ + static Event assign (const char* name); + + /** + * If you forgot to store the + * result of assign, here is + * how to retrieve it knowing + * the name. + */ +// static Event find(const std::string& name); + + /** + * Constructor, generating a + * clear Event. + */ + Event (); + + /** + * Clear all flags + */ + void clear(); + + /** + * Set all flags + */ + void all(); + + /** + * Add the flags of the other event + */ + Event& operator += (const Event& event); + + /** + * Clear the flags of the other event + */ + Event& operator -= (const Event& event); + + /** + * Test whether all the flags + * set in the other Event are + * also set in this one. + */ + bool test (const Event& event) const; + + /** + * Return true if any + * event is set. + */ + bool any () const; + + /** + * List the flags to a stream. + */ + template + void print (OS& os) const; + + /** + * List all assigned events. + */ + template + static void print_assigned (OS& os); + + private: + /** + * Sometimes, actions have to + * be taken by all + * means. Therefore, if this + * value is true, test() always + * returns true. + */ + bool all_true; + + /** + * The actual list of events + */ + std::vector flags; + + /** + * The names of registered events + */ + static std::vector names; + }; + +/// Events used by library operators + namespace Events + { +/// 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; + } + + +//----------------------------------------------------------------------// + + + inline + bool + Event::any () const + { + if (all_true) return true; + for (std::vector::const_iterator i=flags.begin(); + i != flags.end(); ++i) + if (*i) return true; + return false; + } + + + inline + bool + Event::test (const Event& event) const + { + + // First, test all_true in this + if (all_true) return true; + + const unsigned int n = flags.size(); + const unsigned int m = event.flags.size(); + const unsigned int n_min = std::min(n, m); + + // Now, if all_true set in the + // other, then all must be true + // in this + if (event.all_true) + { + // Non existing flags are + // always assumed false + if (m > n) + return false; + + // Test all flags separately + // and return false if one is + // not set + for (std::vector::const_iterator i=flags.begin(); + i != flags.end(); ++i) + if (!*i) return false; + // All flags are set + return true; + } + + // Finally, compare each flag + // separately + for (unsigned int i=0;i::iterator i=flags.begin(); + i != flags.end(); ++i) + *i = false; + return *this; + } + + if (flags.size() < event.flags.size()) + flags.resize(event.flags.size()); + for (unsigned int i=0;i + inline + void + Event::print (OS& os) const + { + if (all_true) + os << " ALL"; + + for (unsigned int i=0;i + inline + void + Event::print_assigned (OS& os) + { + for (unsigned int i=0;i + OS& operator << (OS& o, const Event& e) + { + e.print(o); + return o; + } +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/deal.II/include/numerics/mesh_worker.h b/deal.II/deal.II/include/numerics/mesh_worker.h index 5f956c2f3f..535a2216a9 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker.h +++ b/deal.II/deal.II/include/numerics/mesh_worker.h @@ -1,7 +1,7 @@ //--------------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2006, 2007, 2008, 2009, 2010 by Guido Kanschat +// Copyright (C) 2006, 2007, 2008, 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer diff --git a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h index c822c195bb..f330e125c7 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h @@ -1,7 +1,7 @@ //--------------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2006, 2007, 2008, 2009, 2010 by Guido Kanschat +// Copyright (C) 2006, 2007, 2008, 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer diff --git a/deal.II/deal.II/include/numerics/mesh_worker_info.h b/deal.II/deal.II/include/numerics/mesh_worker_info.h index 795410bece..853c7496e2 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_info.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_info.h @@ -1,7 +1,7 @@ //--------------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2006, 2007, 2008, 2009, 2010 by Guido Kanschat +// Copyright (C) 2006, 2007, 2008, 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer diff --git a/deal.II/deal.II/include/numerics/mesh_worker_loop.h b/deal.II/deal.II/include/numerics/mesh_worker_loop.h index 2d79b8f10a..909af2f9bd 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_loop.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_loop.h @@ -1,7 +1,7 @@ //--------------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2006, 2007, 2008, 2009, 2010 by Guido Kanschat +// Copyright (C) 2006, 2007, 2008, 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer diff --git a/deal.II/deal.II/include/numerics/newton.h b/deal.II/deal.II/include/numerics/newton.h new file mode 100644 index 0000000000..5bb44dea33 --- /dev/null +++ b/deal.II/deal.II/include/numerics/newton.h @@ -0,0 +1,161 @@ +//--------------------------------------------------------------------------- +// $Id: timestepping.cc 822 2010-01-11 08:57:30Z kanschat $ +// +// 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__newton_step_control_h +#define __deal2__newton_step_control_h + +#include +#include + +DEAL_II_NAMESPACE_OPEN + +class ParameterHandler; + +namespace Algorithms +{ +/** + * Operator class performing Newton's iteration with standard step + * size control and adaptive matrix generation. + * + * This class performes a Newton iteration up to convergence + * determined by #control. If after an update the norm of the residual + * has become larger, then step size control is activated and the + * update is subsequently divided by two until the residual actually + * becomes smaller (or the minimal scaling factor determined by + * #n_stepsize_iterations is reached). + * + * Since assembling matrices, depending on the implementation, tends + * to be costly, this method applies an adaptive reassembling + * strategy. Only if the reduction factor for the residual is more + * than #threshold, the event #bad_derivative is submitted to + * #inverse_derivative. It is up to this object to implement + * reassembling accordingly. + * + * @author Guido Kanschat, 2006, 2010 + */ + template + class Newton : public Operator + { + public: + /** + * Constructor, receiving the + * application computing the + * residual and solving the + * linear problem. + */ + Newton (Operator& residual, Operator& inverse_derivative); + + /** + * Declare the parameter + * applicable to Newton's method. + */ + void declare_parameters (ParameterHandler& param); + + /** + * Read the parameters. + */ + void initialize (ParameterHandler& param); + + /** + * The actual Newton iteration. + */ + virtual void operator() (NamedData& out, const NamedData& in); + + virtual void notify(const Event&); + + /** + * Set the maximal residual + * reduction allowed without + * triggering assembling in the + * next step. Return the + * previous value. + */ + double threshold(double new_value); + + /** + * Control object for the + * Newton iteration. + */ + ReductionControl control; + private: + /** + * The operator computing the residual. + */ + SmartPointer > residual; + + /** + * The operator applying the + * inverse derivative to the residual. + */ + SmartPointer > inverse_derivative; + + /** + * This flag is set by the + * function assemble(), + * indicating that the matrix + * must be assembled anew upon + * start. + */ + bool assemble_now; + + /** + * A flag used to decide how + * many stepsize iteration + * should be made. Default is + * the original value of 21. + * + * Enter zero here to turn of + * stepsize control. + * + * @note Controlled by + * Stepsize iterations in + * parameter file + */ + unsigned int n_stepsize_iterations; + + /** + * Threshold for re-assembling matrix. + * + * If the quotient of two + * consecutive residuals is + * smaller than this threshold, + * the system matrix is not + * assembled in this step. + * + * @note This parameter should be + * adjusted to the residual gain + * of the inner solver. + * + * The default values is zero, + * resulting in reassembling in + * every Newton step. + */ + double assemble_threshold; + /** + * Print residual, update and + * updated solution after each + * step into file + * Newton_NNN? + */ + bool debug_vectors; + /** + * Write debug output to + * #deallog; the higher the + * number, the more output. + */ + unsigned int debug; + }; +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/deal.II/include/numerics/newton.templates.h b/deal.II/deal.II/include/numerics/newton.templates.h new file mode 100644 index 0000000000..1918473341 --- /dev/null +++ b/deal.II/deal.II/include/numerics/newton.templates.h @@ -0,0 +1,169 @@ +//--------------------------------------------------------------------------- +// $Id: newton.h 753 2009-05-27 00:30:14Z kanschat $ +// +// 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 + +DEAL_II_NAMESPACE_OPEN + +namespace Algorithms +{ + template + Newton::Newton(Operator& residual, Operator& inverse_derivative) + : + residual(&residual), inverse_derivative(&inverse_derivative), + assemble_now(false), + n_stepsize_iterations(21), + assemble_threshold(0.), + debug_vectors(false), + debug(0) + {} + + + template + void + Newton::declare_parameters(ParameterHandler& param) + { + param.enter_subsection("Newton"); + ReductionControl::declare_parameters (param); + param.declare_entry("Assemble threshold", "0.", Patterns::Double()); + param.declare_entry("Stepsize iterations", "21", Patterns::Integer()); + param.declare_entry("Debug level", "0", Patterns::Integer()); + param.declare_entry("Debug vectors", "false", Patterns::Bool()); + param.leave_subsection(); + } + + template + void + Newton::initialize (ParameterHandler& param) + { + param.enter_subsection("Newton"); + control.parse_parameters (param); + assemble_threshold = param.get_double("Assemble threshold"); + n_stepsize_iterations = param.get_double("Stepsize iterations"); + debug_vectors = param.get_bool("Debug vectors"); + param.leave_subsection (); + } + + + template + void + Newton::notify(const Event& e) + { + residual->notify(e); + inverse_derivative->notify(e); + } + + + template + void + Newton::operator() (NamedData& out, const NamedData& in) + { + Assert (out.size() == 1, ExcNotImplemented()); + deallog.push ("Newton"); + + VECTOR& u = *out(0); + + if (debug>2) + deallog << "u: " << u.l2_norm() << std::endl; + + GrowingVectorMemory mem; + typename VectorMemory::Pointer Du(mem); + typename VectorMemory::Pointer res(mem); + + res->reinit(u); + NamedData src1; + NamedData src2; + VECTOR* p = &u; + src1.add(p, "Newton iterate"); + src1.merge(in); + p = res; + src2.add(p, "Newton residual"); + src2.merge(src1); + NamedData out1; + out1.add(p, "Residual"); + p = Du; + NamedData out2; + out2.add(p, "Update"); + + unsigned int step = 0; + // fill res with (f(u), v) + (*residual)(out1, src1); + double resnorm = res->l2_norm(); + double old_residual = resnorm / assemble_threshold + 1; + + while (control.check(step++, resnorm) == SolverControl::iterate) + { + // assemble (Df(u), v) + if (resnorm/old_residual >= assemble_threshold) + inverse_derivative->notify (Events::bad_derivative); + + Du->reinit(u); + try + { + (*inverse_derivative)(out2, src2); + } + catch (SolverControl::NoConvergence& e) + { + deallog << "Inner iteration failed after " + << e.last_step << " steps with residual " + << e.last_residual << std::endl; + } + + u.add(-1., *Du); + old_residual = resnorm; + (*residual)(out1, src1); + resnorm = res->l2_norm(); + + // Step size control + unsigned int step_size = 0; + while (resnorm >= old_residual) + { + ++step_size; + if (step_size > n_stepsize_iterations) + { + deallog << "No smaller stepsize allowed!"; + break; + } + if (control.log_history()) + deallog << "Trying step size: 1/" << (1<l2_norm(); + } + + if (debug_vectors) + { + std::ostringstream streamOut; + streamOut << "Newton_" << std::setw(3) << std::setfill('0') << step; + std::string name = streamOut.str(); + NamedData out; + out.add(&u, "solution"); + out.add(Du, "update"); + out.add(res, "residual"); +//TODO: Implement! +// app->data_out(name, out); + } + } + deallog.pop(); + + // in case of failure: throw + // exception + if (control.last_check() != SolverControl::success) + throw SolverControl::NoConvergence (control.last_step(), + control.last_value()); + // otherwise exit as normal + } +} + + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/include/numerics/operator.h b/deal.II/deal.II/include/numerics/operator.h new file mode 100644 index 0000000000..e6d70f7211 --- /dev/null +++ b/deal.II/deal.II/include/numerics/operator.h @@ -0,0 +1,79 @@ +//--------------------------------------------------------------------------- +// $Id: timestepping.cc 822 2010-01-11 08:57:30Z kanschat $ +// +// 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__operator_h +#define __deal2__operator_h + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Algorithms +{ +/** + * The abstract base class of all algorithms in this library. An + * operator is an object with an #operator(), which transforms a set + * of named vectors into another set of named vectors. + * + * Furthermore, an operator can be notified of parameter changes by + * the calling routine. The outer iteration can notify() the Operator + * of an Event, which could be for instance a change of mesh, a + * different time step size or too slow convergence of Newton's + * method, which would then trigger reassembling of a matrix or + * similar things. + * + * @author Guido Kanschat, 2010 + */ + template + class Operator : public Subscriptor + { + public: + /** + * The virtual destructor. + */ + ~Operator(); + + /** + * The actual operation, which + * is implemented in a derived class. + */ + virtual void operator() (NamedData& out, const NamedData& in) = 0; + + /** + * Register an event triggered + * by an outer iteration. + */ + virtual void notify(const Event&); + /** + * Clear all #notifications. + */ + void clear_events(); + protected: + /** + * Accumulate reasons for + * reassembling here. If any of + * those is set, the function + * solve() of a terminal + * application must take care + * of reassembling the matrix. + */ + Event notifications; + }; + +} + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/deal.II/include/numerics/operator.templates.h b/deal.II/deal.II/include/numerics/operator.templates.h new file mode 100644 index 0000000000..fc9095bd17 --- /dev/null +++ b/deal.II/deal.II/include/numerics/operator.templates.h @@ -0,0 +1,39 @@ +//--------------------------------------------------------------------------- +// $Id: newton.h 753 2009-05-27 00:30:14Z kanschat $ +// +// 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 + +DEAL_II_NAMESPACE_OPEN + +namespace Algorithms +{ + template + Operator::~Operator() + {} + + + template + void Operator::notify(const Event& e) + { + notifications += e; + } + + + template + void + Operator::clear_events () + { + notifications.clear(); + } +} + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/numerics/event.all_dimensions.cc b/deal.II/deal.II/source/numerics/event.all_dimensions.cc new file mode 100644 index 0000000000..ad45b3aac2 --- /dev/null +++ b/deal.II/deal.II/source/numerics/event.all_dimensions.cc @@ -0,0 +1,69 @@ +//--------------------------------------------------------------------------- +// $Id: vectors.all_dimensions.cc 14038 2006-10-23 02:46:34Z bangerth $ +// Version: $Name$ +// +// 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. +// +//--------------------------------------------------------------------------- + +#include + +DEAL_II_NAMESPACE_OPEN + +//TODO: Thread safety + +namespace Algorithms +{ + std::vector Event::names; + + Event + Event::assign(const char* name) + { + unsigned int index = names.size(); + names.push_back(name); + + Event result; + // The constructor generated an + // object with all flags equal + // zero. Now we set the new one. + result.flags[index] = true; + + return result; + } + + + Event::Event () + : + all_true(false), + flags(names.size(), false) + {} + + + void + Event::clear () + { + all_true = false; + std::fill(flags.begin(), flags.end(), false); + } + + + void + Event::all () + { + all_true = true; + } + + namespace Events + { + const Event bad_derivative = Event::assign("Bad Derivative"); + const Event new_time = Event::assign("New Time"); + const Event new_time_step_size = Event::assign("New Time Step Size"); + } +} + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/numerics/mesh_worker.cc b/deal.II/deal.II/source/numerics/mesh_worker.cc index a89e80acc9..1ccad5b062 100644 --- a/deal.II/deal.II/source/numerics/mesh_worker.cc +++ b/deal.II/deal.II/source/numerics/mesh_worker.cc @@ -1,7 +1,7 @@ //--------------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2006, 2007, 2008, 2009, 2010 by Guido Kanschat +// Copyright (C) 2006, 2007, 2008, 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer diff --git a/deal.II/deal.II/source/numerics/operators.cc b/deal.II/deal.II/source/numerics/operators.cc new file mode 100644 index 0000000000..0e09fb4f5b --- /dev/null +++ b/deal.II/deal.II/source/numerics/operators.cc @@ -0,0 +1,30 @@ +//--------------------------------------------------------------------------- +// $Id: vectors.cc 18474 2009-03-10 16:55:43Z kronbichler $ +// Version: $Name$ +// +// 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. +// +//--------------------------------------------------------------------------- + +#include +#include +#include + +#include +#include + +#include +#include + +DEAL_II_NAMESPACE_OPEN + +using namespace Algorithms; + +#include "operators.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/numerics/operators.inst.in b/deal.II/deal.II/source/numerics/operators.inst.in new file mode 100644 index 0000000000..75979cd3cb --- /dev/null +++ b/deal.II/deal.II/source/numerics/operators.inst.in @@ -0,0 +1,18 @@ +//--------------------------------------------------------------------------- +// $Id: vectors.inst.in 19046 2009-07-08 19:30:23Z bangerth $ +// Version: $Name$ +// +// 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. +// +//--------------------------------------------------------------------------- + +for (VEC : SERIAL_VECTORS) +{ + template class Operator; + template class Newton; +}