Lots of changes to the documentation to reflect the recent changes. More

to be done obviously.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@140 2a26948c-0e4f-0410-aee8-f1d3e353619c

1 files changed, 168 insertions, 137 deletions

diff --git a/doc/documentation.tex b/doc/documentation.tex
index 95c214a..ba33704 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -84,7 +84,7 @@
 \title{Thorn Guide for the {\tt EHFinder} Thorn}
 
 % the date your document was last changed, if your document is in CVS, 
-% please us:
+% please use:
 %    \date{$ $Date$ $}
 \date{$ $Date$ $}
 
@@ -132,15 +132,18 @@ numerical spacetime by evolving a null surface backwards in time. This
 method depends on the fact that, except in cases where the coordinates are
 adapted to outgoing null geodesics, an outgoing null surface started close
 to the EH, when evolved forward in time, is diverging exponentially from the
-EH. Reversing the time evolutions then means that an outgoing null surface
+EH. Reversing the time evolution then means that an outgoing null surface
 will converge exponentially to the EH. The level set function, $f$, is
 evolved according to
-\begin{equation}
-\partial_{t}f = \frac{-g^{ti}\partial_{i}f+\sqrt{(g^{ti}\partial_{i}f)^{2}-
-g^{tt}g^{ij}\partial_{i}f\partial_{j}f}}{g^{tt}}.
+\begin{eqnarray}
+\partial_{t}f & = & \frac{-g^{ti}\partial_{i}f+\sqrt{(g^{ti}\partial_{i}f)^{2}-
+g^{tt}g^{ij}\partial_{i}f\partial_{j}f}}{g^{tt}} \nonumber \\
+ & = & \beta^{i}\partial_{i}f-
+                \sqrt{\alpha^{2}\gamma^{ij}\partial_{i}f\partial_{j}f},
 \label{AEIDevelopment_EHFinder_evolve}
-\end{equation}
-For more details on the theory and implementation
+\end{eqnarray}
+where in the second equation the lapse, shift and 3-metric has been substituted
+for the 4-metric. For more details on the theory and implementation
 see~\cite{AEIDevelopment_EHFinder_Diener02}.
 
 This thorn uses a level set description of the null surface, \ie the surface
@@ -159,15 +162,18 @@ consideration. Since the null surface has to be evolved backwards in time, the
 is necessary to evolve the initial data forward in time while outputting
 enough 3D data, that the full 4-metric can be recovered at each timestep.
 The 3D data is then read back in, in reverse order, while the level-set
-function is evolved backwards in time. More details about the actual use
-of the thorn in section~\ref{AEIDevelopment_EHFinder_UseThorn}
+function is evolved backwards in time. The thorn can evolve more than one
+level set function at a time using different initial guesses for the surfaces
+(NOTE: this is a quite recent feature and has not yet been tested extensively).
+More details about the actual use of the thorn in
+ section~\ref{AEIDevelopment_EHFinder_UseThorn}
 
-\section{Re-parametrization}
+\section{Re-initialization}
 \label{AEIDevelopment_EHFinder_reparam}
 The evolution equation for $f$, 
 equation~(\ref{AEIDevelopment_EHFinder_evolve}), causes steepening of
 the gradient of $f$, which is undesireble numerically. For that reason, $f$
-is periodically re-parametrized to a distance function. That is, the values of
+is periodically re-initialized to a distance function. That is, the values of
 $f$ are changed so that the the value of $f$ in a grid point is equal to the
 (signed) distance from the grid point to the surface $f=0$. This is done by
 evolving $f$ according to the following evolution equation (in the parameter
@@ -177,19 +183,15 @@ $\lambda$)
 \label{AEIDevelopment_EHFinder_reinit}
 \end{equation}
 until a steady state is achieved. This method is called the {\tt pde}-method
-since it is basically evolving a pde. Sometimes the $f=0$ can be moved
-slightly during the re-parametrization procedure. This happens when
+since it is basically evolving a pde. Sometimes the $f=0$ surface can be
+moved slightly during the re-initialization procedure. This happens when
 the surface develops a narrow neck just before a topology change. For
 that reason, there is an option to detect when this is about to happen and
-undo the re-parametrization.
-
-There is also another approximate way of performing the re-parametrization,
-however, even though it is faster, it is not recommended for use, since it
-seems to cause some instabilities and may cause larger movements of the $f=0$
-surface. This method is called {\tt approx}.
+undo the re-initialization.
 
-The two methods can be mixed, but it is not clear how well this works before
-some more experiments have been done.
+There used to be another method doing the re-initialization, but it proved
+inferior to the {\tt pde}-method and was removed. Other methods may be
+implemented in the future.
 
 \section{The initial shape of the surface}
 \label{AEIDevelopment_EHFinder_initial}
@@ -220,16 +222,18 @@ f = (x^{2}+y^{2}+z^{2})^{2} + a^{4} - 2 a^{2} (x^{2}-y^{2}-z^{2})-b^{4}.
 \end{equation}
 There are no translation or rotations implemented for the ovaloid of Cassini.
 
+Different initial shapes can be used for the different level set functions
+used in the same run.
 \section{Excision}
 \label{AEIDevelopment_EHFinder_excise}
 Even though the level set function, $f$, in principle can be defined
 everywhere it is often advantageous to only evolve it in a certain region
-around the surface $f=0$. Since $f$ is re-parametrized regularly to become
+around the surface $f=0$. Since $f$ is re-initialized regularly to become
 a distance function, $f$ itself can be used as a measure of the distance
 from a certain grid point to the surface $f=0$. The parameter
 {\tt ehfinder::shell\_width} specifies the size of the active region
-around $f=0$. However the interior and exterior excision is handled
-differently. The interior excision is most simply, since here all grid points
+around $f=0$. However the interior and exterior excisions are handled
+differently. The interior excision is most simple, since here all grid points
 with $f<-{\mbox{\tt ehfinder::shell\_width}}$ are marked as inactive. This
 should work in all cases when the excised region is everywhere concave, since
 then all points on the boundary of the excised region will have enough
@@ -250,24 +254,24 @@ Figure~\ref{AEIDevelopment_EHFinder_excisefig}, for the case
 \end{figure}
 {\tt ehfinder::shell\_width = 4}, where the excised regions are hashed.
 
-Changes to the excision regions are only done after re-parametrization,
+Changes to the excision regions are only done after re-initialization,
 since it is only at this time that $f$ is a distance function. The excision
-regions can move across the grid, following the surface $f=0$.
+regions can move across the grid, tracking the surface $f=0$.
 
-At present there is no support for using the {\tt SpaceMask} excision
-mask, but this should be straightforward to implement and will be done
-soon.
+Also if the numerical was run done with excision using {\tt SpaceMask},
+the {\tt EHFinder} excision region is guaranteed to completely cover the 
+numerical excision region.
 \section{Upwinding}
 \label{AEIDevelopment_EHFinder_upwind}
 All finite differences used in the evolution of the null surface are second
 order one sided differences. For that reason a {\tt ghost\_size} larger or
-equal to 2 should always be used. It is possible to choose between different
-upwinding schemes depending on whether the shift is non zero or not. This is
-done by setting the parameter {\tt ehfinder::upwind\_type} to either
-{\tt intrinsic} (for no shift) or {\tt shift} (for non zero shift).
+equal to 2 should always be used. It is possible to choose between three
+different upwinding schemes. This is done by setting the parameter
+{\tt ehfinder::upwind\_type} to either {\tt intrinsic}, {\tt shift} or
+{\tt characteristic} {\em [Note: get rid of the shift upwinding]}.
 
 The {\tt intrinsic} scheme, looks at the values of $f$ itself, to determine
-the direction of the stencil. This is basically to be able to handle
+the direction of the stencil. This is basically in order to be able to handle
 situations like the one illustrated in 1D in
 Figure~\ref{AEIDevelopment_EHFinder_upwindfig}.
 \begin{figure}[ht]
@@ -279,44 +283,54 @@ Figure~\ref{AEIDevelopment_EHFinder_upwindfig}.
   \label{AEIDevelopment_EHFinder_upwindfig}
 \end{figure}
 If the stencil for calculating derivatives in the point labeled 1 is taken
-to consist of the points 1, 2 and 3', the non differentiablility of $f$ will
-cause problems. The algoithm detects this and instead uses the points 1, 2, 3
+to consist of the points 1, 2 and 3', the non differentiablility of $f$ may
+cause problems. The algorithm detects this and instead uses the points 1, 2, 3
 as the stencil. This ensures that a non differentiable feature can be
 maintained in the evolution.
 
-When there is a shift present, it can happen that the direction chosen
-by the {\tt intrinsic} scheme is inconsistent with the shift direction,
-causing instabilities. Since re-parametrization is done, these
-instabilities are normally not allowed to develop fully, but they can
-still cause distortions of the surface $f=0$ and are therefore still
-damaging. For this reason it is possible to choose the upwinding direction
-based solely on the shift (since $f$ is evolved backwards in time, it is
-necessary to perform the upwinding in the opposite direction of the shift).
-
-It might happen that the upwinding direction based on the shift results
-the stencil to consist of the points 1, 2, 3' in 
-Figure~\ref{AEIDevelopment_EHFinder_upwindfig}. This might be fixed
-by increasing the frequency of re-parametrizations or by excising a
-larger region inside the surface. However if these workarounds doesn't
-work, it might be necessary to implement a hybrid upwinding scheme using
-the shift, except at points with this problem.
-
-For the re-parametrization the default is to use the {\tt intrinsic} 
-second order scheme (the re-parametrization doesn't depend on the shift,
-so {\tt shift} upwinding is not applicable. It is possible to use a
-first order {\tt intrinsic} scheme, but this is, in my experience, not 
-accurate enough. A centered differencing scheme is also available, but
-is only there for testing purposes and should never be used. These 
-alternative schemes will probably be removed in the future.
+However this scheme only works in a few simple cases. For more general cases
+it is important to use characteristic information in order to ensure that
+the numerical stencil contains the domain of dependence. Therefore the
+{\tt characteristic} scheme was introduced. Fortunately, by using the
+information contained in the level set functions the characteristic of
+the level set function can be calculated using
+\begin{equation}
+\frac{dx^{i}}{dt}=-\beta^{i}+\frac{\alpha^{2}\gamma^{ij}\partial_{j}f}
+{\sqrt{\alpha^{2}\gamma^{kl}\partial_{k}f\partial_{l}f}}. 
+\label{eq:char}
+\end{equation}
+The method estimates the characteristic direction using centered finite
+differences in equation~(\ref{eq:char}) and then recalculates the one
+sided finite differences in the appropriate direction.
+
+It might happen that the upwinding direction based on the characteristic
+direction results in the stencil to consist of the points 1, 2, 3' in 
+Figure~\ref{AEIDevelopment_EHFinder_upwindfig}. However, if the 
+re-initialization is done often enough, this turns out not to cause
+any problems.
+
+For the re-initialization the default is to use the {\tt intrinsic} 
+second order scheme (the re-initialization doesn't depend on where the
+surface is moving), so {\tt characteristic} upwinding is not applicable.
+It is possible to use a first order {\tt intrinsic} scheme, but this is, in my
+experience, not accurate enough. A centered differencing scheme is also
+available, but is only there for testing purposes and should never be used.
+These alternative schemes will probably be removed in the future.
 
 \section{The most important parameters}
 Here the most important parameters are described.
 \begin{itemize}
 \item {\tt ehfinder::mode} \\  
-  The mode can either be set to {\tt normal} (normal event horizon finder mode)
-  or {\tt test\_reparam} (mode to test the re-parametrization routine). The
-  default is {\tt normal} and should normally not be changed. This parameter
-  may be removed in the future.
+  The mode can either be set to {\tt normal} (normal event horizon finder
+  mode), {\tt test\_reparam} (mode to test the re-initialization routine),
+  {\tt analysis} (compare a previously calculated level set function to
+  a small number of analytic spacetimes) and {\tt generator} (only evolve
+  the generators while keeping the level set function fixed). The
+  default is {\tt normal} and should normally not be changed. The other modes
+  are only useful for debugging and testing purposes.
+\item{\tt ehfinder::eh\_number\_level\_sets} \\
+  An integer parameter specifying how many individual level set functions to
+  evolve at a time.
 \item {\tt ehfinder::eh\_metric\_type} \\
   The metric type can either be set to {\tt numerical} or {\tt analytic}. If
   it is set to {\tt numerical} the {\tt EHFinder} will attempt to read in the
@@ -330,45 +344,49 @@ Here the most important parameters are described.
   metric. It is possible to only set the metric on the initial slice, but it
   is also possible to have a thorn (like thorn {\tt Exact}) set the metric at
   {\tt CCTK\_PRESTEP} if the analytic metric is time dependent.
-  The default is at present {\tt analytic}. This should probably be changed.
+  The default is {\tt numerical}.
 \item {\tt ehfinder::eh\_lapse\_type} \\
   The same for the lapse.
 \item {\tt ehfinder::eh\_shift\_type} \\
   The same for the shift.
-\item {\tt initial\_f} \\
-  The initial shape of the null surface can currently be chosen from
-  {\tt sphere}, {\tt ellipsoid} and {\tt cassini} as described in
+\item {\tt initial\_f[i]} \\
+  A vector parameter specifying the initial shape of the null surface for
+  the individual level set functions. The initial shape can currently be
+  chosen from {\tt sphere}, {\tt ellipsoid} and {\tt cassini} as described in
   section~\ref{AEIDevelopment_EHFinder_initial}. The default is {\tt sphere}.
-\item {\tt initial\_rad} \\
-  The radius of the initial sphere ($r_{0}$ in
+\item {\tt initial\_rad[i]} \\
+  A vector parameter specifying the radius of the initial sphere ($r_{0}$ in
   equation~\ref{AEIDevelopment_EHFinder_sphere}). The deafault is 1.
-\item {\tt translate\_x} \\
-  How much to translate the initial surface in the $x$-direction ($x_{0}$ in
-  equation~\ref{AEIDevelopment_EHFinder_sphere}). Also used for the initial
-  ellipsoid. The default is 0.
-\item {\tt translate\_y} \\
-  How much to translate the initial surface in the $y$-direction ($y_{0}$ in
-  equation~\ref{AEIDevelopment_EHFinder_sphere}). Also used for the initial
-  ellipsoid. The default is 0.
+\item {\tt translate\_x[i]} \\
+  A vector parameter specifying how much to translate the initial surface in
+  the $x$-direction ($x_{0}$ in equation~\ref{AEIDevelopment_EHFinder_sphere}).
+  Also used for the initial ellipsoid. The default is 0.
+\item {\tt translate\_y[i]} \\
+  A vector parameter specifying how much to translate the initial surface in
+  the $y$-direction ($y_{0}$ in equation~\ref{AEIDevelopment_EHFinder_sphere}).
+  Also used for the initial ellipsoid. The default is 0.
 \item {\tt translate\_z} \\
-  How much to translate the initial surface in the $z$-direction ($z_{0}$ in
-  equation~\ref{AEIDevelopment_EHFinder_sphere}). Also used for the initial
-  ellipsoid. The default is 0.
-\item {\tt initial\_a} \\
-  $a$ in equation~\ref{AEIDevelopment_EHFinder_ellipsoid}. The default is 1.
-\item {\tt initial\_b} \\
-  $b$ in equation~\ref{AEIDevelopment_EHFinder_ellipsoid}. The default is 1.
-\item {\tt initial\_c} \\
-  $c$ in equation~\ref{AEIDevelopment_EHFinder_ellipsoid}. The default is 1.
-\item {\tt rotation\_alpha} \\
-  Rotation angle $\alpha$ for the ellipsoid around the $z$-axis. The default
-  is 0.
-\item {\tt rotation\_beta} \\
-  Rotation angle $\beta$ for the ellipsoid around the $y$-axis. The default
-  is 0.
-\item {\tt rotation\_gamma} \\
-  Rotation angle $\gamma$ for the ellipsoid around the $x$-axis. The default
-  is 0.
+  A vector parameter specifying how much to translate the initial surface in
+  the $z$-direction ($z_{0}$ in equation~\ref{AEIDevelopment_EHFinder_sphere}).
+  Also used for the initial ellipsoid. The default is 0.
+\item {\tt initial\_a[i]} \\
+  A vector parameter specifying $a$ in
+  equation~\ref{AEIDevelopment_EHFinder_ellipsoid}. The default is 1.
+\item {\tt initial\_b[i]} \\
+  A vector parameter specifying $b$ in
+  equation~\ref{AEIDevelopment_EHFinder_ellipsoid}. The default is 1.
+\item {\tt initial\_c[i]} \\
+  A vector parameter specifying $c$ in
+  equation~\ref{AEIDevelopment_EHFinder_ellipsoid}. The default is 1.
+\item {\tt rotation\_alpha[i]} \\
+  A vector parameter specifying the rotation angle $\alpha$ for the ellipsoid
+  around the $z$-axis. The default is 0.
+\item {\tt rotation\_beta[i]} \\
+  A vector parameter specifying the rotation angle $\beta$ for the ellipsoid
+  around the $y$-axis. The default is 0.
+\item {\tt rotation\_gamma[i]} \\
+  A vector parameter specifying the rotation angle $\gamma$ for the ellipsoid
+  around the $x$-axis. The default is 0.
 \item {\tt shell\_width} \\
   The width of the active evolution region. Grid points more than
   {\tt shell\_width} gridspacings away from the $f=0$ surface are marked
@@ -376,41 +394,42 @@ Here the most important parameters are described.
   section~\ref{AEIDevelopment_EHFinder_excise}. The default is $7$
   gridspacings.
 \item {\tt upwind\_type} \\
-  The type of upwinding to be used (either {\tt intrinsic} or {\tt shift}). See
-  the detailed description of the upwinding types in
-  section~\ref{AEIDevelopment_EHFinder_upwind}. The default is {\tt intrinsic}.
+  The type of upwinding to be used (either {\tt intrinsic} or 
+  {\tt characteristic}). See the detailed description of the upwinding types
+  in section~\ref{AEIDevelopment_EHFinder_upwind}. The default is 
+  {\tt characteristic}.
 \item {\tt reparam\_undo} \\
-  Should the re-parametrization be undone just before pinch-off or not as
+  Should the re-initialization be undone just before pinch-off or not as
   described in section~\ref{AEIDevelopment_EHFinder_reparam}. The
   default is {\tt "no"}.
 \item {\tt re\_param\_method} \\
-  Choose the re-parametrization method. The choices are {\tt approx},
+  Choose the re-initialization method. The choices are {\tt approx},
   {\tt pde} or {\tt mixed}. As described in
   section~\ref{AEIDevelopment_EHFinder_reparam} the recommended method is
   {\tt pde}. Only use {\tt approx} if you really know what you are doing. The
   default is {\tt pde}.
 \item {\tt re\_param\_int\_method} \\
-  Choose the integration method in the {\tt pde}-re-parametrization method.
+  Choose the integration method in the {\tt pde}-re-initialization method.
   Choose either a simple Euler ({\tt euler}) integration scheme or a second
   order Runge-Kutta ({\tt rk2}) scheme. Since a pde is evolved to steady state,
   it seems that the Euler scheme works just fine and is faster than the
   Runge-Kutta scheme. The default is {\tt euler}.
 \item {\tt re\_param\_max\_iter} \\
-  The maximum number of iterations in the {\tt pde}-re-parametrization scheme,
+  The maximum number of iterations in the {\tt pde}-re-initialization scheme,
   before giving up. Unless you are working at high resolution the default
   should be enough. The default is 800.
 \item {\tt pde\_differences} \\
   Choose the type of finite differencing used in the {\tt pde}
-  re-parametrization. Don't ever use anything else than second order
-  uwinding ({\tt upwind2}). The other choices ({\tt centered} and {\tt upwind})
+  re-initialization. Don't ever use anything else than second order
+  upwinding ({\tt upwind2}). The other choices ({\tt centered} and {\tt upwind})
   are there only for testing purposes and might be removed.
 \item {\tt reparametrize\_every\_pde} \\
-  How often to re-parametrize using the {\tt pde} re-parametrization. This
+  How often to re-initialize using the {\tt pde} re-initialization. This
   depends on the problem. For some problems it is necessary to do it more
   often than for other problems. You'll have to experiment to figure out
-  what works best. The default (too high) is 100.
+  what works best. The default is 10.
 \item {\tt reparametrize\_every\_approx} \\
-  How often to re-parametrize using the {\tt approx} re-parametrization.
+  How often to re-initialize using the {\tt approx} re-initialization.
   Again it depends on the problem. Can be used together with
   {\tt reparametrize\_every\_pde} if {\tt re\_param\_method} is set to
   {\tt mixed}. The default is 10.
@@ -441,8 +460,13 @@ done with zero shift and/or lapse equal to one, it is not necessary to output
 these grid functions, as long as storage is turned on and they are set to the
 correct value initially when {\tt EHFinder} is run. If the evolution was done
 with a conformal factor then {\tt StaticConformal::psi} has to be output as
-well, since it is necessary in order to reconstruct the 4-metric. It is not
-necessary to output the derivatives of the conformal factor.
+well, since it is necessary in order to reconstruct the 4-metric. However it
+is only necessary to output is once since it is constant during the evolution.
+It is not necessary to output the derivatives of the conformal factor. If
+excision was used in the numerical run then it is also necessary to output 
+{\tt SpaceMask::emask}\footnote{When more thorns have been converted to use
+{\tt SpaceMask::space\_mask} this option will also be supported.}, since
+{\tt EHFinder} has to make sure that its internal mask covers the space mask.
 
 At present it is necessary to output all timesteps into the same file (use
 {\tt IO::out\_timesteps\_per\_file = -1}, which is the default). In principle
@@ -451,13 +475,13 @@ HDF5 output has been tested. Since {\tt EHFinder} can be run on a lot less
 processors compared to the spacetime evolution, it is often advantageous to
 either do unchuncked output or to recombine the output files, since it is then
 possible to read the data onto a smaller number of processors (use
-{\tt IO::out\_unchunked = "yes"}). If the numerical
+{\tt IO::out\_unchunked = "yes"} to write unchunked data). If the numerical
 run is larger than the EH containing region (hopefully that is the case;
 otherwise the boundaries are definitely to close in), it is possible to
 use hyperslabbing to just output the EH containing region (see for example
 {\tt CactusPUGHIO/IOHDF5} for details on this). If hyperslabbing is used it
-is definitely necessary to do the output unchunked. An example parameter
-file can be seen in the {\tt par/Misner\_2.2\_80\_3D.par}.
+is definitely necessary to do the output unchunked or recombine afterwards.
+An example parameter file can be seen in the {\tt par/Misner\_2.2\_80\_3D.par}.
 
 In principle {\tt EHFinder} should also work for downsampled (in both space
 and time) data, but no experiments have been done so far to estimate the loss
@@ -475,38 +499,45 @@ be used as a black box. But still I can give some guidelines and advice on
 how to proceed.
 
 The first concern is to setup the initial guess for the surface. Ideally one
-would like to make at least two runs with {\tt EHFinder}. One run with an
-initial guess for the surface completely inside the EH and one run with an
-initial guess completely outside the EH. The easiest way to get an initial
-surface inside the EH, is to set up the initial guess to be completely inside
-the
+would like to use at least two different initial surfaces. One surface 
+completely inside the EH and one surface completely outside the EH. In
+practice it is often a good idea to use more than two different initial
+surfaces, since then the initial surfaces closest to the EH can be identified.
+Using the feature of evolving multiple level set functions at a time, this
+can be done while reading in the numerical data only once. The easiest way
+to get an initial surface inside the EH, is to set up the initial guess to
+be completely inside the
 apparent horizon (AH). To get an initial guess that is outside of the EH
 is not as easy. One way is to choose a surface, that starts to contract
 everywhere when evolved according to
-equation~(\ref{AEIDevelopment_EHFinder_evolve}. This of course means to do
-it by trial and error. Set up some initial guess evolve it for a little while,
-look at 3D output to determine if the surface is contracting everywhere and
-change the initial surface if necessary.
-
-Then comes the question of how often to do the re-parametrization and how
+equation~(\ref{AEIDevelopment_EHFinder_evolve}). However this is not a
+guarantee, since the EH can be expanding in the numerical coordinates.
+This of course means that it is necessary to do it by trial and error.
+Set up some initial guess evolve it for a little while, look at 3D output
+to determine if the surface is contracting everywhere and change the
+initial surface if necessary.
+
+Then comes the question of how often to do the re-initialization and how
 much to excise. These parameters depend on the numerical data. In principle,
-since the re-parametrization can move the surface, one wants to do it as
-rarely as possible. On the other hand, re-parametrization is necessary
+since the re-initialization can move the surface, one wants to do it as
+rarely as possible. On the other hand, re-initialization is necessary
 in order to keep the evolution nicely controlled (by avoiding large gradients),
 so a compromise has to be found. This might require some experimentation.
-Because movements of the surface during re-parametrization, usually only
+Because movements of the surface during re-initialization, usually only
 occurs close to moments of topology change, it might be necessary to evolve
 all the way beyond the change of topology and look at 3D output to see if
-any problems occured. How often to do the re-parametrization also depends
-of the width of the active region. If the active region around the surface
-is narrow, it might be necessary to re-parametrize more often, since in this
-case the boundaries of the active region is closer to the surface. At the
-boundaries the stencil direction is dictated by the geometry and not
-$f$ itself or the shift, which might cause instabilities if it is not
-re-parametrized. Good initial guesses for
-{\tt ehfinder::reparametrize\_every\_pde} seems to in the range 5--10, however
-I have seen cases where less were required and others where more were
-possible. For {\tt ehfinder::shell\_width} I normally use at least 7.
+any problems occured. If significant movement of the surface during
+re-initialization near the change of topology is observed, then try again
+with re-initialization undo activated. How often to do the re-initialization
+also depends on the width of the active region. If the active region around
+the surface is narrow, it might be necessary to re-initialize more often,
+since in this case the boundaries of the active region is closer to the
+surface. At the boundaries the stencil direction is dictated by the geometry
+and not $f$ itself or the characteristics, which might cause instabilities
+if it is not re-initialized. Good values guesses for
+{\tt ehfinder::reinitialize\_every\_pde} seems to in the range 5--10 at
+low or medium resolutions but can usually be increased for higher
+resolutions. For {\tt ehfinder::shell\_width} I normally use at least 7.
 
 This documentation will be updated, as input comes in from users.