Introduction

In this example, we use the Helmholtz solver of Ballfun, available through the helmholtz command, to solve the advection-diffusion equation in the unit ball. We also use some of the vector calculus and visualization capabilities of Ballfun.

Advection-diffusion in the ball

The advection-diffusion equation in the ball is given by $$ \frac{\partial c}{\partial t}=D\nabla^2c-v\cdot\nabla c, $$ where $D$ is the diffusion coefficient and $v$ is a divergence-free vector field.

We choose $D=1/5000$ and $v = \nabla\times[ze^{-5(x^2+y^2+z^2)}(x,y,z)]$ to satisfy the no-slip condition $v\cdot \vec{n}=0$.

w = ballfunv( @(x,y,z) z.*exp(-5*(x.^2+y.^2+z.^2)).*x,...
              @(x,y,z) z.*exp(-5*(x.^2+y.^2+z.^2)).*y,...
              @(x,y,z) z.*exp(-5*(x.^2+y.^2+z.^2)).*z);
v = curl( w );
quiver(v, 4, 'numpts',30), axis('off'), colorbar

We verify that $v$ is divergence-free:

norm(div(v))
ans =
     1.269097204341824e-17

The vector field $v$ also satisfies the no-slip boundary condition $v\cdot\vec{n}=0$, as shown by the following command:

vn = dot(v(1,:,:,'spherical'),spherefunv.unormal);
norm(vn)
ans =
     1.296776305849438e-33

Vizualisation of functions

We impose the initial condition $c = -xe^{-5(x^2+y^2+z^2)}$ to solve the advection-diffusion equation.

c = ballfun(@(x,y,z) -x.*exp(-5*(x.^2+y.^2+z.^2)));

The function $c$ can be vizualised by the different plotting commands in Ballfun;

subplot(2,2,1)
plot(c), caxis([-0.19 0.19]), axis off
title("Plot")

subplot(2,2,2)
slice(c), caxis([-0.19 0.19]), axis off
title("Slice")

subplot(2,2,3)
plot(c, 'WedgeAz'), caxis([-0.19 0.19]), axis off
title("WedgeAz")

subplot(2,2,4)
plot(c, 'WedgePol'), caxis([-0.19 0.19]), axis off
title("WedgePol")

Time discretization

Now we solve the advection-diffusion equation numerically using the implicit-explicit order 1 backward differentiation time-stepping scheme (IMEX-BDF1). This yields a Helmholtz equation at each time step: $$ \nabla^2c^{n+1}+K^2c^{n+1}=K^2c^n+\frac{1}{D}v\cdot\nabla c^n,\quad \left.\frac{\partial c}{\partial \vec{n}}\right|_{\partial B(0,1)} = 0, $$ where $c_n$ denotes the solution at time $t = n\Delta t$, $\Delta t = 0.1$ is the time step, and $K^2 = -1/(D\Delta t)$. This equation can be solved by using the Ballfun command helmholtz.

The following code solves the advection-diffusion numerically to time $t=15$ and plots the solution $c$ at different times.

D = 1/5000;                                     % Diffusion constant
dt = 0.1;                                      % Time step
K = 1i*sqrt(1/(dt*D));                          % Helmholtz frequency
T = 15;                                         % Stopping time
nsteps = ceil(T/dt);                            % Number of time steps
m = 100;                                        % Spatial discretization

for n = 0:nsteps
    if mod(n,50) == 0
        clf, slice(c), caxis([-0.2,0.2])
        title(sprintf('Time %d',n*dt)), colorbar, axis('off'), snapnow
    end

    rhs = K^2*c+dot(v,grad(c))/D;
    c = helmholtz(rhs, K, @(x,y,z)0, m, 'neumann');
end