--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <base/config.h>
+
+#include <vector>
+#include <string>
+#include <iostream>
+
+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
+ *
+ * <pre>
+ * class A
+ * {
+ * static Event event;
+ * };
+ *
+ * Event A::event = Event::assign("Event for A");
+ * </pre>
+ */
+ 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 <tt>true</tt> if any
+ * event is set.
+ */
+ bool any () const;
+
+ /**
+ * List the flags to a stream.
+ */
+ template <class OS>
+ void print (OS& os) const;
+
+ /**
+ * List all assigned events.
+ */
+ template <class OS>
+ 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<bool> flags;
+
+ /**
+ * The names of registered events
+ */
+ static std::vector<std::string> 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<bool>::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<bool>::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<n_min;++i)
+ if (event.flags[i] && !flags[i])
+ return false;
+ for (unsigned int i=n_min;i<m;++i)
+ if (event.flags[i])
+ return false;
+ return true;
+ }
+
+
+
+ inline
+ Event& Event::operator += (const Event& event)
+ {
+ all_true |= event.all_true;
+ if (all_true) return *this;
+
+ if (flags.size() < event.flags.size())
+ flags.resize(event.flags.size());
+ for (unsigned int i=0;i<flags.size();++i)
+ flags[i] = flags[i] || event.flags[i];
+
+ return *this;
+ }
+
+
+ inline
+ Event& Event::operator -= (const Event& event)
+ {
+ if (!event.any()) return *this;
+
+ all_true = false;
+ if(event.all_true)
+ {
+ for (std::vector<bool>::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<flags.size();++i)
+ if (event.flags[i]) flags[i] = false;
+
+ return *this;
+ }
+
+
+ template <class OS>
+ inline
+ void
+ Event::print (OS& os) const
+ {
+ if (all_true)
+ os << " ALL";
+
+ for (unsigned int i=0;i<flags.size();++i)
+ if (flags[i])
+ os << ' ' << names[i];
+ }
+
+
+ template <class OS>
+ inline
+ void
+ Event::print_assigned (OS& os)
+ {
+ for (unsigned int i=0;i<names.size();++i)
+ os << i << '\t' << names[i] << std::endl;
+ }
+
+
+ /**
+ * Output shift operator for
+ * events. Calls Event::print().
+ *
+ * @relates Event
+ */
+ template <class OS>
+ OS& operator << (OS& o, const Event& e)
+ {
+ e.print(o);
+ return o;
+ }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
//---------------------------------------------------------------------------
// $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
//---------------------------------------------------------------------------
// $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
//---------------------------------------------------------------------------
// $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
//---------------------------------------------------------------------------
// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <base/smartpointer.h>
+#include <lac/solver_control.h>
+
+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 VECTOR>
+ class Newton : public Operator<VECTOR>
+ {
+ public:
+ /**
+ * Constructor, receiving the
+ * application computing the
+ * residual and solving the
+ * linear problem.
+ */
+ Newton (Operator<VECTOR>& residual, Operator<VECTOR>& 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<VECTOR*>& out, const NamedData<VECTOR*>& 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<Operator<VECTOR> > residual;
+
+ /**
+ * The operator applying the
+ * inverse derivative to the residual.
+ */
+ SmartPointer<Operator<VECTOR> > 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
+ * <tt>Stepsize iterations</tt> 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
+ * <tt>Newton_NNN</tt>?
+ */
+ bool debug_vectors;
+ /**
+ * Write debug output to
+ * #deallog; the higher the
+ * number, the more output.
+ */
+ unsigned int debug;
+ };
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <numerics/newton.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Algorithms
+{
+ template <class VECTOR>
+ Newton<VECTOR>::Newton(Operator<VECTOR>& residual, Operator<VECTOR>& inverse_derivative)
+ :
+ residual(&residual), inverse_derivative(&inverse_derivative),
+ assemble_now(false),
+ n_stepsize_iterations(21),
+ assemble_threshold(0.),
+ debug_vectors(false),
+ debug(0)
+ {}
+
+
+ template <class VECTOR>
+ void
+ Newton<VECTOR>::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 <class VECTOR>
+ void
+ Newton<VECTOR>::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 <class VECTOR>
+ void
+ Newton<VECTOR>::notify(const Event& e)
+ {
+ residual->notify(e);
+ inverse_derivative->notify(e);
+ }
+
+
+ template <class VECTOR>
+ void
+ Newton<VECTOR>::operator() (NamedData<VECTOR*>& out, const NamedData<VECTOR*>& in)
+ {
+ Assert (out.size() == 1, ExcNotImplemented());
+ deallog.push ("Newton");
+
+ VECTOR& u = *out(0);
+
+ if (debug>2)
+ deallog << "u: " << u.l2_norm() << std::endl;
+
+ GrowingVectorMemory<VECTOR> mem;
+ typename VectorMemory<VECTOR>::Pointer Du(mem);
+ typename VectorMemory<VECTOR>::Pointer res(mem);
+
+ res->reinit(u);
+ NamedData<VECTOR*> src1;
+ NamedData<VECTOR*> src2;
+ VECTOR* p = &u;
+ src1.add(p, "Newton iterate");
+ src1.merge(in);
+ p = res;
+ src2.add(p, "Newton residual");
+ src2.merge(src1);
+ NamedData<VECTOR*> out1;
+ out1.add(p, "Residual");
+ p = Du;
+ NamedData<VECTOR*> 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<<step_size)
+ << " since residual was " << residual << std::endl;
+ u.add(1./(1<<step_size), *Du);
+ (*residual)(out1, src1);
+ resnorm = res->l2_norm();
+ }
+
+ if (debug_vectors)
+ {
+ std::ostringstream streamOut;
+ streamOut << "Newton_" << std::setw(3) << std::setfill('0') << step;
+ std::string name = streamOut.str();
+ NamedData<VECTOR*> 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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <base/config.h>
+#include <base/named_data.h>
+#include <numerics/event.h>
+
+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 VECTOR>
+ class Operator : public Subscriptor
+ {
+ public:
+ /**
+ * The virtual destructor.
+ */
+ ~Operator();
+
+ /**
+ * The actual operation, which
+ * is implemented in a derived class.
+ */
+ virtual void operator() (NamedData<VECTOR*>& out, const NamedData<VECTOR*>& 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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <numerics/operator.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Algorithms
+{
+ template <class VECTOR>
+ Operator<VECTOR>::~Operator()
+ {}
+
+
+ template <class VECTOR>
+ void Operator<VECTOR>::notify(const Event& e)
+ {
+ notifications += e;
+ }
+
+
+ template <class VECTOR>
+ void
+ Operator<VECTOR>::clear_events ()
+ {
+ notifications.clear();
+ }
+}
+
+DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <numerics/event.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+//TODO: Thread safety
+
+namespace Algorithms
+{
+ std::vector<std::string> 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
//---------------------------------------------------------------------------
// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <base/parameter_handler.h>
+#include <base/logstream.h>
+#include <lac/vector_memory.h>
+
+#include <numerics/operator.templates.h>
+#include <numerics/newton.templates.h>
+
+#include <lac/vector.h>
+#include <lac/block_vector.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+using namespace Algorithms;
+
+#include "operators.inst"
+
+DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+//---------------------------------------------------------------------------
+// $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<VEC>;
+ template class Newton<VEC>;
+}