# Diffusion in heterogeneous media
<div id="diffu:varcoeff"></div>

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

<!-- Equation labels as ordinary links -->
<div id="diffu:pde2"></div>

$$
\begin{equation}
\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],
\label{diffu:pde2} \tag{1}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="diffu:pde2:ic:u"></div>

$$
\begin{equation} 
u(x,0) = I(x), \quad  x\in [0,L],
\label{diffu:pde2:ic:u} \tag{2}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="diffu:pde2:bc:0"></div>

$$
\begin{equation} 
u(0,t)  = U_0, \quad  t>0,
\label{diffu:pde2:bc:0} \tag{3}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="diffu:pde2:bc:L"></div>

$$
\begin{equation} 
u(L,t)  = U_L, \quad  t>0.
\label{diffu:pde2:bc:L} \tag{4}
\end{equation}
$$

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

## Discretization
<div id="diffu:varcoeff:discr"></div>

We can discretize ([1](#diffu:pde2)) by a $\theta$-rule in time
and centered differences in space:

$$
\lbrack D_t u\rbrack^{n+\frac{1}{2}}_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\thinspace .
$$

Written out, this becomes

$$
\begin{align*}
\frac{u^{n+1}_i-u^{n}_i}{\Delta t} &=
\theta\frac{1}{\Delta x^2}
(\dfc_{i+\frac{1}{2}}(u^{n+1}_{i+1} - u^{n+1}_{i})
- \dfc_{i-\frac{1}{2}}(u^{n+1}_i - u^{n+1}_{i-1})) +\\
&\quad (1-\theta)\frac{1}{\Delta x^2}
(\dfc_{i+\frac{1}{2}}(u^{n}_{i+1} - u^{n}_{i})
- \dfc_{i-\frac{1}{2}}(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+\frac{1}{2}}$:

$$
\dfc_{i+\frac{1}{2}} = \frac{1}{2}(\dfc_i + \dfc_{i+1})\thinspace .
$$

## Implementation
<div id="diffu:varcoeff:impl"></div>

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)

    # Representation of sparse matrix and right-hand side
    diagonal = zeros(Nx+1)
    lower    = zeros(Nx)
    upper    = zeros(Nx)
    b        = zeros(Nx+1)

    # Precompute sparse matrix (scipy format)
    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])
    # Insert boundary conditions
    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')

    # Set initial condition
    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)

    # Time loop
    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]))
        # Boundary conditions
        b[0]  = u_L(t[n+1])
        b[-1] = u_R(t[n+1])
        # Solve
        u[:] = scipy.sparse.linalg.spsolve(A, b)

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

        # Switch variables before next step
        u_n, u = u, u_n
```

The code is found in the file [`diffu1D_vc.py`](${src_diffu}/diffu1D_vc.py).

## Stationary solution
<div id="diffu:varcoeff:stationary"></div>

As $t\rightarrow\infty$, the solution of the
problem ([1](#diffu:pde2))-([4](#diffu:pde2:bc:L))
will approach
a stationary limit where $\partial u/\partial t=0$. The governing
equation is then

<!-- Equation labels as ordinary links -->
<div id="diffu:fd2:pde:st"></div>

$$
\begin{equation}
\frac{d}{dx}\left(\alpha\frac{du}{dx}\right) =0,
\label{diffu:fd2:pde:st} \tag{5}
\end{equation}
$$

with boundary conditions $u(0)=U_0$ and $u(L)=U_L$.
mathcal{I}_t is possible to obtain an exact solution of ([5](#diffu:fd2:pde:st))
for any $\alpha$. Integrating twice and applying the boundary conditions
to determine the integration constants gives

<!-- Equation labels as ordinary links -->
<div id="diffu:fd2:pde:st:sol"></div>

$$
\begin{equation}
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}
\thinspace .
\label{diffu:fd2:pde:st:sol} \tag{6}
\end{equation}
$$

## Piecewise constant medium
<div id="diffu:varcoeff:piecewise"></div>

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

<!-- Equation labels as ordinary links -->
<div id="diffu:fd2:pde:st:pc:alpha"></div>

$$
\begin{equation}
\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.
\end{equation}
\label{diffu:fd2:pde:st:pc:alpha} \tag{7}
$$

The exact solution ([6](#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:

<!-- Equation labels as ordinary links -->
<div id="diffu:fd2:pde:st:sol:pc"></div>

$$
\begin{equation}
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)}
\label{diffu:fd2:pde:st:sol:pc} \tag{8}
\end{equation}
$$

**Remark.**
mathcal{I}_t 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
<div id="diffu:varcoeff:impl:piecewise"></div>

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
([8](#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`](${src_diffu}/Heaviside.py)
module were made to simplify programming with
functions like ([7](#diffu:fd2:pde:st:pc:alpha)) and expressions like
([8](#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 = 1/a is needed in formulas
        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)
        # Denominator in the exact formula is constant
        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
%matplotlib inline

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](#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
```

<!-- dom:FIGURE: [fig-diffu/flow_in_layers_case1.png, width=400 frac=1.0] Solution of the stationary diffusion equation corresponding to a piecewise constant diffusion coefficient. <div id="diffu:fd2:pde:st:sol:pc:fig1"></div> -->
<!-- begin figure -->
<div id="diffu:fd2:pde:st:sol:pc:fig1"></div>

<p>Solution of the stationary diffusion equation corresponding to a piecewise constant diffusion coefficient.</p>


<!-- end figure -->


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](#diffu:fd2:pde:st:sol:pc:fig2)
shows the result.

<!-- dom:FIGURE: [fig-diffu/flow_in_layers_case1_eps.png, width=400 frac=1.0] Solution of the stationary diffusion equation corresponding to a *smoothed* piecewise constant diffusion coefficient. <div id="diffu:fd2:pde:st:sol:pc:fig2"></div> -->
<!-- begin figure -->
<div id="diffu:fd2:pde:st:sol:pc:fig2"></div>

<p>Solution of the stationary diffusion equation corresponding to a <em>smoothed</em> piecewise constant diffusion coefficient.</p>


<!-- end figure -->


## Axi-symmetric diffusion
<div id="diffu:fd2:radial"></div>


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)\thinspace .
$$

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

<!-- Equation labels as ordinary links -->
<div id="diffu:fd2:radial:PDE"></div>

$$
\begin{equation}
\frac{\partial u}{\partial t} = \frac{1}{r}\frac{\partial}{\partial r}
\left(r\dfc(r)\frac{\partial u}{\partial r}\right) + f(t),  r\in (0,R),\ t\in (0,T],
\label{diffu:fd2:radial:PDE} \tag{9}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="diffu:fd2:radial:symmr0"></div>

$$
\begin{equation} 
\frac{\partial u}{\partial r}(0,t) = 0,  t\in (0,T],
\label{diffu:fd2:radial:symmr0} \tag{10}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="diffu:fd2:radial:uR"></div>

$$
\begin{equation} 
u(R,t) = 0,  t\in (0,T],
\label{diffu:fd2:radial:uR} \tag{11}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="diffu:fd2:radial:initial"></div>

$$
\begin{equation} 
u(r,0) = I(r),  r\in [0,R].
\label{diffu:fd2:radial:initial} \tag{12}
\end{equation}
$$

The condition ([10](#diffu:fd2:radial:symmr0)) is a necessary symmetry condition
at $r=0$, while ([11](#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}\thinspace .
$$

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'Hosptial's rule can be used:

$$
\lim_{r\rightarrow 0} \frac{1}{r}\frac{\partial u}{\partial r}
= \frac{\partial^2 u}{\partial r^2}\thinspace .
$$

The PDE at $r=0$ therefore becomes

<!-- Equation labels as ordinary links -->
<div id="diffu:fd2:radial:eq_PDEr0:aconst"></div>

$$
\begin{equation}
\frac{\partial u}{\partial t} = 2\dfc\frac{\partial^2 u}{\partial r^2}
+ f(t)\thinspace .
\label{diffu:fd2:radial:eq_PDEr0:aconst} \tag{13}
\end{equation}
$$

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}\thinspace .
$$

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}\thinspace .
$$

The PDE at $r=0$ now looks like

<!-- Equation labels as ordinary links -->
<div id="diffu:fd2:radial:eq_PDEr0:avar"></div>

$$
\begin{equation}
\frac{\partial u}{\partial t} = 2\dfc(0)
\frac{\partial^2 u}{\partial r^2}
+ f(t),
\label{diffu:fd2:radial:eq_PDEr0:avar} \tag{14}
\end{equation}
$$

so there is no essential difference between the constant coefficient
and variable coefficient cases.

The second-order derivative in ([13](#diffu:fd2:radial:eq_PDEr0:aconst))
and ([14](#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}\thinspace .
$$

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$:

<!-- Equation labels as ordinary links -->
<div id="_auto1"></div>

$$
\begin{equation}
4\dfc \frac{u^{n}_{1} - u^{n}_0}{\Delta r^2}\thinspace .
\label{_auto1} \tag{15}
\end{equation}
$$

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+\frac{1}{2}}\dfc_{i+\frac{1}{2}}(u_{i+1}^n - u_i^n) - r_{i-\frac{1}{2}}\dfc_{i-\frac{1}{2}}(u_{i}^n - u_{i-1}^n)\right)\thinspace .
\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
<div id="diffu:fd2:spherical"></div>

### Discretization in spherical coordinates

Let us now pose the problem from the section [Axi-symmetric diffusion](#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

<!-- Equation labels as ordinary links -->
<div id="_auto2"></div>

$$
\begin{equation}
\frac{\partial u}{\partial t} = \frac{\dfc}{r^\gamma}\frac{\partial}{\partial r}
\left(r^\gamma\frac{\partial u}{\partial r}\right) + f(t),
\label{_auto2} \tag{16}
\end{equation}
$$

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 ([9](#diffu:fd2:radial:PDE))-([12](#diffu:fd2:radial:initial)).

Since the PDE in spherical coordinates has the same form as the PDE
in the section [Axi-symmetric diffusion](#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'Hosptial's rule shows that this term equals $\gamma\partial^2 u/
\partial r^2$, and the PDE at $r=0$ becomes

<!-- Equation labels as ordinary links -->
<div id="_auto3"></div>

$$
\begin{equation}
\frac{\partial u}{\partial t} = (\gamma+1)\dfc\frac{\partial^2 u}{\partial r^2}
+ f(t)\thinspace .
\label{_auto3} \tag{17}
\end{equation}
$$

The associated discrete form is then

<!-- Equation labels as ordinary links -->
<div id="_auto4"></div>

$$
\begin{equation}
[D_t u = \frac{1}{2} (\gamma+1)\dfc D_rD_r \overline{u}^t + \overline{f}^t]^{n+\frac{1}{2}}_i,
\label{_auto4} \tag{18}
\end{equation}
$$

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)\thinspace .
$$

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
\thinspace .
$$

The two terms in the parenthesis can be combined to

$$
r\frac{\partial}{\partial r}\left( \dfc\frac{\partial v}{\partial r}\right)\thinspace .
$$

The PDE for $v$ takes the form

<!-- Equation labels as ordinary links -->
<div id="_auto5"></div>

$$
\begin{equation}
\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]\thinspace .
\label{_auto5} \tag{19}
\end{equation}
$$

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\thinspace .
$$

2aa10fcc098e235d7026d09b14f330dd0349edc8
ipynb
Jupyter Notebook
fdm-devito-notebooks/03_diffu/diffu_fd2.ipynb
devitocodes/devito_book
30405c3d440a1f89df69594fd0704f69650c1ded
CC-BY-4.0
2020-07-17T13:19:15.000Z
2021-03-27T05:21:09.000Z
fdm-devito-notebooks/03_diffu/diffu_fd2.ipynb
devitocodes/devito_book
30405c3d440a1f89df69594fd0704f69650c1ded
CC-BY-4.0
2020-07-14T15:38:52.000Z
2020-09-25T11:54:59.000Z
fdm-devito-notebooks/03_diffu/diffu_fd2.ipynb
devitocodes/devito_book
30405c3d440a1f89df69594fd0704f69650c1ded
CC-BY-4.0
2021-03-27T05:21:14.000Z
2021-03-27T05:21:14.000Z
Qwen/Qwen-72B
1. YES
2. YES