10. Utils#

10.1. Optimization#

st_optim_cg(varargin)#

Nonlinear conjugate gradient.

The algorithm uses usual methods (for the notation see [Hager2006]). The line search is a naive backtracking algorithm that just check the Armijo conditions, not the full Wolfe conditions.

Parameters:: cg – a struct containing the problem definitions.
Returns:: an array x with the solution.
Returns:: a struct with information about the iterations.

st_optim_options(varargin)#

Default options for st_optim_cg().

Required fields:

x0: initial guess.
[f, df] = fn(x): function handle that evaluates at a given x. If two output arguments are requested, the second one should be the gradient at x, of the same size and shape as x, i.e. \(f: \mathbb{R}^n \to \mathbb{R}\).

Optional fields include

display: display per-iteration information, e.g. function values.
maxiter: maximum number of iterations.
restart_maxit: maximum number of iterations before a restart. A restart resets the descent direction to the one given by the gradient. This can also happen automatically if the descent direction points in an opposite direction to the gradient, i.e. it is no longer a descent direction.
ftol: relative tolerance between consecutive values of function evaluations.
xtol: relative tolerance between consecutive values of the argument.
gtol: relative tolerance for the gradient.
beta: method used to compute the cg update parameter.
hz_eta: threshold in Hager-Zhang algorithm, see Equation 7.1 in [Hager2006].
hz_theta: weight in Hager-Zhang algorithm, see Equation 7.1 in [Hager2006].
alpha_max: maximum allowable step size. Can also be a user defined function handle that takes (x, d) as arguments.
alpha_min: minimum allowable step size, which can also be a user defined function handle.
armijo_beta: parameter in the Armijo condition.
armijo_ratio: decrease in alpha at each iteration.
output(o): user-defined function handle. Takes a struct with the fields x, f, d, iteration.
outputfreq: frequency of output per number of iterations.
project(x, d): user-defined function that can be used to project the result of an update to an allowable region.

Additional parameters, that remain unchecked, can be passed in a separate 'user' field. The default values can be obtained by calling this function without any arguments.

The CG update methods that are implemented are (see [Hager2006]):

'Steepest': \(\beta \equiv 0\).
'HestenesStiefl'.
'FletcherReeves'.
'PolakRibiere'.
'HagerZhang': like 'CGDescent' with \(\eta = \infty\) and \(\theta = 2\).
'CGDescent': uses the parameters hz_eta and hz_theta. This is not well-named, since it is not the algorithm from [Hager2006].

Returns:: a struct with default and updated options.

st_optim_geometry_alpha_max(x, d, beta)#

Compute maximum CG step size for shape optimization.

We have two heuristic requirements for the step size:

first, it should be sufficiently small so that the resulting shape deformation is small.
second, due to axisymmetry we need to ensure that \(\rho > 0\). This imposes the condition that

\[\rho + \alpha d_\rho > 0 \implies \alpha < \min -\rho_i / d_{\rho, i}.\]

Parameters:

x – interface coordinates.
d – descent direction.
beta (default 1.0.) – additional factor to scale alpha.

Returns:

a maximum step size.

10.2. Elliptic Integrals#

st_ellipke(x, x0)#

Compute elliptic integrals.

Matches the interface and results of the default MATLAB ellipke, but is a lot faster (and sufficient for our needs). The function essentially computes:

K2 = 4 * Y * Y0 / ((X - X0)^2 + (Y + Y0)^2);
[F, E] = ellipke(K2);

Parameters:

x – source point.
x0 – target point.

Returns:

a tuple (F, E) of elliptic integrals of the first and second kind.

Seealso:

the 05_stokes/drop_ax/ell_int.m function in the BEMLIB sources.

st_ellint_order3(x, x0, F, E)#

Compute elliptic integrals in the Stokeslet.

Parameters:

x – source point.
x0 – target point.
F – elliptic integral of the first kind, see st_ellipke().
E – elliptic integral of the second kind.

Returns:

elliptic integrals described in Section 2.4 of [Pozrikidis1992] used to compute the Stokeslet. Namely, it computes \(I_{10}, I_{11}, I_{30}, I_{31}\) and \(I_{32}\).

st_ellint_order5(x, x0, F, E)#

Compute elliptic integrals in the Stresslet.

Parameters:

x – source point.
x0 – target point.
F – elliptic integral of the first kind, see st_ellipke().
E – elliptic integral of the second kind.

Returns:

elliptic integrals described in Section 2.4 of [Pozrikidis1992] used to compute the Stresslet. Namely, it computes \(I_{50}, I_{51}, I_{52}\) and \(I_{53}\).

st_ellint_order7(x, x0, F, E)#

Compute elliptic integrals for the Stresslet derivative.

Parameters:

x – source point.
x0 – target point.
F – elliptic integral of the first kind, see st_ellipke().
E – elliptic integral of the second kind.

Returns:

elliptic integrals described in Section 2.4 of [Pozrikidis1992] used to compute the Stresslet derivative. Namely, it computes \(I_{70}, I_{71}, I_{72}, I_{73}\) and \(I_{54}\).

10.3. Struct Enhancements#

st_struct_field(s, name, default_value)#

Get the field from a struct.

This functions similarly to Python’s getter on dictionaries.

Parameters:

s – struct.
name – field name.
default_value – default value, if the field does not exist.

Returns:

field value.

st_struct_hash(varargin)#

Hash a given struct.

This function does not hash all the contents of a struct. We just use the scalar quantities and fieldnames. There is also no randomness, so the same struct will hash to the same value if no values were modified.

Parameters:: varargin – any number of structs.
Returns:: a string of the hexadecimal representation of the MD5 hash.

st_struct_merge(s, s0, overwrite)#

Merge fields from two structs.

This functions similarly to Python’s update on dictionaries when overwrite is true.

Parameters:

s – struct.
s0 – other struct.
overwrite (default false) – if true, the fields in s are overwritten with those in s0.

Returns:

new struct with fields in the union of the two structs.

st_struct_update(s, s0)#

Update a struct with fields from a another struct.

This function is function is meant to only update the existing fields in s with fields from s0. Not matching fields are simply ignored.

This functions similarly to Python’s update on dictionaries.

Parameters:

s – struct.
s0 – other struct.

Returns:

new struct with updated fields.

st_struct_varargin(varargin)#

Convert a cell array to a struct.

The first argument can be a struct, whose fields will be overwritten and the remaining parameters must be pairs of ("Name", value), in the usual MATLAB way. The same value can appear multiple times, but the last occurrence will overwrite any previous ones.

Returns:: a struct with appropriate fields.

10.4. Misc#

st_print(fig, filename, varargin)#

Wrapper around the default print command.

Default values can be set by the global variables:

PrintResolution
PrintFormat

Parameters:

fig – figure handle.
filename – filename.
varargin – additional arguments.

st_sym_handle(f, x, force_simplify)#

Wrapper around matlabFunction.

This function creates function handles of the form f(z) where z is either the source or target geometry information. This function constructs two types of handles

if x contains two symbolic variables, we assume they are the \((x, y)\) coordinates of the geometry.
if x is a single symbolic variable, we assume it refers to the parametrization \(\xi\) of the curve.

Other cases should be handled by directly calling matlabFunction as

f = simplify(f, 'Steps', force_simplify);
fn = matlabFunction(f, 'Vars', x);

Parameters:

f – symbolic expression for the function.
x – symbolic function variables.
force_simplify (default 0) – simplify the given expressions. This is meant to be an integer denoting the number of simplification steps (see sym/simplify).

Returns:

a function handle that takes the geometry as an argument.

st_ax_filter(f, method, p)#

Apply a filter to some data.

There are a few implemented filters at the moment. They are

'fft': smoothing by convolving with a Gaussian in Fourier space. Note that in this case, f must be periodic. Its smoothing parameter is \(p\), with a smaller \(p\) denoting less smoothing.
'lele': spectral-like filter based on [Lele1992] formulas (C.2.4) and (C.2.11). The strength of the smoothing is determined by the number of smoothing passes \(p\).
'lsq': least-squares method from [Eilers2003]. The parameter determines the smoothing strength.

Parameters:

f – real or complex array.
method – smoothing method.
p – smoothing parameter.

Returns:

filtered data.

st_fftfreq(m, dx)#

Compute FFT sample frequencies like numpy.

Parameters:

m – window length.
dx – domain size.

Returns:

array containing sample frequencies.

st_dot(x, y)#

Compute element-wise dot product of two complex arrays.

We represent a 2D real array \((u, v)\) as a complex number \(u + \imath v\). To compute the element-wise dot product of two such complex arrays, we have:

\[x \cdot y \triangleq \Re\{x\} \cdot \Re\{y\} + \Im\{x\} \cdot \Im\{y\}\]

Parameters:

x – complex array.
y – complex array.

Returns:

element-wise dot product (real array).

st_array_stack(w)#

Transform a complex array into a stacked real array.

We represent two-dimensional vector fields as complex one component arrays. This function can be used to stack the real and imaginary component so as to make a real array of size \(2 n\) of all the DOFs.

Note that if a real array is passed in, the imaginary component is assumed to be zero and the stacking proceeds in the same way.

Parameters:

w – complex representation of a vector field.
shape (If not given, the result will be a row or column vector the same as w.) – shape of the resulting array.

Returns:

stacked array.

st_array_unstack(v)#

Transform a real array into a complex array half the size.

This is a counterpart to st_array_stack() and satisfies

st_array_unstack(st_array_stack(x)) == x;

Parameters:: v – real array.
Returns:: a complex array.

st_array_ax(v, cclosed)#

Complete an axisymmetric array.

We represent an axisymmetric scalar or vector field only by its components in the top half plane. This function returns an array in the whole \(x-y\) plane by flipping the components accordingly. Mostly equivalent to

w = [v, fliplr(conj(v:2:end - 1))];

Parameters:

v – array representing an axisymmetric quantity.
cclosed – if true, the first and last element will be the same.

Returns:

complete vector field.

st_cache_create(datadir, cache_id, prefix, cachedir, overwrite)#

Simple structure for handling caches and general checkpointing.

If datadir is not provided, the constructed directory has the form:

$PWD/{cachedir}/{prefix}_{cache_id}

The returned struct has attached a few file generation routines:

filename = filename(cache, suffix, ext) can be used to generate filename in the cache directory with a canonical name.
[next, prev] = filenext(cache, suffix, ext) is similar to filename, but generates files with {suffix}_%05d in the cache directory and returns the next file that does not exist on disk.
[cache, filename] = fileseries(cache, suffix, ext) is similar to filenext, but generates index filenames without checking on disk. All the generated files can also be accessed in the file_series_cache field. Only one series can be active at the time and must have a matching suffix.

Parameters:

datadir (default []) – a full directory name for the cache; if empty, it gets created from the remaining data.
cache_id (default 'now') – a unique identifier for the cache directory; the special value 'now' creates a timestamp for the identifier.
prefix (default 'cache') – a non-unique prefix for the cache subfolder.
cachedir – name of the cache directory.

Returns:

a struct containing cache details.

st_util_timeit(fn, nruns, alpha)#

Time a given function.

Parameters:

fn – function handle with no arguments.
nruns (default \(128\)) – number of runs.
alpha (default \(0.95\)) – confidence level.

Returns:

a tuple (t_avg, t_std, t_total), with the average time, confidence interval and total time (all in seconds). The confidence interval is computed under the assumption that the timing samples have a Gaussian distribution.

st_util_eoc(eoc, capacity, window, verbose)#

Compute step by step estimated order of convergence (EOC).

The returned struct will have the following information:

index: current in the error array, i.e. values after this are not initialized.
errors: error arrays for each variable.
orders: order slopes with a window of 1 for each variable.
h: array containing grid sizes.
npoints: array containing number of grid points (optional).

Intended use is as follows:

npoints = [8, 16, 32, 64];
eoc = st_util_eoc({'var1', 'var2'}, length(npoints));

for n = 1:length(npoints)
    % ... perform computations ...

    % update errors with latest results
    eoc.errors.var1(n) = norm(var1 - var1_exact);
    eoc.errors.var2(n) = norm(var2 - var2_exact);
    eoc.h(n) = hmax;

    % update and print eoc
    eoc = st_util_eoc(eoc);
end

Parameters:

eoc – struct containing information about previous errors.
capacity – excepted number of values, only used at initialization.
window – window to compute eoc.
verbose – if true, also print errors at each call.

Returns:

updated struct.

st_util_eoc_show(eoc, filename, varargin)#

Plot errors and print orders of convergence

Parameters:

filename – name of output file.
errors – a struct with fields containing errors arrays.
h – mesh sizes corresponding to each errors level.

Returns:

convergence orders.

st_util_latex_table(data, fid)#

Print a table in latex format.

The printed table depends on the booktabs package. The data argument contains a set of fields, that corresponds to columns in the resulting table. Each field can be an array or a struct with the following fields:

data: data array to be written out.
header: column name.
fmt: formatting string for the column. By default, numbers are printed in a fixed precision and strings printed as is.

Not all columns need to have the same number of rows and can contain any type of information. The only requirement is that each column contains the same type of data, i.e. they are not cell arrays.

The data struct can also contain a fmt field that describes the column layout, passed as an argument to the tabular environment.

Parameters:

data – structured data to write out.
fid (default stdout.) – file identifier, from fopen.

Utils

Contents

10. Utils#

10.1. Optimization#

10.2. Elliptic Integrals#

10.3. Struct Enhancements#

10.4. Misc#