File:NewtonHook.f90: Difference between revisions
mNo edit summary |
m (→The code) |
||
Line 18: | Line 18: | ||
== The code == | == The code == | ||
To download, click the link above. | To download, click the link above. The code also requires [[File:GMRESm.f90]] and Lapack. | ||
* Prior to calling the subroutine <tt>newtonhook(...)</tt>, the arrays <tt>new_x(:)</tt>$\equiv\vec{x}_i$ and <tt>new_fx(:)</tt>$\equiv\vec{F}(\vec{x}_i)$ need to be allocated the dimension $n$. | * Prior to calling the subroutine <tt>newtonhook(...)</tt>, the arrays <tt>new_x(:)</tt>$\equiv\vec{x}_i$ and <tt>new_fx(:)</tt>$\equiv\vec{F}(\vec{x}_i)$ need to be allocated the dimension $n$. | ||
* <tt>new_x(:)</tt> | * Also prior to calling <tt>newtonhook(...)</tt>, <tt>new_x(:)</tt> must be assigned the initial guess $\vec{x}_0$. | ||
* In addition to scalar and array variables, <tt>newtonhook(...)</tt> needs to be passed | * In addition to scalar and array variables, <tt>newtonhook(...)</tt> needs to be passed | ||
** a function that calculates the dot product of two vectors, | ** a function that calculates the dot product of two vectors, | ||
Line 28: | Line 28: | ||
** a subroutine that replaces a vector $\vec{x}$ with the solution of $M\vec{x}'=\vec{x}$ for $\vec{x}'$, where $M$ is a preconditioner matrix. This may simply be an empty subroutine if no preconditioner is required, i.e. $M=I$. | ** a subroutine that replaces a vector $\vec{x}$ with the solution of $M\vec{x}'=\vec{x}$ for $\vec{x}'$, where $M$ is a preconditioner matrix. This may simply be an empty subroutine if no preconditioner is required, i.e. $M=I$. | ||
** a subroutine that is called at the end of each iteration. This may be used to save the current state after each iteration, if desired. | ** a subroutine that is called at the end of each iteration. This may be used to save the current state after each iteration, if desired. | ||
* As the calculation proceeds, the vector <tt>new_x(:)</tt> is updated with improved vectors | * The functions above may require auxiliary data to $\vec{x}$. Place this data in a module and access via '<tt>use mymodule</tt>' in the function/subroutine. | ||
* As the calculation proceeds, the vector <tt>new_x(:)</tt> is updated with improved vectors, having smaller corresponding <tt>new_f(:)</tt>. | |||
== Parallel use == | == Parallel use == | ||
No special adaptations are required for parallel (MPI) use -- let each thread pass its subsection for the vector <tt>sv</tt>, and let the <tt>dotprod</tt> function <tt>allreduce</tt> the result of the dot product. | No special adaptations are required for parallel (MPI) use -- let each thread pass its subsection for the vector <tt>sv</tt>, and let the <tt>dotprod</tt> function <tt>allreduce</tt> the result of the dot product. |
Revision as of 02:45, 16 December 2016
$ \renewcommand{\vec}[1]{ {\bf #1} } \newcommand{\bnabla}{ \vec{\nabla} } \newcommand{\Rey}{Re} \def\vechat#1{ \hat{ \vec{#1} } } \def\mat#1{#1} $
The Newton-Krylov method
The code above implements the Newton-Raphson-Krylov method for solving \[\vec{F}(\vec{x})=\vec{0},\] where $\vec{x}$ and $\vec{F}$ are $n$-vectors, using a Hookstep-Trust-region approach.
This is a powerful method that can solve for $\vec{x}$ for complex $\vec{F}(\vec{x})$. E.g., for a periodic orbit, let $\vec{F}(\vec{x}) = \vec{X}(\vec{x})-\vec{x}$, where $\vec{X}(\vec{x})$ is the result of time integration of an initial condition $\vec{x}$.
For each iteration, $\vec{x}_{i+1} = \vec{x}_i + \vec{\delta x}$, the linear problem \[\left.\frac{\vec{dF}}{\vec{dx}}\right|_{\vec{x}_i}\,\vec{\delta x}=-\vec{F}(\vec{x}_i) \] needs to be solved for the step $\vec{\delta x}$ subject to the condition that it be smaller than some magnitude $\delta$, where $\delta$ is the size of the trust region. This problem is in the form $A\vec{x}=\vec{b}$, where $A$ is an $n\times n$ matrix, and is solved using the Krylov-subspace method, GMRES(m). See File:GMRESm.f90.
To perform the GMRES(m) operations, for a given $\vec{\delta x}$, the product $\left.\frac{\vec{dF}}{\vec{dx}}\right|_{\vec{x}_i}\,\vec{\delta x}$ may be approximated, by e.g. $\frac{1}{\epsilon}(\vec{F}(\vec{x}_i+\epsilon\,\vec{\delta x})-\vec{F}(\vec{x}_i))$ for some small value $\epsilon$. Thus, the method does need to know the matrix $\left.\frac{\vec{dF}}{\vec{dx}}\right|_{\vec{x}_i}$ itself, only a routine for calculating $\vec{F}(\vec{x})$ is required.
Note that provided that each step $\vec{\delta x}$ of the Newton method is essentially in the correct direction, the method is expected to converge. Thus the tolerance specified in the accuracy of the solution for $\vec{\delta x}$, calculated via the GMRES method, typically need not be so stringent as the tolerance placed on the Newton method for the solution $\vec{x}$.
For further details of the method, see e.g. Newton-Krylov-Hookstep (channelflow.org).
The code
To download, click the link above. The code also requires File:GMRESm.f90 and Lapack.
- Prior to calling the subroutine newtonhook(...), the arrays new_x(:)$\equiv\vec{x}_i$ and new_fx(:)$\equiv\vec{F}(\vec{x}_i)$ need to be allocated the dimension $n$.
- Also prior to calling newtonhook(...), new_x(:) must be assigned the initial guess $\vec{x}_0$.
- In addition to scalar and array variables, newtonhook(...) needs to be passed
- a function that calculates the dot product of two vectors,
- a function that calculates $\vec{F}(\vec{x})$ for a given $\vec{\delta x}$,
- a function that calculates $\left.\frac{\vec{dF}}{\vec{dx}}\right|_{\vec{x}_i}\,\vec{\delta x}$ for a given $\vec{\delta x}$. See possible approximation above.
- a subroutine that replaces a vector $\vec{x}$ with the solution of $M\vec{x}'=\vec{x}$ for $\vec{x}'$, where $M$ is a preconditioner matrix. This may simply be an empty subroutine if no preconditioner is required, i.e. $M=I$.
- a subroutine that is called at the end of each iteration. This may be used to save the current state after each iteration, if desired.
- The functions above may require auxiliary data to $\vec{x}$. Place this data in a module and access via 'use mymodule' in the function/subroutine.
- As the calculation proceeds, the vector new_x(:) is updated with improved vectors, having smaller corresponding new_f(:).
Parallel use
No special adaptations are required for parallel (MPI) use -- let each thread pass its subsection for the vector sv, and let the dotprod function allreduce the result of the dot product.
File history
Click on a date/time to view the file as it appeared at that time.
Date/Time | Dimensions | User | Comment | |
---|---|---|---|---|
current | 04:50, 26 June 2019 | (6 KB) | Apwillis (talk | contribs) | Minor edits in comments only. |
05:32, 13 December 2016 | (6 KB) | Apwillis (talk | contribs) | == The method == Implements the Newton-Raphson method for solving \[\vec{F}(\vec{x})=\vec{0}\] with a Hookstep-Trust-region approach. At each iteration, the linear problem for the step in $\vec{x}$ is solved subject to the condition that the step be ... |
You cannot overwrite this file.
File usage
The following 4 pages use this file: