From 20c0aa0ea5a7b85403523e2565e8a5e0ebee5c45 Mon Sep 17 00:00:00 2001 From: vishalkenchan Date: Wed, 12 Jul 2017 10:30:53 +0200 Subject: [PATCH] updated documentation --- .../{minor => major}/20170711VishalBoddu | 2 +- include/deal.II/lac/solver_fire.h | 176 ++++++++++-------- tests/lac/solver.cc | 5 +- tests/lac/solver_fire_01.cc | 26 ++- tests/lac/solver_fire_rosenbrock_01.cc | 15 +- 5 files changed, 138 insertions(+), 86 deletions(-) rename doc/news/changes/{minor => major}/20170711VishalBoddu (77%) diff --git a/doc/news/changes/minor/20170711VishalBoddu b/doc/news/changes/major/20170711VishalBoddu similarity index 77% rename from doc/news/changes/minor/20170711VishalBoddu rename to doc/news/changes/major/20170711VishalBoddu index c31f544988..1b44033266 100644 --- a/doc/news/changes/minor/20170711VishalBoddu +++ b/doc/news/changes/major/20170711VishalBoddu @@ -1,4 +1,4 @@ New: SolverFIRE implements FIRE (Fast Inertial Relaxation Engine) for solving the problem of minimization of a given objective function.
-(Vishal Boddu, 2017/07/11) +(Vishal Boddu, Denis Davydov 2017/07/11) diff --git a/include/deal.II/lac/solver_fire.h b/include/deal.II/lac/solver_fire.h index 76e5b021de..f6b7b0d6e9 100644 --- a/include/deal.II/lac/solver_fire.h +++ b/include/deal.II/lac/solver_fire.h @@ -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$.
- * Objective function: $E(U)$.
- * Rate of change of unknowns: $V$.
- * Gradient of the objective - * function w.r.t unknowns: $G$.
- * Mass matrix: $M$.
- * Initial guess of unknowns: $U_0$.
- * Time step: $\Delta t$.
+ * - 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,
- * $U = U + \Delta t V$,
- * $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, + *
+ * $\mathbf x = \mathbf x + \Delta t \mathbf v$,
+ * $\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 * * Energy-Minimization in Atomic-to-Continuum Scale-Bridging Methods by * Eidel et al. 2011. @@ -77,25 +91,33 @@ class SolverFIRE : public Solver 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 &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> void solve - (std::function compute, - VectorType &u, - const PreconditionerType &inverse_masses); + (const std::function &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 special case + * when $E(\mathbf x) + * = \frac{1}{2} \mathbf x^{T} \mathbf A \mathbf x - \mathbf x^{T} \mathbf b$. */ template 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 template void SolverFIRE::solve -(std::function compute, - VectorType &u, - const PreconditionerType &inverse_masses) +(const std::function &compute, + VectorType &x, + const PreconditionerType &inverse_mass_matrix) { deallog.push("FIRE"); @@ -242,20 +269,20 @@ SolverFIRE::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::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::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::solve (const MatrixType &A, { std::function 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::solve (const MatrixType &A, template void -SolverFIRE:: -print_vectors (const unsigned int, - const VectorType &, - const VectorType &, - const VectorType &) const +SolverFIRE::print_vectors (const unsigned int, + const VectorType &, + const VectorType &, + const VectorType &) const {} diff --git a/tests/lac/solver.cc b/tests/lac/solver.cc index ef88d74f65..ae42de9b9c 100644 --- a/tests/lac/solver.cc +++ b/tests/lac/solver.cc @@ -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) { diff --git a/tests/lac/solver_fire_01.cc b/tests/lac/solver_fire_01.cc index f2f1268826..24991c75a9 100644 --- a/tests/lac/solver_fire_01.cc +++ b/tests/lac/solver_fire_01.cc @@ -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 @@ -8,12 +28,6 @@ #include #include -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; diff --git a/tests/lac/solver_fire_rosenbrock_01.cc b/tests/lac/solver_fire_rosenbrock_01.cc index 51c8b3fefc..1e189f832c 100644 --- a/tests/lac/solver_fire_rosenbrock_01.cc +++ b/tests/lac/solver_fire_rosenbrock_01.cc @@ -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 @@ -9,7 +23,6 @@ #include #include -using namespace dealii; // Test to verify correctness of SolverFIRE::sovle() -- 2.39.5