diffu_fd2.qmd

## Diffusion with variable coefficient {#sec-diffu-varcoeff}

Diffusion in heterogeneous media normally implies a non-constant
diffusion coefficient $\alpha = \alpha (x)$.
A 1D diffusion model with such a variable diffusion coefficient reads

```{=latex}
\begin{alignat}{2}
\frac{\partial u}{\partial t} &=
\frac{\partial}{\partial x}\left( \alpha (x) \frac{\partial u}{\partial x}
\right) + f(x,t), \quad &x\in (0,L),\ t\in (0,T],
\\
u(x,0) &= I(x), \quad  &x\in [0,L],
\\
u(0,t) & = U_0, \quad  &t>0,
\\
u(L,t) & = U_L, \quad  &t>0.
\end{alignat}
```

A short form of the diffusion equation with variable coefficients is
$u_t = (\alpha u_x)_x + f$.

## Discretization {#sec-diffu-varcoeff-discr}

We can discretize the diffusion equation $u_t = (\alpha u_x)_x + f$ by a $\theta$-rule in time
and centered differences in space:
$$
\lbrack D_t u\rbrack^{n+\half}_i = \theta\lbrack D_x(\overline{\dfc}^x
D_x u) + f\rbrack^{n+1}_i +
(1-\theta)\lbrack D_x(\overline{\dfc}^x
D_x u) + f\rbrack^{n}_i\tp
$$
Written out, this becomes

\begin{align*}
\frac{u^{n+1}_i-u^{n}_i}{\Delta t} &=
\theta\frac{1}{\Delta x^2}
(\dfc_{i+\half}(u^{n+1}**{i+1} - u^{n+1}**{i})
- \dfc_{i-\half}(u^{n+1}**i - u^{n+1}**{i-1})) +\\
&\quad (1-\theta)\frac{1}{\Delta x^2}
(\dfc_{i+\half}(u^{n}**{i+1} - u^{n}**{i})
- \dfc_{i-\half}(u^{n}**i - u^{n}**{i-1})) +\\
&\quad \theta f_i^{n+1} + (1-\theta)f_i^{n},
\end{align*}
where, e.g., an arithmetic mean can to be used for $\dfc_{i+\half}$:
$$
\dfc_{i+\half} = \half(\dfc_i + \dfc_{i+1})\tp
$$
## Implementation {#sec-diffu-varcoeff-impl}

Suitable code for solving the discrete equations is very similar to
what we created for a constant $\dfc$.
Since the Fourier number has no meaning for varying
$\dfc$, we introduce a related parameter $D=\Delta t /\Delta x^2$.

```python
def solver_theta(I, a, L, Nx, D, T, theta=0.5, u_L=1, u_R=0,
                 user_action=None):
    x = linspace(0, L, Nx+1)   # mesh points in space
    dx = x[1] - x[0]
    dt = D*dx**2
    Nt = int(round(T/float(dt)))
    t = linspace(0, T, Nt+1)   # mesh points in time

    u   = zeros(Nx+1)   # solution array at t[n+1]
    u_n = zeros(Nx+1)   # solution at t[n]

    Dl = 0.5*D*theta
    Dr = 0.5*D*(1-theta)

    diagonal = zeros(Nx+1)
    lower    = zeros(Nx)
    upper    = zeros(Nx)
    b        = zeros(Nx+1)

    diagonal[1:-1] = 1 + Dl*(a[2:] + 2*a[1:-1] + a[:-2])
    lower[:-1] = -Dl*(a[1:-1] + a[:-2])
    upper[1:]  = -Dl*(a[2:] + a[1:-1])
    diagonal[0] = 1
    upper[0] = 0
    diagonal[Nx] = 1
    lower[-1] = 0

    A = scipy.sparse.diags(
        diagonals=[diagonal, lower, upper],
        offsets=[0, -1, 1],
        shape=(Nx+1, Nx+1),
        format='csr')

    for i in range(0,Nx+1):
        u_n[i] = I(x[i])

    if user_action is not None:
        user_action(u_n, x, t, 0)

    for n in range(0, Nt):
        b[1:-1] = u_n[1:-1] + Dr*(
            (a[2:] + a[1:-1])*(u_n[2:] - u_n[1:-1]) -
            (a[1:-1] + a[0:-2])*(u_n[1:-1] - u_n[:-2]))
        b[0]  = u_L(t[n+1])
        b[-1] = u_R(t[n+1])
        u[:] = scipy.sparse.linalg.spsolve(A, b)

        if user_action is not None:
            user_action(u, x, t, n+1)

        u_n, u = u, u_n
```

The code is found in the file [`diffu1D_vc.py`](https://github.com/devitocodes/devito_book/tree/main/src/diffu/diffu1D_vc.py).

## Stationary solution {#sec-diffu-varcoeff-stationary}

As $t\rightarrow\infty$, the solution of the
variable-coefficient diffusion problem
will approach
a stationary limit where $\partial u/\partial t=0$. The governing
equation is then
$$
\frac{d}{dx}\left(\alpha\frac{du}{dx}\right) =0,
$$ {#eq-diffu-fd2-pde-st}
with boundary conditions $u(0)=U_0$ and $u(L)=U_L$.
It is possible to obtain an exact solution of (@eq-diffu-fd2-pde-st)
for any $\alpha$. Integrating twice and applying the boundary conditions
to determine the integration constants gives
$$
u(x) = U_0 + (U_L-U_0)\frac{\int_0^x (\alpha(\xi))^{-1}d\xi}{\int_0^L (\alpha(\xi))^{-1}d\xi} \tp
$$ {#eq-diffu-fd2-pde-st-sol}

## Piecewise constant medium {#sec-diffu-varcoeff-piecewise}

Consider a medium built of $M$ layers. The layer boundaries
are denoted $b_0, \ldots, b_M$,
where $b_0=0$ and $b_M=L$.
If the layers potentially have different material properties, but
these properties are constant within each layer, we can express $\alpha$ as a
*piecewise constant function* according to
$$
\alpha (x) = \left\lbrace\begin{array}{ll}
\alpha_0,& b_0 \leq x < b_1,\\
\vdots &\\
\alpha_i,& b_i \leq x < b_{i+1},\\
\vdots &\\
\alpha_{M-1},& b_{M-1} \leq x \leq b_M.
\end{array}\right.
$$ {#eq-diffu-fd2-pde-st-pc-alpha}

The exact solution (@eq-diffu-fd2-pde-st-sol) in case of such a
piecewise constant $\alpha$ function is easy to derive. Assume that
$x$ is in the $m$-th layer: $x\in [b_m, b_{m+1}]$. In the integral
$\int_0^x (a(\xi))^{-1}d\xi$ we must integrate through the first
$m-1$ layers and then add the contribution from the remaining part
$x-b_m$ into the $m$-th layer:
$$
u(x) = U_0 + (U_L-U_0)
\frac{\sum_{j=0}^{m-1} (b_{j+1}-b_j)/\alpha(b_j) + (x-b_m)/\alpha(b_m)}{\sum_{j=0}^{M-1} (b_{j+1}-b_j)/\alpha(b_j)}
$$ {#eq-diffu-fd2-pde-st-sol-pc}

__Remark.__
It may sound strange to have a discontinuous $\alpha$ in a differential
equation where one is to differentiate, but a discontinuous $\alpha$
is compensated by a discontinuous $u_x$ such that $\alpha u_x$ is
continuous and therefore can be differentiated as $(\alpha u_x)_x$.

## Implementation of diffusion in a piecewise constant medium {#sec-diffu-varcoeff-impl-piecewise}

Programming with piecewise function definitions quickly becomes
cumbersome as the most naive approach is to test for which interval
$x$ lies, and then start evaluating a formula like
(@eq-diffu-fd2-pde-st-sol-pc). In Python, vectorized expressions may
help to speed up the computations.
The convenience classes `PiecewiseConstant` and
`IntegratedPiecewiseConstant` in the [`Heaviside`](https://github.com/devitocodes/devito_book/tree/main/src/diffu/Heaviside.py)
module were made to simplify programming with
functions like (@eq-diffu-fd2-pde-st-pc-alpha) and expressions like
(@eq-diffu-fd2-pde-st-sol-pc). These utilities not only represent
piecewise constant functions, but also *smoothed* versions of them
where the discontinuities can be smoothed out in a controlled fashion.

The `PiecewiseConstant` class is created by sending in the domain as a
2-tuple or 2-list and a `data` object describing the boundaries
$b_0,\ldots,b_M$ and the corresponding function values
$\alpha_0,\ldots,\alpha_{M-1}$. More precisely, `data` is a nested
list, where `data[i][0]` holds $b_i$ and `data[i][1]` holds the
corresponding value $\alpha_i$, for $i=0,\ldots,M-1$. Given $b_i$ and
$\alpha_i$ in arrays `b` and `a`, it is easy to fill out the nested
list `data`.

In our application, we want to represent $\alpha$ and $1/\alpha$
as piecewise constant functions, in addition to the $u(x)$ function
which involves the integrals of $1/\alpha$. A class creating the
functions we need and a method for evaluating $u$, can take the
form

```python
class SerialLayers:
    """
    b: coordinates of boundaries of layers, b[0] is left boundary
    and b[-1] is right boundary of the domain [0,L].
    a: values of the functions in each layer (len(a) = len(b)-1).
    U_0: u(x) value at left boundary x=0=b[0].
    U_L: u(x) value at right boundary x=L=b[0].
    """

    def __init__(self, a, b, U_0, U_L, eps=0):
        self.a, self.b = np.asarray(a), np.asarray(b)
        self.eps = eps  # smoothing parameter for smoothed a
        self.U_0, self.U_L = U_0, U_L

        a_data = [[bi, ai] for bi, ai in zip(self.b, self.a)]
        domain = [b[0], b[-1]]
        self.a_func = PiecewiseConstant(domain, a_data, eps)

        inv_a_data = [[bi, 1./ai] for bi, ai in zip(self.b, self.a)]
        self.inv_a_func = \
             PiecewiseConstant(domain, inv_a_data, eps)
        self.integral_of_inv_a_func = \
             IntegratedPiecewiseConstant(domain, inv_a_data, eps)
        self.inv_a_0L = self.integral_of_inv_a_func(b[-1])

    def __call__(self, x):
        solution = self.U_0 + (self.U_L-self.U_0)*\
                   self.integral_of_inv_a_func(x)/self.inv_a_0L
        return solution
```

A visualization method is also convenient to have. Below we plot $u(x)$
along with $\alpha (x)$ (which works well as long as $\max \alpha(x)$
is of the same size as $\max u = \max(U_0,U_L)$).

```python
class SerialLayers:
    ...

    def plot(self):
        x, y_a = self.a_func.plot()
        x = np.asarray(x); y_a = np.asarray(y_a)
        y_u = self.u_exact(x)
        import matplotlib.pyplot as plt
        plt.figure()
        plt.plot(x, y_u, 'b')
        plt.hold('on')  # Matlab style
        plt.plot(x, y_a, 'r')
        ymin = -0.1
        ymax = 1.2*max(y_u.max(), y_a.max())
        plt.axis([x[0], x[-1], ymin, ymax])
        plt.legend(['solution $u$', 'coefficient $a$'], loc='upper left')
        if self.eps > 0:
            plt.title('Smoothing eps: %s' % self.eps)
        plt.savefig('tmp.pdf')
        plt.savefig('tmp.png')
        plt.show()
```

Figure @fig-diffu-fd2-pde-st-sol-pc-fig1 shows the case where

```python
b = [0, 0.25, 0.5, 1]   # material boundaries
a = [0.2, 0.4, 4]       # material values
U_0 = 0.5;  U_L = 5     # boundary conditions
```

![Solution of the stationary diffusion equation corresponding to a piecewise constant diffusion coefficient.](fig/flow_in_layers_case1){#fig-diffu-fd2-pde-st-sol-pc-fig1 width="400px"}

By adding the `eps` parameter to the constructor of the `SerialLayers`
class, we can experiment with smoothed versions of $\alpha$ and see
the (small) impact on $u$. Figure @fig-diffu-fd2-pde-st-sol-pc-fig2
shows the result.

![Solution of the stationary diffusion equation corresponding to a *smoothed* piecewise constant diffusion coefficient.](fig/flow_in_layers_case1_eps){#fig-diffu-fd2-pde-st-sol-pc-fig2 width="400px"}

## Axi-symmetric diffusion {#sec-diffu-fd2-radial}

Suppose we have a diffusion process taking place in a straight tube
with radius $R$. We assume axi-symmetry such that $u$ is just a
function of $r$ and $t$, with $r$ being the radial distance from the center
axis of the tube to a point. With such axi-symmetry it is
advantageous to introduce *cylindrical coordinates* $r$, $\theta$, and
$z$, where $z$ is in the direction of the tube and $(r,\theta)$ are
polar coordinates in a cross section. Axi-symmetry means that all
quantities are independent of $\theta$. From the relations $x=\cos\theta$,
$y=\sin\theta$, and $z=z$, between Cartesian and cylindrical coordinates,
one can (with some effort) derive the diffusion equation in cylindrical
coordinates, which with axi-symmetry takes the form
$$
\frac{\partial u}{\partial t} = \frac{1}{r}\frac{\partial}{\partial r}
\left(r\dfc(r,z)\frac{\partial u}{\partial r}\right) + \frac{\partial}{\partial z}
\left(\alpha(r,z)\frac{\partial u}{\partial z}\right) + f(r,z,t)\tp
$$
Let us assume that $u$ does not change along the tube axis so it
suffices to compute variations in a cross section. Then $\partial u/\partial
z = 0$ and we have a 1D diffusion equation in the radial coordinate
$r$ and time $t$. In particular, we shall address the initial-boundary
value problem

$$
\frac{\partial u}{\partial t} = \frac{1}{r}\frac{\partial}{\partial r}
\left(r\dfc(r)\frac{\partial u}{\partial r}\right) + f(t), \quad r\in (0,R),\ t\in (0,T],
$$ {#eq-diffu-fd2-radial-PDE}
$$
\frac{\partial u}{\partial r}(0,t) = 0, \quad t\in (0,T],
$$ {#eq-diffu-fd2-radial-symmr0}
$$
u(R,t) = 0, \quad t\in (0,T],
$$ {#eq-diffu-fd2-radial-uR}
$$
u(r,0) = I(r), \quad r\in [0,R].
$$ {#eq-diffu-fd2-radial-initial}

The condition (@eq-diffu-fd2-radial-symmr0) is a necessary symmetry condition
at $r=0$, while (@eq-diffu-fd2-radial-uR) could be any Dirichlet
or Neumann condition (or Robin condition in case of cooling or heating).

The finite difference approximation will need the discretized version
of the PDE for $r=0$ (just as we use the PDE at the boundary when
implementing Neumann conditions). However, discretizing the PDE at
$r=0$ poses a problem because of the $1/r$ factor. We therefore need
to work out the PDE for discretization at $r=0$ with care.
Let us, for the case of constant $\dfc$, expand the spatial derivative term to
$$
\alpha\frac{\partial^2 u}{\partial r^2} + \alpha\frac{1}{r}\frac{\partial u}{\partial r}\tp
$$
The last term faces a difficulty at $r=0$, since it becomes a $0/0$ expression
caused by the symmetry condition at $r=0$.
However, L'Hospital's rule can be used:
$$
\lim_{r\rightarrow 0} \frac{1}{r}\frac{\partial u}{\partial r}
= \frac{\partial^2 u}{\partial r^2}\tp
$$
The PDE at $r=0$ therefore becomes
$$
\frac{\partial u}{\partial t} = 2\dfc\frac{\partial^2 u}{\partial r^2}
- f(t)\tp
$$ {#eq-diffu-fd2-radial-eq-PDEr0-aconst}
For a variable coefficient $\dfc(r)$ the expanded spatial derivative term reads
$$
\dfc(r)\frac{\partial^2 u}{\partial r^2} +
\frac{1}{r}(\dfc(r) + r\dfc'(r))\frac{\partial u}{\partial r}\tp
$$
We are interested in this expression for $r=0$. A necessary condition
for $u$ to be axi-symmetric is that all input data, including $\alpha$,
must also be axi-symmetric, implying that $\alpha'(0)=0$ (the second
term vanishes anyway because of $r=0$). The limit of interest is
$$
\lim_{r\rightarrow 0}
\frac{1}{r}\dfc(r)\frac{\partial u}{\partial r} =
\dfc(0)\frac{\partial^2 u}{\partial r^2}\tp
$$
The PDE at $r=0$ now looks like
$$
\frac{\partial u}{\partial t} = 2\dfc(0)
\frac{\partial^2 u}{\partial r^2}
- f(t),
$$ {#eq-diffu-fd2-radial-eq-PDEr0-avar}
so there is no essential difference between the constant coefficient
and variable coefficient cases.

The second-order derivative in (@eq-diffu-fd2-radial-eq-PDEr0-aconst)
and (@eq-diffu-fd2-radial-eq-PDEr0-avar)
is discretized in the usual way.
$$
2\dfc\frac{\partial^2}{\partial r^2}u(r_0,t_n) \approx
[2\dfc D_rD_r u]^n_0 =
2\dfc \frac{u^{n}_{1} - 2u^{n}**0 + u^n**{-1}}{\Delta r^2}\tp
$$
The fictitious value $u^n_{-1}$ can be eliminated using the discrete
symmetry condition
$$
[D_{2r} u =0]^n_0 \quad\Rightarrow\quad u^n_{-1} = u^n_1,
$$
which then gives the modified approximation to the term with the second-order derivative
of $u$ in $r$ at $r=0$:
$$
4\dfc \frac{u^{n}_{1} - u^{n}_0}{\Delta r^2}\tp
$$
The discretization of the term with the second-order derivative in $r$ at any
internal mesh point is straightforward:

\begin{align*}
\left[\frac{1}{r}\frac{\partial}{\partial r}
\left(r\dfc\frac{\partial u}{\partial r}\right)\right]_i^n
& \approx [r^{-1} D_r (r \dfc D_r u)]_i^n\\
&= \frac{1}{r_i}\frac{1}{\Delta r^2}\left(
r_{i+\half}\dfc_{i+\half}(u_{i+1}^n - u_i^n) - r_{i-\half}\dfc_{i-\half}(u_{i}^n - u_{i-1}^n)\right)\tp
\end{align*}

To complete the discretization, we need a scheme in time, but that can
be done as before and does not interfere with the discretization in space.

## Spherically-symmetric diffusion
### Discretization in spherical coordinates {#sec-diffu-fd2-spherical}
Let us now pose the problem from Section @sec-diffu-fd2-radial
in spherical coordinates, where $u$ only depends on the radial coordinate
$r$ and time $t$. That is, we have spherical symmetry.
For simplicity we restrict the diffusion coefficient $\dfc$ to be
a constant. The PDE reads
$$
\frac{\partial u}{\partial t} = \frac{\dfc}{r^\gamma}\frac{\partial}{\partial r}
\left(r^\gamma\frac{\partial u}{\partial r}\right) + f(t),
$$
for $r\in (0,R)$ and $t\in (0,T]$. The parameter $\gamma$ is 2 for
spherically-symmetric problems and 1 for axi-symmetric problems.
The boundary and initial conditions
have the same mathematical form as
in (@eq-diffu-fd2-radial-PDE)-(@eq-diffu-fd2-radial-initial).

Since the PDE in spherical coordinates has the same form as the PDE
in Section @sec-diffu-fd2-radial, just with the $\gamma$ parameter
being different, we can use the same discretization approach.
At the origin $r=0$ we get problems with the term
$$
\frac{\gamma}{r}\frac{\partial u}{\partial t},
$$
but L'Hospital's rule shows that this term equals $\gamma\partial^2 u/
\partial r^2$, and the PDE at $r=0$ becomes
$$
\frac{\partial u}{\partial t} = (\gamma+1)\dfc\frac{\partial^2 u}{\partial r^2}
- f(t)\tp
$$
The associated discrete form is then
$$
[D_t u = \half (\gamma+1)\dfc D_rD_r \overline{u}^t + \overline{f}^t]^{n+\frac{1}{2}}_i,
$$
for a Crank-Nicolson scheme.

### Discretization in Cartesian coordinates
The spherically-symmetric spatial derivative can be transformed to
the Cartesian counterpart by introducing
$$
v(r,t) = ru(r,t)\tp
$$
Inserting $u=v/r$ in
$$
\frac{1}{r^2}\frac{\partial}{\partial r}
\left(\dfc(r)r^2\frac{\partial u}{\partial r}\right),
$$
yields
$$
r\left(\frac{d \dfc}{dr}\frac{\partial v}{\partial r} +
\dfc\frac{\partial^2 v}{\partial r^2}\right) - \frac{d \dfc}{dr}v \tp
$$
The two terms in the parenthesis can be combined to
$$
r\frac{\partial}{\partial r}\left( \dfc\frac{\partial v}{\partial r}\right)\tp
$$
The PDE for $v$ takes the form
$$
\frac{\partial v}{\partial t} = \frac{\partial}{\partial r}\left( \dfc
\frac{\partial v}{\partial r}\right) - \frac{1}{r}\frac{d\dfc}{dr}v + rf(r,t),
\quad r\in (0,R),\ t\in (0,T]\tp
$$
For $\alpha$ constant we immediately realize that we can reuse a
solver in Cartesian coordinates to compute $v$. With variable $\alpha$,
a "reaction" term $v/r$ needs to be added to the Cartesian solver.
The boundary condition $\partial u/\partial r=0$ at $r=0$, implied
by symmetry, forces $v(0,t)=0$, because
$$
\frac{\partial u}{\partial r} = \frac{1}{r^2}\left(
r\frac{\partial v}{\partial r} - v\right) = 0,\quad r=0\tp
$$