]> https://gitweb.dealii.org/ - dealii.git/commitdiff
updated documentation
authorvishalkenchan <vishal.boddu@fau.de>
Wed, 12 Jul 2017 08:30:53 +0000 (10:30 +0200)
committerVishal Boddu <vishal.boddu@fau.de>
Fri, 14 Jul 2017 22:19:11 +0000 (00:19 +0200)
doc/news/changes/major/20170711VishalBoddu [moved from doc/news/changes/minor/20170711VishalBoddu with 77% similarity]
include/deal.II/lac/solver_fire.h
tests/lac/solver.cc
tests/lac/solver_fire_01.cc
tests/lac/solver_fire_rosenbrock_01.cc

similarity index 77%
rename from doc/news/changes/minor/20170711VishalBoddu
rename to doc/news/changes/major/20170711VishalBoddu
index c31f544988a66557c4dc318f08bd09d05ed325cc..1b4403326696146becbeb6ab9526c863fd5928c6 100644 (file)
@@ -1,4 +1,4 @@
 New: SolverFIRE implements FIRE (Fast Inertial Relaxation Engine) for solving
 the problem of minimization of a given objective function.
 <br>
-(Vishal Boddu, 2017/07/11)
+(Vishal Boddu, Denis Davydov 2017/07/11)
index 76e5b021de173df9b2b0c76591be3a98ab9a26d5..f6b7b0d6e97570e35f391602f8978c923c9a6c50 100644 (file)
@@ -1,3 +1,18 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2017 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.
+//
+// ---------------------------------------------------------------------
+
 #ifndef dealii__solver_fire_h
 #define dealii__solver_fire_h
 
@@ -18,8 +33,8 @@ DEAL_II_NAMESPACE_OPEN
 
 /**
  * FIRE (Fast Inertial Relaxation Engine) for minimization of (potentially
- * non-linear) objective function $E(U)$, $U$ is a vector of $n$ variables
- * ($n$ is the number of variables of the objective function).
+ * non-linear) objective function $E(\mathbf x)$, $\mathbf x$ is a vector of
+ * $n$ variables ($n$ is the number of variables of the objective function).
  * Like all other solver classes, it can work on any kind of vector and matrix
  * as long as they satisfy certain requirements (for the requirements on
  * matrices and vectors in order to work with this class, see the documentation
@@ -34,37 +49,36 @@ DEAL_II_NAMESPACE_OPEN
  * atomistic system, the algorithm relies on inertia to obtain (nearest)
  * configuration with least potential energy.
  *
- * The problem of solving the system of linear equations $AU = F$ for $U$
- * can be re-casted into the problem of finding $U$ for which the quadratic
- * function $\frac{1}{2} U^{T}AU - U^{T}F$ is minimized.
- *
  * Notation:
- * The global vector of unknowns:   $U$.        <BR>
- * Objective function:              $E(U)$.     <BR>
- * Rate of change of unknowns:      $V$.        <BR>
- * Gradient of the objective
- * function w.r.t unknowns:         $G$.        <BR>
- * Mass matrix:                     $M$.        <BR>
- * Initial guess of unknowns:       $U_0$.      <BR>
- * Time step:                       $\Delta t$. <BR>
+ *  - The global vector of unknown variables: $\mathbf x$.
+ *  - Objective function:                     $E(\mathbf x)$.
+ *  - Rate of change of unknowns:             $\mathbf v$.
+ *  - Gradient of the objective
+ *    function w.r.t unknowns:                $\mathbf g = \nabla E(\mathbf x)$.
+ *  - Mass matrix:                            $\mathbf M$.
+ *  - Initial guess of unknowns:              $\mathbf x_0$.
+ *  - Time step:                              $\Delta t$.
  *
- * Given initial values for $\Delta t$, $\alpha = \alpha_0$, $U = U_0$ and $V=0$
- * along with a given mass matrix $M$,
- * FIRE algorithm is as follows,
- * 1. Calculate $G = \nabla E(U)$ and check for convergence.
- *    Calculate $U$ and $V$ using any common Molecular Dynamics time integrator.
- *    Here we use simple Euler integration step, <BR>
- *        $U = U + \Delta t V$,                  <BR>
- *        $V = V + \Delta t M^{-1} \cdot G$.
- * 2. Calculate $P = G \cdot V$.
- * 3. Set $V = (1-\alpha) V + \alpha \frac{|V|}{|G|} G$.
- * 4. If $P \leq 0$ and number of steps since P was non-negative is larger than
- *    certain value, then increase time step $\Delta t$ and decrease $\alpha$.
- * 5. If $P>0$, then decrease the time step, freeze the system i.e., V = 0 and
- *    reset $\alpha = \alpha_0$.
- * 6. Return to 1.
+ * Given initial values for $\Delta t$, $\alpha = \alpha_0$, $\epsilon$,
+ * $\mathbf x = \mathbf x_0$ and $\mathbf v= \mathbf 0$ along with a given mass
+ * matrix $\mathbf M$, FIRE algorithm is as follows,
+ * 1. Calculate $\mathbf g = \nabla E(\mathbf x)$ and check for convergence
+ *    ($\mathbf g \cdot \mathbf g < \epsilon^2 $).
+ * 2. Update $\mathbf x$ and $V$ using simple (forward) Euler integration step,
+ *    <BR>
+ *        $\mathbf x = \mathbf x + \Delta t \mathbf v$,                 <BR>
+ *        $\mathbf v = \mathbf v + \Delta t \mathbf M^{-1} \cdot \mathbf g$.
+ * 3. Calculate $p = \mathbf g \cdot \mathbf v$.
+ * 4. Set $\mathbf v = (1-\alpha) \mathbf v
+ *                   + \alpha \frac{|\mathbf v|}{|\mathbf g|} \mathbf g$.
+ * 5. If $p<0$ and number of steps since $p$ was last negative is larger
+ *    than certain value, then increase time step $\Delta t$ and decrease
+ *    $\alpha$.
+ * 6. If $p>0$, then decrease the time step, freeze the system i.e.,
+ *    $\mathbf v = \mathbf 0$ and reset $\alpha = \alpha_0$.
+ * 7. Return to 1.
  *
- * Alse see
+ * Also see
  * <a href="http://onlinelibrary.wiley.com/doi/10.1002/pamm.201110246/full">
  * Energy-Minimization in Atomic-to-Continuum Scale-Bridging Methods </a> by
  * Eidel et al. 2011.
@@ -77,25 +91,33 @@ class SolverFIRE : public Solver<VectorType>
 
 public:
 
+  /**
+   * Standardized data struct to pipe additional data to the solver.
+   */
   struct AdditionalData
   {
+    /**
+     * Constructor. By default, set the initial time step for the (forward)
+     * Euler integration step to 0.1, the maximum time step to 1 and the
+     * maximum change allowed in any variable (per iteration) to 1.
+     */
     explicit
     AdditionalData (const double  initial_timestep    = 0.1,
                     const double  maximum_timestep    = 1,
                     const double  maximum_linfty_norm = 1);
 
     /**
-     * Initial time step.
+     * Initial time step for the (forward) Euler integration step.
      */
     const double initial_timestep;
 
     /**
-     * Maximum time step.
+     * Maximum time step for the (forward) Euler integration step.
      */
     const double maximum_timestep;
 
     /**
-     * Maximum change allowed in any degree of freedom.
+     * Maximum change allowed in any variable of the objective function.
      */
     const double maximum_linfty_norm;
 
@@ -106,7 +128,7 @@ public:
    */
   SolverFIRE (SolverControl            &solver_control,
               VectorMemory<VectorType> &vector_memory,
-              const AdditionalData     &data          );
+              const AdditionalData     &data = AdditionalData());
 
   /**
    * Constructor. Use an object of type GrowingVectorMemory as a default to
@@ -121,22 +143,26 @@ public:
   virtual ~SolverFIRE();
 
   /**
-   * Obtain a set of #p u (variables) that minimize an objective function
+   * Obtain a set of variables @p x that minimize an objective function
    * described by the polymorphic function wrapper @p compute, with a given
-   * preconditioner @p inverse_masses and initial @p u values.
-   * The function @p compute takes in the state of the (u) variables as
-   * argument and returns a pair of objective function's value and
-   * objective function's gradient (with respect to the variables).
+   * preconditioner @p inverse_mass_matrix and initial @p x values.
+   * The function @p compute returns the objective function's value and updates
+   * the objective function's gradient (with respect to the variables) when
+   * passed in as first argument based on the second argument-- the state of
+   * variables.
+   *
+   * @author Vishal Boddu, Denis Davydov, 2017
    */
   template<typename PreconditionerType = DiagonalMatrix<VectorType>>
   void solve
-  (std::function<double(VectorType &, const VectorType &)>  compute,
-   VectorType                                              &u,
-   const PreconditionerType                                &inverse_masses);
+  (const std::function<double(VectorType &, const VectorType &)> &compute,
+   VectorType                                                    &x,
+   const PreconditionerType                                      &inverse_mass_matrix);
 
   /**
-   * Solve for x that minimizes the quadratic function
-   * $E(x) = \frac{1}{2} x^{T}Ax - x^{T}b$.
+   * Solve for x that minimizes $E(\mathbf x)$ for the <EM>special case</EM>
+   * when $E(\mathbf x)
+   * = \frac{1}{2} \mathbf x^{T} \mathbf A \mathbf x - \mathbf x^{T} \mathbf b$.
    */
   template<typename MatrixType, typename PreconditionerType>
   void solve (const MatrixType         &A,
@@ -148,16 +174,17 @@ protected:
 
   /**
    * Interface for derived class. This function gets the current iteration
-   * u, u's time derivative and the gradient in each step. It can be used
-   * for a graphical output of the convergence history.
+   * @p x (variables), @p v (x's time derivative) and @p g (the gradient) in
+   * each step.
+   * It can be used for a graphical output of the convergence history.
    */
-  void print_vectors (const unsigned int,
-                      const VectorType &,
-                      const VectorType &,
-                      const VectorType &) const;
+  virtual void print_vectors (const unsigned int,
+                              const VectorType &x,
+                              const VectorType &v,
+                              const VectorType &g) const;
 
   /**
-   * Additional parameters.
+   * Additional data to the solver.
    */
   const AdditionalData additional_data;
 
@@ -222,9 +249,9 @@ template<typename VectorType>
 template<typename PreconditionerType>
 void
 SolverFIRE<VectorType>::solve
-(std::function<double(VectorType &, const VectorType &)>  compute,
- VectorType                                              &u,
- const PreconditionerType                                &inverse_masses)
+(const std::function<double(VectorType &, const VectorType &)> &compute,
+ VectorType                                                    &x,
+ const PreconditionerType                                      &inverse_mass_matrix)
 {
   deallog.push("FIRE");
 
@@ -242,20 +269,20 @@ SolverFIRE<VectorType>::solve
 
   // Set velocities to zero but not gradients
   // as we are going to compute them soon.
-  v->reinit(u,false);
-  g->reinit(u,true);
+  v->reinit(x,false);
+  g->reinit(x,true);
 
   // Refer to v and g with some readable names.
   VectorType &velocities = *v;
   VectorType &gradients  = *g;
 
-  // Update gradients for the new u.
-  compute(gradients, u);
+  // Update gradients for the new x.
+  compute(gradients, x);
 
   unsigned int iter = 0;
 
   SolverControl::State conv = SolverControl::iterate;
-  conv = this->iteration_status (iter, gradients * gradients, u);
+  conv = this->iteration_status (iter, gradients * gradients, x);
   if (conv != SolverControl::iterate)
     return;
 
@@ -272,15 +299,15 @@ SolverFIRE<VectorType>::solve
     {
       ++iter;
       // Euler integration step.
-      u.add (timestep, velocities);                // U += dt     * V
-      inverse_masses.vmult(gradients, gradients);  // G  = M^{-1} * G
-      velocities.add (-timestep, gradients);       // V -= dt     * G
+      x.add (timestep, velocities);                    // x += dt     * v
+      inverse_mass_matrix.vmult(gradients, gradients); // g  = M^{-1} * g
+      velocities.add (-timestep, gradients);           // v -= dt     * h
 
-      // Compute gradients for the new u.
-      compute(gradients, u);
+      // Compute gradients for the new x.
+      compute(gradients, x);
 
       const real_type gradient_norm_squared = gradients * gradients;
-      conv = this->iteration_status(iter, gradient_norm_squared, u);
+      conv = this->iteration_status(iter, gradient_norm_squared, x);
       if (conv != SolverControl::iterate)
         break;
 
@@ -332,7 +359,7 @@ SolverFIRE<VectorType>::solve
             timestep = minimal_timestep;
         }
 
-      print_vectors(iter, u, velocities, gradients);
+      print_vectors(iter, x, velocities, gradients);
 
     } // While we need to iterate.
 
@@ -356,16 +383,16 @@ void SolverFIRE<VectorType>::solve (const MatrixType         &A,
 {
 
   std::function<double(VectorType &,  const VectorType &)> compute_func =
-    [&]               (decltype(x) &G, decltype(b)       &x) -> double
+    [&]               (VectorType &g, const VectorType &x) -> double
   {
     // Residual of the quadratic form $ \frac{1}{2} xAx - xb $.
     // G = b - Ax
-    A.residual(G, x, b);
+    A.residual(g, x, b);
 
     // Gradient G = Ax -b.
-    G *= -1.;
+    g *= -1.;
 
-    // The quadratic form $ \frac{1}{2} x^t A x - x^{t} b $.
+    // The quadratic form $\frac{1}{2} xAx - xb $.
     return 0.5*A.matrix_norm_square(x) - x*b;
   };
 
@@ -376,11 +403,10 @@ void SolverFIRE<VectorType>::solve (const MatrixType         &A,
 
 template <typename VectorType>
 void
-SolverFIRE<VectorType>::
-print_vectors (const unsigned int,
-               const VectorType &,
-               const VectorType &,
-               const VectorType &) const
+SolverFIRE<VectorType>::print_vectors (const unsigned int,
+                                       const VectorType &,
+                                       const VectorType &,
+                                       const VectorType &) const
 {}
 
 
index ef88d74f6537c5e89621bfd83a34bfd3b81a8758..ae42de9b9ca904422d3bebef2c526a14cecc05a0 100644 (file)
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 1998 - 2016 by the deal.II authors
+// Copyright (C) 1998 - 2017 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -94,8 +94,7 @@ int main()
   SolverBicgstab<> bicgstab(control, mem);
   SolverRichardson<> rich(control, mem);
   SolverQMRS<> qmrs(control, mem);
-  SolverFIRE<>::AdditionalData fire_data(0.1, 1, 1);
-  SolverFIRE<> fire(control, mem, fire_data);
+  SolverFIRE<> fire(control, mem);
 
   for (unsigned int size=4; size <= 30; size *= 3)
     {
index f2f126882634cadfc8e5e7cc8986c9fa5ebc8c3f..24991c75a948dbd5cd5bf511e075c750742fcc12 100644 (file)
@@ -1,3 +1,23 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2017 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test to verify correctness of SolverFIRE::solve()
+// The objective function is f(x,y) = x^2 + y^2.
+
 
 #include "../tests.h"
 #include <fstream>
@@ -8,12 +28,6 @@
 #include <deal.II/lac/solver_control.h>
 #include <deal.II/lac/solver_fire.h>
 
-using namespace dealii;
-
-
-// Test to verify correctness of SolverFIRE::solve()
-// The objective function is f(x,y) = x^2 + y^2.
-
 
 using vector_t = typename dealii::Vector<double>;
 
index 51c8b3fefc856a0f36ff2b8173f19207d6f64f19..1e189f832ccddf23e6fe5f1cd155b6a2d7d08c71 100644 (file)
@@ -1,3 +1,17 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2017 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.
+//
+// ---------------------------------------------------------------------
 
 #include "../tests.h"
 #include <fstream>
@@ -9,7 +23,6 @@
 #include <deal.II/lac/solver_control.h>
 #include <deal.II/lac/solver_fire.h>
 
-using namespace dealii;
 
 
 // Test to verify correctness of SolverFIRE::sovle()

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.