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)

(**D** is the diffusion coefficient) is such that we ask for what is the
value of the *field* (wave)
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 1-D, 2-D, or 3-D. 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.

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,
is obtained from ,
which can not be parallelized.

We will model heat diffusion through a 2-D plate. The parabolic PDE to solve is

∂
u(x,y,t) / ∂ t = ∂^{2} u(x,y,t) / ∂x^{2} +
∂^{2} u(x,yt) / ∂y^{2} + s(x,y,t)

Dirichlet boundary conditions are assumed, the temperature being fixed at the
top and bottom of the plate, u_{top} and u_{bot} , and on the
left and right sides, the latter being proportional to distance

u_{right} = u_{left}
= u_{bot }+ delta_{t}
y
where delta_{t} = ( u_{top }- u_{bot }) /
L and L=height of plate

A source of heat within the domain is assumed to exist. This source of heat may depend on time.

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} ~ 2 N_{y} to acomodate 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

(u
^{k+1} _{ij} - u^{k}_{ij} ) / dt =

1/2 ( δ^{2}_{x} u ^{k+1} _{ij} +
δ^{2}_{y} u ^{k+1} _{ij} + δ^{2}_{x}
u ^{k} _{ij} + δ^{2}_{y} u ^{k} _{ij}
) + 1/2 ( S ^{k+1}_{ij} + S^{k}_{ij}
)

where

δ^{2}_{x}
u ^{k+1} _{ij} =
( u ^{k+1} _{i+1
j} - 2 u ^{k+1} _{ij} + u ^{k+1} _{i-1 j}
) / h^{2}

and
the same for the time slice specified by k. Let **r = dt / ( 2 h ^{2} )**.
The terms at time k+1 are separated from those at time k to yield

(
1+4r ) u ^{k+1} _{ij} -

r ( u ^{k+1}
_{i+1 j} + u ^{k+1}
_{i-1 j} + u ^{k+1} _{i j+1} + u ^{k+1} _{i
j-1}) = d ^{k} _{ij}

where

d ^{k}
_{ij} = r ( u ^{k} _{i+1
j} + u ^{k} _{i-1 j}
+ u ^{k} _{i j+1} + u ^{k} _{i j-1}) +

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

u ^{k+1}
_{ij} = ( r
( u ^{k+1} _{i+1 j}
+ u ^{k+1} _{i-1 j} +
u ^{k+1} _{i j+1} + u ^{k+1} _{i j-1}) +

Download the *tar* file Diffusion_2d.tar
containing a directory structure on any serial computer
where you can use *VisIt*.
Expand it via:

*tar -xvf Diffusion_2d.tar*

The directory *Difusion_2d* contains two directories: *diffusion* and *source*.

The code *2d_diffusion.c* in the directory *diffusion*
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

The code is compiled and executed via

gcc 2d_diffusion.c -lm -o 2d_diffusion

./2d_diffusion N_x N_y

where N_x and N_y 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 *2d_diffusion.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 about the code 2d_diffusion.c:*

- 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 to 0 solves the diffusion equation with no source
- When a source is included, the code displays the source term

*Better code:*

The code Diffusion_2d_pipe_python.tar makes the movie via our python utilities. This makes the movie in real time!

*The source:*

The code *2d_diffusion.c* needs the code
*2d_source.c* found in the
sub-directory *source*. This code generates the source term to
include in the equations. The sub-directory *source* also contains
*2d_source_main.c* to see the source. It
is compiled and executed via

gcc 2d_source_main.c -lm -o 2d_source

./2d_source

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 (function 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: