Surface Evolver Documentation

Back to top of Surface Evolver documentation.       Index.

Gradient Descent Iteration

Basic surface evolution is done by gradient descent iteration. Relevant topics are:

Gradient descent iteration step

Each g command iteration does the following:

Scale factor

Once a direction of motion is found, the direction must be multiplied by a scale factor to compute the actual motion. If one interprets the direction of motion as a velocity (as in motion by mean curvature), then the scale factor becomes a time step. The scale factor may be fixed with the m command, or it may be in optimizing mode. The default is to start in optimizing mode.

Optimizing scale

In using gradient descent to seek a minimum energy, one finds a direction of motion and does a line search along that direction to find the minimum energy. Evolver will do that in optimizing scale mode. The line search consists of halving or doubling the current scale factor until an energy minimum is bracketed; then quadratic interpolation is used to estimate the optimum scale. Optimizing scale is the default; it also may be turned on with the m command or the optimizing command.

For safety, there is an upper bound to the scale; it defaults to 1 but may be changed with the optimizing command. There is also a lower bound; if Evolver gets a scale below 1e-12 of the scale bound when attempting to find a minimum, it gives up and just uses scale 0. Scale 0 is not a null operation since it still projects to constraints, if they are not exactly satisfied.

In general, a good scale factor depends on the type of energy being minimized and the level of refinement. However, for minimizing area, when the triangulation is well-behaved and area normalization is off, the best scale factor is usually around 0.2, independent of refinement. In optimizing mode, a scale factor getting small, say below 0.01, indicates triangulation problems. Too large a fixed scale factor will show up as total energy increasing. If you have motion by area normalization ON use a small scale factor, like 0.001, until you get a feel for what works.

If check_increase is toggled on, then the motion is not done if it would increase energy. But be aware that energy sometimes may have to increase in order to satisfy constraints.

Conjugate gradient

"Conjugate gradient" is a method of accelerating gradient descent. In ordinary gradient descent, one uses the gradient of energy to find the steepest downhill direction, then moves along that line to the minimum energy in that direction. Hence successive steps are at right angles. However, this can be very inefficient, as you can spend a lot of time zigzagging across an energy "valley" without making much progress "downstream". With conjugate gradient, the search direction is chosen to be in a "conjugate" direction to the previous direction. For a mathematical explanation, see any decent book in numerical analysis, such as [P]. In practice, the conjugate gradient method remembers a cumulative "history vector", which it combines with the ordinary gradient to figure out the conjugate gradient direction. The upshot is that conjugate gradient can converge much faster than ordinary gradient descent.

Conjugate gradient can be toggled with the U command, or with the conj_grad toggle. It should always be used with optimizing scale.

Notes: The conjugate gradient method is designed for quadratic energy functions. As long as the energy function is nearly quadratic, as it should be near an energy minimum, conjugate gradient works well. Otherwise, it may misbehave, either by taking too big steps or by getting stalled. Both effects are due to the history vector being misleading. To prevent too big steps, one should iterate without conjugate gradient for a few steps whenever significant changes are made to the surface (refining, changing a constraint, etc.). On the other hand, if it looks like conjugate gradient is converging, it may have simply become confused by its own history. See the catenoid example for a case in point. A danger signal is the scale factor going to zero. If you are suspicious, toggle conjugate gradient off and on ("U 2" does nicely) to erase the history vector and start over.

Diffusion

The Evolver can simulate the real-life phenomenon of gas diffusion between neighboring bubbles. This diffusion is driven by the pressure difference across a surface. This is invoked by the keyword DIFFUSION in the first part of the datafile, followed by the value of the diffusion constant. The amount diffused across a facet during an iteration is calculated as scale*diffusion_constant*facet_area*pressure_difference. The scale factor is included as the time step of an iteration. The amount is added to or subtracted from the prescribed volumes of the bodies on either side of the facet. The diffusion constant can be accessible at runtime through the read-write internal variable diffusion_coeff.

If you want finer control over the rate of diffusion across various surfaces, you can define the edge_diffusion edge attribute in the string model or the facet_diffusion facet attribute in the soapfilm model and give individual values for edges or facets as you desire. If the attribute is defined, then its value is used instead of the global diffusion constant.

Stability

The timestep of an iteration should not be so large as to amplify perturbations of the surface. Short wavelength perturbations are most prone to amplification. This section contains a sketch of the stability characteristics of the various mobility modes, enough to let the user relate the maximum timestep to the minimum facet or edge size. Two examples are discussed: a zigzag string and a nearly flat surface with equilateral triangulation. Effective area is not included, as it is an insignificant correction for nearly flat surfaces. The general moral of this section is that the maximum time step in iteration is limited by the length of the shortest edge or the area of the smallest facet, except in one case.

Zigzag string

Let the amplitude of the perturbation about the midline be Y and the edge length L. Then the force on a vertex is F = -4Y/L for small Y. Let the timestep (the Evolver scale factor) be \Delta t. Let V be the vertex velocity. Then the critical timestep for amplification of the perturbation is given by V\Delta t = -2Y, or \Delta t = -2Y/V.

Vertex mobility. Here V = F, so \Delta t = L/2.

Area normalization. Here the vertex star has length 2L, so V = F/L and \Delta t = L^2/2.

Approximate curvature. It turns out that the zigzag is an eigenvector of the mobility matrix M for the largest eigenvalue 3/L, so V = 3F/L and \Delta t = L^2/6. This is a major disadvantage of approximate curvature. If perturbation instability is the limitation on the timestep, it will take three times as many iterations as with area normalization to do the same evolution.

Perturbed sheet with equilateral triangulation.

Consider a plane surface triangulated with equilateral triangles of area A. The perturbation consists of a tiling of hexagonal dimples with their centers of height Y above their peripheries. The force at a central vertex turns out to be -sqrt(3)Y and the force at a peripheral vertex sqrt(3)Y/2.

Vertex mobility. The critical time step is given by
(\sqrt(3)Y + sqrt(3)Y/2)\Delta t = 2Y,
so \Delta t = 4/3*sqrt(3). Note that this is independent of the triangle size. This is consistent with experience in evolving with optimizing scale factor, where the optimum time step is in the range 0.2 - 0.3 indenpendent of triangle size. This is a definite advantage of this version of mobility, since different parts of the surface can have different size triangulations and one size time step can work for all.

Area normalization. The star area of each vertex is 6A, so the velocities become -sqrt(3)Y/2A and sqrt(3)Y/4A, and the critical time step \Delta t = 4/6*sqrt(3)A. Hence on a surface with varying size triangles, the timestep is limited by the area of the smallest triangles.

Approximate area. This force turns out to be an eigenvector of the mobility matrix with eigenvalue 2/A. Hence the velocity is four times that of the area normalization, and the critical time step four times shorter.

Back to top of Surface Evolver documentation.       Index.