Diffusion Equation
Parabolic Equations
PDEs are traditionally classified into three broad categories: elliptic, parabolic and hyperbolic equations. However, from the numerics point of view, PDEs divide into static equations (elliptic) and time dependent equations (parabolic & hyperbolic).
The hyperbolic and parabolic equations represent initial value problems. The diffusion equation (parabolic)
\[ \frac{\delta u}{\delta t} = -\frac{\delta}{\delta x} \left( D \frac{\delta u}{\delta x} \right), \]
where $D$ is the diffusion coefficient, are such that we ask for what is the value of the field (wave) $u(x,t)$ at a later time $t$ knowing the field at an initial time $t=0$ and subject to some specific boundary conditions at all times.
The domain may be 1D, 2D, or 3D. This equation might represent a heat diffusion problem in which the temperature within the domain is sought at all times, starting from a given temperature field at the initial time, and preserving the temperature at the boundary of the domain. In addition, a source term $S(x,t)$ could be added to the right hand side of the equation. If $u(x,t)$ represents a temperature field, then $S(x,t)$ would represent a heat source.
In this case, the domain decomposition applies to the space variables. The time evolution is essentially sequential, $u(t + d t, x)$ is obtained from $u(t,x)$ which can not be parallelized. This is also true in obtaining a solution of ODEs.
The parallelization of the solution of a PDE of parabolic or hyperbolic PDEs is obtained through domain decomposition. However, in this case, the domain decomposition applies only to the space variables. The time evolution is essentially sequential, $u(t + d t, x)$ is obtained from $u(t,x)$ which can not be parallelized.
Model
We will model heat diffusion through a 2D plate. The parabolic PDE to solve is
\[ \frac{\delta u(x,y,t)}{\delta t} = \frac{\delta^2 u(x,y,t)}{\delta x^2} + \frac{\delta^2 u(x,y,t)}{\delta y^2} + S(x,y,t) \]
Dirichlet boundary conditions are assumed, the temperature being fixed at the top and bottom of the plate, $u_\text{top}$ and $u_\text{bot}$, and on the left and right sides, the latter being proportional to distance
\[ u_\text{right} = u_\text{left} = u_\text{bot} + (u_\text{top} - u_\text{bot})\frac{y}{L} \]
where $L$ is the height of the plate.
A source of heat within the domain is assumed to exist. This source of heat may depend on time.
Parabolic PDE Solutions
A finite difference scheme is applied to solve the diffusion equation. First, discretization of the domain is performed. We use $h_x = h_y = h = 0.01$ and $dt = 0.0002$ to perform this discretization. We use a rectangular domain with $N_x \sim 2 N_y$ to accommodate the shape of the source. In what follows, $i$ and $j$ will specify $x$ and $y$ respectively, while the superscript $k$ will specify the time slice.
This diffusion PDE is written in finite difference form as
\[ \frac{u_{i,j}^{k+1} - u_{i,j}^k}{dt} = \frac{1}{2} \left( \partial_x^2 u_{i,j}^{k+1} + \partial_y^2 u_{i,j}^{k+1} + \partial_x^2 u_{i,j}^k + \partial_y^2 u_{i,j}^k + S_{i,j}^{k+1} + S_{i,j}^k \right) \]
where
\[ \begin{aligned} \partial_x^2 u_{i,j}^k = \frac{u_{i+1,j}^k - 2 u_{i,j}^k + u_{i-1,j}^k}{h^2} \\ \partial_y^2 u_{i,j}^k = \frac{u_{i,j+1}^k - 2 u_{i,j}^k + u_{i,j-1}^k}{h^2} \end{aligned} \]
and the same for the time slice specified by $k+1$. Separate the terms time $k+1$ from those at time $k$ to yield
\[ (1+4 r) u_{i,j}^{k+1} - r( u_{i+1,j}^{k+1} + u_{i-1,j}^{k+1} + u_{i,j+1}^{k+1} + u_{i,j-1}^{k+1}) = d_{i,j}^k \]
where $r \equiv \frac{dt}{2 h^2}$ and.
\[ d_{i,j}^k = r ( u_{i+1,j}^k + u_{i-1,j}^k + u_{i,j+1}^k + u_{i,j-1}^k ) + \frac{dt}{2}(S_{i,j}^{k+1} + S_{j,j}^k) \]
These equations form a linear system of equations for $u^{k+1}$. These are solved iteratively for each time slice. The lattice is swept repetitively, solving for $u_{i,j}^{k+1}$ until convergence is obtained.
\[ u_{i,j}^{k+1} = \frac{r d_{i,j}^k}{1 + 4 r} ( u_{i+1,j}^k + u_{i-1,j}^k + u_{i,j+1}^k + u_{i,j-1}^k ) \]
Serial Implementation
Download the tarball diffusion_2d.tar.gz containing a directory structure on any serial computer where you can use VisIt. Expand it with
tar -xzvf diffusion_2d.tar.gz
The directory diffusion_2d
contains two
directories: diffusion
and source
.
The code diffusion/diffusion_2d.c
implements the above
scheme to solve a 2D diffusion problem. The steps in the code are:
- Initialize the numerical grid and set the boundary values.
- Provide the initial field values.
- Set the source term.
- Loop over time.
- Iterate the solution until convergence at each time step.
- Store the field data every so often in time, whereby allowing to make a movie post-mortem.
Compile and execute the code with
gcc diffusion_2d.c -lm -o diffusion_2d ./diffusion_2d Nx Ny
where Nx
and Ny
are the (arbitrary)
number of grid points (image size). A ratio 2 to 1 is recommended for
the grid sizes in $x$ and $y$ directions. Good numbers to use are
(700, 400). The code defaults to scan over 3500 time steps.
The code diffusion_2d.c requires visit_writer.c and visit_writer.h with which it links. This package is used to store the field at each specified time slice in VTK file format in a sub-directory movie. VisIt can form a movie out of these files labeled by a frame number.
Notes:
- The iterative solution is obtained via double sweeps of the lattice for easier numerical solution.
- The code is written in the simplest possible way. The main program is modular to show the steps in the algorithm.
- The code outputs the resulting temperature field in files
labeled by a frame number. These files are written in a
sub-directory. The code
visit_writer.c
writes the data in VTK format. - The image data is written every so many (10) iterations through the time evolution of the temperature field, thereby allowing the production of a movie post-mortem.
- Setting
SOURCE
to0
solves the diffusion equation with no source. - When a source is included, the code displays the source term.
The code diffusion_2d.c
needs the
code source_2d.c
found in the
sub-directory source
. This code generates the source
term to include in the equations. The
sub-directory source
also
contains source_2d_main.c
to view the source. Compile and
execute with
gcc source_2d_main.c -lm -o source_2d ./source_2d
The shape of the source is chosen to represent the word "DREXEL".
This word is discretized in the function drexel_name()
to
occupy a good fraction of the domain. The shape of each letter is 1
pixel large. This shape is then enlarged via a smoothing procedure
(smooth_text()
). This is done via a simple sum of the
neighboring points adjacent to each pixel in the interior domain.
This produces the following source geometry (in a 200 x 100 grid
case):
Smoothing passes are repeated to obtain the desired widths for the
source (4 times in the code). The function see_image()
displays the source geometry (image above).
A time dependence is introduced in the code. A graphical representation (through Maple) follows: