% ex_reaction_diffusion.m — A time-dependent PDE on a surface.
%
% The complex Ginzburg-Landau equation
%
%     du/dt = delta (1 + i*alpha) Delta_Gamma u + u - (1 + i*beta) u |u|^2
%
% is integrated with a semi-implicit (backward-Euler) scheme: the
% diffusion term is treated implicitly, the reaction term explicitly. The
% diffusion operator is time-independent, so one surfaceop factorization
% is reused at every step. We use a sphere and a short run to stay quick.

mip load --install flatironinstitute/flatironinstitute/surfacefun

% Parameters.
dt = 0.1;
nsteps = 60;
alpha = 0;
beta = 1.5;
delta = 5e-3;

% Nonlinear (reaction) term.
N = @(u) u - (1 + beta*1i)*u .* (abs(u).^2);

% Mesh and the reusable implicit solver for (I - dt*delta*(1+i*alpha) Lap).
p = 8;
nref = 1;
dom = surfacemesh.sphere(p + 1, nref);
pdo = struct('lap', -dt*delta*(1 + alpha*1i), 'c', 1);
L = surfaceop(dom, pdo);
L.build();

% Random initial condition.
rng(1)
f = randnfun3(0.5, boundingbox(dom));
u = surfacefun(@(x, y, z) f(x, y, z), dom);

figure(1);
plot(real(u)), colorbar
title('Initial condition');

% Time stepping: each step is one solve with a fresh right-hand side.
for k = 1:nsteps
    L.rhs = u + dt*N(u);
    u = L.solve();
end

figure(2);
plot(real(u)), colorbar
title(sprintf('Re(u) after %d steps (t = %.1f)', nsteps, nsteps*dt));

fprintf('done: %d steps, ||Re(u)|| = %.4f\n', nsteps, norm(real(u)));