Function for calculating the Laplacian of the concentration in cylindrical coordinates via convolutions. More...

#include <stdio.h>
#include <stdlib.h>

Include dependency graph for convo.c:

Macros
#define	A(i, j) a[((i)*N+(j))]
	Retrieves the appropriate element from the 1D 'a' array given pseudo-2D indices i (z index) and j (r index), for the convolutions; uses column-major ordering.

#define	OUT(i, j) out[((i)*N+(j))]
	Retrieves the appropriate element from the 1D 'out' array given pseudo-2D indices i (z index) and j (r index); uses column-major ordering.

Functions
void	convolve3 (int M, int N, double a, double scale1, double scale2, double invr, double *out)
	Calculates updates to the concentration in a layer by applying the Laplacian in cylindrical coordinates; also scales by $\theta D_{free}$ . More...

Detailed Description

Function for calculating the Laplacian of the concentration in cylindrical coordinates via convolutions.

The change in the concentration with time depends on the Laplacian of the concentration. In cylindrical coordinates, the Laplacian of $ c $ is

$\nabla^2{c} = \frac{\partial^2 c}{\partial z^2} + \frac{\partial^2 c}{\partial r^2} + \frac{1}{r} \frac{\partial c} {\partial r} + \frac{1}{r^2} \frac{\partial^2 c}{\partial \phi^2}$

Because of circular symmetry, $\partial^2c/\partial \phi^2 = 0$ .

Note that in this program $\Delta z = \Delta r$ , which simplifies the kernels below.

The first 2 terms in the equation for the Laplacian are implemented numerically as if they were a 2D Laplacian in cartesian coordinates, by convolving the concentration array $ c $ with the Laplacian kernel

$\mathbf{L} = \left( \begin{array}{ccc} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0 \end{array} \right)$

and dividing by $\Delta r^2 \; (= \Delta z^2)$ (via a scaling factor).

The first derivative with respect to $ r $ in the third term is implemented numerically by convolving the $ c $ array with the derivative kernel

$\mathbf{D} = \left( \begin{array}{ccc} 0 & 0 & 0 \\ -1 & 0 & 1 \\ 0 & 0 & 0 \end{array} \right)$

and dividing by $2 \Delta r$ (via another scaling factor). Rows of the concentration array correspond to $ r $ and columns correspond to $ z $ .

In the third term, $ 1/r $ is a problem numerically for $ r = 0 $ . But by L'Hopital's Rule,

$\lim_{r \to 0} (\partial c/\partial r) / r = \partial ^2c/\partial r^2$

So for the case $ r = 0 $ we use

$\nabla^2{c} = \partial ^2c/\partial z^2 + 2 \partial ^2c/\partial r^2 \qquad (r = 0)$

which is implemented by convolving $ c $ with the kernel

$\mathbf{L0} = \left( \begin{array}{ccc} 0 & 1 & 0 \\ 2 & -6 & 2 \\ 0 & 1 & 0 \end{array} \right)$

Also the Laplacian is scaled by $\theta D_{free}$ via the scaling factors.

Author: David Lewis, CABI, NKI

Copyright: GNU Public License

Date: 2012-2013

Definition in file convo.c.

Function Documentation

void convolve3	(	int	M,
		int	N,
		double *	a,
		double	scale1,
		double	scale2,
		double *	invr,
		double *	out
	)

Calculates updates to the concentration in a layer by applying the Laplacian in cylindrical coordinates; also scales by $\theta D_{free}$ .

Calculates the update terms

$\mathbf{out} = s_1 \mathbf{L} * \mathbf{a} + \frac{s_2}{r} \mathbf{D} * \mathbf{a} \qquad (r \neq 0)$

or

$\mathbf{out} = s_1 \mathbf{L0} * \mathbf{a} \qquad (r = 0)$

where

$\mathbf{out}$ = output array

$ s_1, s_2 $ = scaling factors

$\mathbf{a}$ = input array (concentration)

$ * $ means convolution

$\mathbf{L}$ = 3x3 Laplace operator (see above)

$\mathbf{L0}$ = 3x3 Laplace operator modified for $ r = 0 $

$\mathbf{D}$ = 3x3 derivative operator for $ r $

Note that the index $ j = 1 $ corresponds to $ r = 0 $ .

Parameters

[in]	M	Number of columns of input matrix (z)
[in]	N	Number of rows of input matrix (r)
[in]	a	Input matrix (concentration, c)
[in]	scale1	Scaling factor #1
[in]	scale2	Scaling factor #2
[in]	invr	Vector of 1/r values
[out]	out	Output matrix

Definition at line 133 of file convo.c.

Macros

Functions

Detailed Description

Function Documentation