File:NewtonHook.f90: Difference between revisions

From openpipeflow.org
Jump to navigation Jump to search
mNo edit summary
mNo edit summary
 
(50 intermediate revisions by the same user not shown)
Line 1: Line 1:
{{latexPreamble}}
{{latexPreamble}}
Click link just above to download.  The code '''requires''' [[File:GMRESm.f90]] and LAPACK.


== The method ==


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.
'''See also Template/Example for FORTRAN / MATLAB at [[Newton-Krylov_method]]'''


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}$.
'''*** THE DESCRIPTION OF THE METHOD HAS MOVED TO [[Newton-Krylov_method]] ***'''


For each iteration, $\vec{x}_{i+1} = \vec{x}_i + \vec{\delta x}$,
== The Jacobian-free Newton-Krylov (JFNK) Method ==
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]].
 
The code above implements the '''Jacobian-free Newton-Krylov (JFNK) method''' for solving
\begin{equation}\label{eq:ROOTF}\vec{F}(\vec{x})=\vec{0},\end{equation}  
where $\vec{x}$ and $\vec{F}$ are $n$-vectors, '''supplemented with a Hookstep--Trust-region approach'''.


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.


== The code ==
== The code ==


To download, click the link above.
To download, click the links at top of page.  A template for calling the routine: [[File:Newton_template.f90]].  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> needs to be assigned the initial guess $\vec{x}_0$ prior to calling <tt>newtonhook(...)</tt>.
* Also prior to calling <tt>newtonhook(...)</tt>, <tt>new_x(:)</tt> must be assigned the initial guess $\vec{x}_0$. (As the calculation proceeds, the vector <tt>new_x(:)</tt> is updated with improved vectors, having smaller corresponding <tt>new_fx(:)</tt>.)
 
* 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,
** a function that calculates $\vec{F}(\vec{x})$ for a given $\vec{\delta x}$,
** a function that calculates $\vec{F}(\vec{x})$ for a given $\vec{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 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, where $\vec{F}(\vec{x}_i)$ is already stored in <tt>new_fx(:)</tt>.
** 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 (smaller corresponding <tt>new_f(:)</tt>).
 
* The functions above may require auxiliary data, in addition to the given $\vec{x}$ or $\vec{\delta x}$.  A few ways to access the data:
** place the required data in a module and access via '<tt>use mymodule</tt>'.
** place the data in a <tt>common</tt> block.
** write the vector to disk, call a script that runs an external programme, then load the result.
 
* Extra constraints may be necessary, e.g. that determine an update to the period of an orbit.  The function that evaluates $\left.\frac{\vec{dF}}{\vec{dx}}\right|_{\vec{x}_i}\,\vec{\delta x}$ for a given $\vec{\delta x}$ should append to the result the evaluations of the constraints.  Correspondingly, for each constraint an extra <tt>0</tt> should be appended to $\vec{F}(\vec{x}_i)$ (the evaluated constraint should equal zero when converged). 
 
== Parallel use ==
 
It is NOT necessary to edit this code for parallel (MPI) use:
* let each thread pass its subsection of the vector $\vec{x}$,
* make the dot product  function <tt>mpi_allreduce</tt> the result of the dot product.
* to avoid multiple outputs to the terminal, set <tt>info=1</tt> on rank 0 and <tt>info=0</tt> for the other ranks.

Latest revision as of 04:12, 1 February 2023

$ \renewcommand{\vec}[1]{ {\bf #1} } \newcommand{\bnabla}{ \vec{\nabla} } \newcommand{\Rey}{Re} \def\vechat#1{ \hat{ \vec{#1} } } \def\mat#1{#1} $ Click link just above to download. The code requires File:GMRESm.f90 and LAPACK.


See also Template/Example for FORTRAN / MATLAB at Newton-Krylov_method

*** THE DESCRIPTION OF THE METHOD HAS MOVED TO Newton-Krylov_method ***

The Jacobian-free Newton-Krylov (JFNK) Method

The code above implements the Jacobian-free Newton-Krylov (JFNK) method for solving \begin{equation}\label{eq:ROOTF}\vec{F}(\vec{x})=\vec{0},\end{equation} where $\vec{x}$ and $\vec{F}$ are $n$-vectors, supplemented with a Hookstep--Trust-region approach.


The code

To download, click the links at top of page. A template for calling the routine: File:Newton template.f90. 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$. (As the calculation proceeds, the vector new_x(:) is updated with improved vectors, having smaller corresponding new_fx(:).)
  • 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{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, where $\vec{F}(\vec{x}_i)$ is already stored in new_fx(:).
    • 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, in addition to the given $\vec{x}$ or $\vec{\delta x}$. A few ways to access the data:
    • place the required data in a module and access via 'use mymodule'.
    • place the data in a common block.
    • write the vector to disk, call a script that runs an external programme, then load the result.
  • Extra constraints may be necessary, e.g. that determine an update to the period of an orbit. The function that evaluates $\left.\frac{\vec{dF}}{\vec{dx}}\right|_{\vec{x}_i}\,\vec{\delta x}$ for a given $\vec{\delta x}$ should append to the result the evaluations of the constraints. Correspondingly, for each constraint an extra 0 should be appended to $\vec{F}(\vec{x}_i)$ (the evaluated constraint should equal zero when converged).

Parallel use

It is NOT necessary to edit this code for parallel (MPI) use:

  • let each thread pass its subsection of the vector $\vec{x}$,
  • make the dot product function mpi_allreduce the result of the dot product.
  • to avoid multiple outputs to the terminal, set info=1 on rank 0 and info=0 for the other ranks.

File history

Click on a date/time to view the file as it appeared at that time.

Date/TimeDimensionsUserComment
current04: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.

Retrieved from ‘http:///index.php?title=File:NewtonHook.f90&oldid=1054’