2. Tutorial

In this tutorial, we will attempt to solve a simple shape optimization problem, which will hopefully illustrate what this little code can do. The equations we work with are the two-phase static Stokes equations

\[\begin{split}\begin{cases} \nabla \cdot \mathbf{u}_\pm = 0, \\ \nabla \cdot \sigma_\pm = 0, \\ \mathbf{u}_+ = \mathbf{u}_\infty. \end{cases}\end{split}\]

The main usecase for these equations is to describe droplets suspended in an ambient fluid. At the interface, we pose the existence of a constant surface tension. This gives rise to the well-known Young-Laplace law. Specifically, we have the following jump conditions

\[\begin{split}\begin{cases} [\![ \mathbf{u} ]\!] = 0, \\ \displaystyle [\![ \mathbf{n} \cdot \sigma ]\!] = \frac{1}{\mathrm{Ca}} \kappa \mathbf{n}, \end{cases}\end{split}\]

where \(\kappa\) is the sum of the principal curvatures, \(\mathbf{n}\) is the normal exterior to the drop (here denoted by the \(-\)) and \(\mathrm{Ca}\) is the Capillary number.

Using these constraints, we propose to minimize the cost functional

\[\mathcal{J} = \frac{1}{2} \int_\Sigma |\mathbf{u} \cdot \mathbf{n} - u_d| \,\mathrm{d}S.\]

The desired velocity field \(u_d\) is solved on a desired geometry \(\mathbf{X}_d\), which we will see at the end. We now proceed to setting up this problem in code as follows.

% {{{ flow setup

% As a first step, we will initialize a simple extensional flow.

% The actual parameters don't matter that much for this example, but a unit
% viscosity ratio will make it a lot faster (no need to solve a system)

param.Ca = 0.1;
param.lambda = 1.0;

% We use a special function to get all the boundary conditions. There's
% quite a few that need to be setup, so this is faster. The ``'fn'`` key
% can also be a struct with all the boundary conditions we would like to
% have for custom cases.

bc = st_boundary_condition_options(...
    'param', param, ...
    'fn', 'FourRoller');

% Next, we also initialize a special struct for the constant capillary
% number Young-Laplace jump we are going to use. The jump is already defined
% in the boundary conditions, but this struct will also contain expressions
% for shape derivatives of this particular jump, so we don't have to worry
% about it.

bc.jump = st_jump_capillary('param', bc.param);

% }}}

% {{{ cost functional

% For the cost functional, we will attempt to match the normal velocity on
% a given geometry. Here, we just choose a nice UFO shaped geometry for fun.

% As for the jump, calling :func:`st_cost_normal_velocity` also defines all
% the required derivatives for this particular cost functional.

curve = @(xi) st_mesh_curve(xi, 'ufo', 'R', 1.0, 'k', 2);
cost = st_cost_normal_velocity(...
    'param', bc.param, ...
    'desired_curve', curve);

% Furthermore, it also defines all the variables that are required to compute
% the cost and all the derivatives and such. We need to aggregate these, so
% we can construct all the quadratures later on.
variables = unique([...
    cost.variables, ...
    cost.forward_variables, ...
    cost.adjoint_variables]);

% }}}

% {{{ geometry setup

% Now that the flow is setup, we define a simple geometry for our initial
% guess in the optimization algorithm. We use a 4th order scheme here.

% We also pass all the variables in the problem to :func:`st_geometry`, since
% it knows how to construct quadrature rules for all of them. The singular
% quadrature rules will have orders / nodes based on `nnodes` and some useful
% defaults. However, one can also pass a `s_quad` key with specific values
% for each one (this is passed to :func:`st_layer_quadrature`).

geometry = st_geometry(...
    'npanels', 32, ...
    'nnodes', 4, ...
    'curve', @(xi) st_mesh_curve(xi, 'circle', 'r', 1.0), ...
    's_singularity', variables);

x0 = geometry.x;
y0 = geometry.y;

% }}}

% {{{ optimization setup

% Next, we setup the parameters for the optimization routine. We use a
% gradient descent method with a very simple backtracking line search based
% on the Armijo rule. In general, shape optimization works best for small
% deformations, so small step sizes will be common and a fancy optimization
% method will nor fare substantially better.

% There are two parameters below that require some explanation. The `alpha_max`
% parameter gives and upper bound for the step size. We select it dynamically
% in such a way that an update maintains the axisymmetric assumption that all
% the points need to reside on the top half plane. Second, we also pass a
% `filter_strength`, which is used to further smooth the gradient during
% optimization.

cg = st_optim_options(...
    'display', true, ...
    'beta_update', 'flecherreeves', ...
    'maxit', 128, ...
    'gtol', 1.0e-6, ...
    'alpha_min', 1.0e-8, ...
    'alpha_max', @st_optim_geometry_alpha_max, ...
    'user', struct('filter_strength', 1.0e-2));

% }}}

% {{{ desired velocity field

% Now that we have a geometry, we will use the same parameters to compute
% the desired geometry in the cost functional. This will compute the points
% themselves and other variables (e.g. curvature and normal) on the desired
% curve. We will use this information to solve for the velocity on this
% curve and then use that as part of the optimization.

[xd, yd] = st_mesh_geometry_update(x0, y0, cost.desired_curve);
cost.x = xd;
cost.y = yd;

% To solve, we simply call :func:`st_repr_density_solution` and tell it what
% variables to compute and on what grid. This will give us the full velocity
% field on the surface of the desired geometry.

solution = st_repr_density_solution(xd, yd, bc, cost.variables);

cost.user.udx = st_dot(solution.u, xd.n);
cost.user.udy = st_dot(st_interpolate(xd, solution.u, yd.xi), yd.n);

% }}}

% {{{ optimization

% We now have done all the (super tedious) setup for our problem and can just
% pass it on to the optimizer and hope or the best!

% However, before doing that, we'll also set up a cache and have it output data
% at each iteratin for us. This is completely up to user control, so we just
% define an output function for the optimization.

cache = st_cache_create('tutorial', 'now', 'shape_optim');
cg.output = @(info) optim_output(info, cache, geometry, bc, cost);

% And on to the optimization!

[x, stats] = st_optim_geometry(cg, geometry, bc, cost);
[xf, yf] = st_mesh_geometry_update(x0, y0, x);

% }}}

% {{{ results

% We now have a solution, so lets plot it and see how well we did!

z0 = st_array_ax(x0.x);
zd = st_array_ax(xd.x);
zf = st_array_ax(xf.x);

fig = figure();
h1 = plot(real(z0), imag(z0), ':');
h2 = plot(real(zd), imag(zd), 'k-');
h3 = plot(real(zf), imag(zf), '--');
axis('equal');
grid('on');
xlabel('$x$');
ylabel('$\rho$');
legend([h1, h2, h3], {'Original', 'Desired', 'Final'}, ...
    'Location', 'EastOutside');
st_print(fig, 'tutorial_shape_optimization.png')

% }}}

function optim_output(info, cache, geometry, bc, cost)  %#ok
    % TODO
end