Numerical Simulation and Parallel Programming
Notebook
Version 1.02
Haide College, OUC
2023.08
1
0 Information
Author：Fukaya Takeshi, Iwashita Takeshi (Hokkaido University), ZHANG Linjie, FU Kai (Ocean University of China)
E-mails: zhanglinjie@ouc.edu.cn, kfu@ouc.edu.cn
✓Outcomes ✏
1. Learn numerical simulations for simple physical models described by differential equations.
2. Learn the basic of MPI (Message-Passing Interface) parallel programming.
✒ ✑
Keyword: Heat conduction equation, finite difference method, MPI parallelization
1 Schedule This part has 30 credit hours in total.
1-4: How to program on Linux server.
5-8: Learn numerical simulation for 1-D heat conduction.
9-12: Learn numerical simulation program for 2-D heat conduction.
13-16: Learn the basic of parallel programming with MPI.
17-24: Parallelize the 2-D heat conduction simulation program with MPI.
It is recommended that make use of the class time effectively. Keep up with the programming even if
you don’t understand something.
2 How to submit a report Explain at the beginning (or in the middle) of the class.
• Time Limit: Further notice.
• How to submit: Print in A4 paper. Single-sided printing, double-sided printing, black and white
printing, color printing, are acceptable.
3 Academic Integrity It is unacceptable to submit work that is not your own. Assignments will be
checked for similarity and copying to ensure that students work is their own. Information sources must
be paraphrased (put in your OWN words - not copied directly) and the reference must be included. You
are also responsible for complying with any additional rules related to assessments communicated to you
by your instructors. Assignments that are found to breach this may be given zero. Further actions may
be taken.
2
4 Precautions about this part
• Represent numerical values (for example, temperature value) in double decision.
• It will not lose points if your program have to be recompiled every time when the setting changes.
• In the sample program, “ ” means a half-width space.
5 Precautions about report
• When visualizing the computing results, use the actual coordinates or time, instead of the number
of steps (space or time).
• Don’t insert too many pictures or graphs unless it is necessary.
• All pictures and graphs should be mentioned in the report.
• If you cite a paper or a website, identify it as a reference.
• Formulas should be concise, and be explained in words.
3
1 How to program on Linux server
Important information of server
• Server IP address: Further notice
• User Account: Further notice
• Password: Further notice
Basic Linux Commands The server is a Linux system (a type of operating system) which is
considerably different from Windows, and similar to Mac. In order to be more productive in this class,
you need to learn some basic about how to get things done in the Linux system.
1. ls (short for list). To find out what is in your current directory, type
$ ls
To list all files in your directory including those whose names begin with a dot (hidden files), type
$ ls -a
2. mkdir (make directory). Make a subdirectory in your current working directory, for example test.
Type
$ mkdir test
3. cd (change directory)
To change to the directory you have just made, type
$ cd test
Exercise 1a Make another directory inside the test directory called backups
4. The directories . and ..
In Linux, “.” means the current directory, ”..” means the parent of the current directory, so typing
$ cd .. (NOTE: there is a space between cd and the dot)
4
will take you one directory up the hierarchy.
5. pwd (print working directory). To find out the absolute pathname of your working directory, type
$ pwd
6. ∼ (your home directory). To list the contents of your home directory, type:
$ ls ~
7. cp (copy). Use “echo” and “>” to creat a file. Change to directory “test”, and copy the file
NumSim.txt to the current directory, type:
$ echo "Numerical simulation" > NumSim.txt
$ cd ~/test
$ cp ~/NumSim.txt .
8. mv (move). mv file1 file2 moves (or renames) file1 to file2. mv file1 dir moves file1 to directory
dir. For example:
$ cp NumSim.txt NumSim.bak.txt
$ mkdir backups
$ mv NumSim.bak.txt backups
9. Removing files and directories rm (remove), rmdir (remove directory)
As an example, create a copy of the NumSim.txt file then delete it.
$ cp NumSim.txt temp.txt
$ ls (to check if it has created the file)
$ rm temp.txt
$ ls (to check if it has deleted the file)
Summary
ls list files and directories
ls -a list all files and directories
mkdir make a directory
cd directory change to named directory
cd change to home-directory
cd ∼ change to home-directory
cd .. change to parent directory
pwd display the path of the current directory
cp file1 file2 copy file1 and call it file2
mv file1 file2 move or rename file1 to file2
rm file remove a file
rmdir directory remove a directory
Table 1.1: Summary of basic Linux commands
5
Review of C programming While reviewing the C language programming, we will try to compute
C ← AB + C (1.1)
as an example, where A, B, C are n × n matrices. There are 6 steps:
1. Read matrix dimension n from command-line argument.
2. Allocate memory dynamically for A, B, C.
3. Initialize elements of matrix A, B, C.
4. Compute each element of matrix C with Eq. (1.1).
5. Display execution time and FLOPS.
6. Free memory, etc.
In step 1, dimension n of matrix is given from command-line arguments,
$./a.out 1000
where 1000 is the given value of n. In the program, we can use
n = atoi (argv[1]);
to read the command-line argument 1000 and save it to variable n.
In step 2, it is straightforward to treat matrix as 2-D array in C programming. But, it is a little
complicated to create a 2-D array dynamically. So we will treat matrix as a 1-D array with size n2. Let
ai,j be the element belong to column i and row j in matrix A, then the corresponding element in 1-D
array is A[j ∗ n + i] (0 ≤ i, j ≤ n − 1). we use column-major form here. After declaring a double-type
pointer ∗A, we allocate memory space for A with the following code.
A = (double *)malloc (sizeof (double)*n*n);
Note that if we don’t know the size of array until the program is executed, we have to allocate memory
space dynamically.
In step 3, in order to make it easier to verify the correctness of our program, the initial values of matrix
A, B, C are set as
aij = 1
√n
, bij = 1
√n
, cij = 0. (1.2)
(It’s not difficult to find the the values of matrix C after Eq. (1.1).)
Now consider step 4, which is the main part of this program. The calculation of Eq. (1.1) is
cij ← cij +
n
!−1
k=0
aikbkj , (for each i, j). (1.3)
So, we can write the program like
for (i = 0; i < n; i++)
{
for (j = 0; j < n; j++)
{
/* Compute each element of matrix $C$ */
}
}
6
We can confirm whether our program runs correctly by checking the values of matrix C after calculation.
For measuring the execution time, we can write a function get time as following.
#include
double get_time ()
{
struct timeval tv;
gettimeofday (&tv, NULL);
return tv.tv_sec + (double)tv.tv_usec*1e-6;
}
Then we can get the elapsed time (seconds) like
t0 = get_time ();
/* process to be measured */
t1 = get_time ();
time = t1 - t0;
In addition to the execution time, we can calculate FLOPS (Floating-point Operations Per Second) too.
FLOPS = Times of Floating-point Operations
Elapsed Time (second) (1.4)
Since the number of floating-point operations in Eq. (1.1) is 2n3, so flops in our program can be computed
as
flops = 2.0* (double)n* (double)n* (double)n/time;
Then we display the execution time and FLOPS.
printf ("n!=!%d,!!time!(sec)!=!%8.3e,!!FLOPS!=!%8.3e\n", n, time, flops);
Finally, free memory allocated in step 2.
free (a);
7
Task 1.1
Most compilers have optimization option. Depending on this option, our program is optimized
at compile time. Please check how the program performance influenced by optimization option.
Write a program that performs the calculation in Eq. (1.1). Then compile it like
$gcc main.c -OX -lm
The letter after the hyphen − is the uppercase letter O, and the X part is a number from
0(zero) to 3. The higher the number, the stronger the compiler optimization (0 means no optimization). ”-lm” may be necessary if you include math.h in your program.
Set n = 1000, 2000, 3000, change the optimization option from -O0 to -O3 and measure the
execution time and FLOPS for each case.
Task 1.2
The main part of this program (the part that calculates Eq.(1.1)) is a triple loop structure, and
the loop variables are i, j, k. The order of these loops can be changed. For example, j → k →
i does not change the final result, but the memory access pattern (the way of reading/writing
the elements of array) will be different. Set n = 1000, 2000, 3000, change the loop order (6 types
in total), and check the performance. Note that, in order to eliminate the effects of compiler
optimization, please compile your program with optimization option -O0. If the execution time
(performance) differs, please explain it as far as you can understand.
8
2 Numerical simulation for 1-D heat
conduction
Computer simulations are widely used in the fields of engineering and natural sciences. Generally, when
simulating a certain physical phenomenon, four steps are included.
1. Modeling: derive one or more partial differential equations (PDEs) to describe the phenomenon.
2. Discretization: transform partial differential equations into discrete form.
3. Calculation: solve discretized problems with computer, get the numerical solutions of partial differential equations.
4. Visualization: Display calculation results into visible form.
We will experience the above flow through a simple physical phenomenon as an example.
Problem setting: One-dimensional heat conduction Consider the conduction of heat
along a stick of length L as shown in Fig. 2.1. We don’t care about the cross-sectional area since it’s very
small relative to the length. Let the left endpoint be the original point (x = 0), and the x-axis is taken in
the length direction. Heat conduction happens along ±x direction. In addition, heat exchange between
this stick and the outside world occurs only at the end points.
Fig 2.1: A stick
Derivation of partial differential equations Let u(x, t)denote the temperature at position x
at time t, 0 ≤ x ≤ L, 0 ≤ t ≤ T. 1 This 1-D heat conduction phenomenon can be described by a partial
differential equation. (See Task 2.3 for the derivation.)
∂
∂t
u(x, t) = α
∂2
∂x2 u(x, t). (2.1)
1Consider heat conduction up to time t = T.
9
α is temperature diffusivity. Eq. (2.1) is solvable when initial condition u(x, 0) (temperature distribution
of the object at time t = 0) and boundary condition (condition related to temperature at both left and
right ends) are given.
Discretization Most of the PDEs that appear in engineering and the natural sciences are difficult to
solve analytically. Therefore, instead of finding an analytical solution, computer is used to get numerical
solution.
In general, the domains in which the PDE is defined are continuous, and the solution is continuous
too. In case of Eq.(2.1), they are [0, L] × [0, T] and u(x, t) separately. So, they cannot be handled directly
on the computer, therefore “discretization” is required.
Now we will study how to discretize the 1-D heat conduction problem given in Eq.(2.1). As mentioned
before, it is not possible to handle all points (x, t) of 0 ≤ x ≤ L, 0 ≤ t ≤ T on a computer. Therefore, we
separate the spatial domain x by step size ∆x and the time domain t by step size ∆t. (See Fig. 2.2 and
Fig. 2.3)
xi := i∆x (i = 0, 1,..., L
∆x =: Mx), (2.2)
tn := n∆t (n = 0, 1,..., T
∆t =: N). (2.3)
Then we only compute the solutions on the discretized points (xi, tn) only, that is
un
i := u(xi, tn). (2.4)
Fig 2.2: Spatial discretization Fig 2.3: Temporal discretization
Now, we will compute the solutions on (Mx + 1) × (N + 1) discrete points. We approximate the
derivative with forward difference in time and with central difference in space. Then we get
un+1
i − un
i
∆t = α
un
i+1 − 2un
i + un
i−1
(∆x)2 . (2.5)
See Task 2.4 for details. Note that Eq. (2.5) is defined only for the points i = 1, 2,...,Mx − 1, end points
are not included. For end points (i = 0, Mx), we have
un
0 := u(0, tn), un
Mx := u(L, tn) (2.6)
by applying the given boundary conditions.
Calculation of numerical solution Discretization has made it possible to handle the PDE
problem on a computer.
First, Transform Eq. (2.5) to the following from
un+1
i = ··· . (2.7)
10
Now we can see that only the value at time step n appears on the right side. In other words, if un
i is
known, the values of un+1
i can be easily obtained from Eq. (2.7). 2 The values at time step n = 0 is given
by initial condition.
u0
i := u(xi, 0). (2.8)
Therefore, we take the values at time step n = 0 as a starting point, then we compute the values of other
time steps n = 1, 2,... with Eq. (2.7) in order.
See Program 2.1 for the skeleton. Note that, in this program, instead of setting the step sizes ∆x and
∆t, we use the number of discrete points in each dimension Mx, N, which can be specified from run-time
arguments. print data (...) is a function (described later) that outputs the calculation results to file.
Program 2.1: Overview of 1-D heat conduction simulation program by finite difference method
1 /∗ Include required header files ∗/
2 #define L 1.0 //Length of bar
3 #define T 1.0 //Maximun time to simulation
4 #define alpha 1.0 //Thermal diffusivity
5 #define T SPAN 20 //Result output frequency (adjusted appropriately)
6
7 int main(int argc, char ∗argv[])
8 {
9 /∗ Declaration of variables (add each necessary variable) ∗/
10 double ∗u, ∗uu, dx, dt;
11 int Mx, N;
12
13 /∗ Get the values of Mx and N from the arguments ∗/
14
15 /∗ Dynamic allocation of memory, variable setting, etc. ∗/
16
17 for(i = 0; i <= MX; i++)
18 {
19 /∗ Set initial conditions for u ∗/
20 }
21 print data(...); //Output the initial state to a file
22
23 for(n = 1; n <= N; n++)
24 {
25 for(i = 1; i < MX; i++)
26 {
27 /∗ Calculate the value (uu) of the next time step from the current value (u) ∗/
28 }
29 for(i = 1; i < MX; i++)
30 {
31 /∗ Copy the value of uu to u ∗/
32 }
33 if(n%T SPAN == 0)
34 {
35 print data(...); //Output to a file at a regular frequency
36 }
37 }
38
39 /∗ Free memory ∗/
40
41 return 0;
42 }
2At a certain point i, values at time step n + 1 can be (easily) obtained from the values at the previous step n. We call
this explicit method.
11
Visualization When analyzing the calculation results, it is important to show them in a visually easyto-see form, that is visualization. In this experiment, visualization will be performed using a free software
gnuplot, which can graph data stored in file.
First, save the calculation results (temperature in this sample) of a certain time step to file (See
Program 2.2). Function print data creates a file named data ******.dat, and outputs the x-coordinate
and computed temperature (separated by space) to this file. 3
Program 2.2: Function that saves the calculation result to file
1 void print data(int Mx, int n, double dx, double dt, double ∗u)
2 {
3 FILE ∗fp;
4 char sfile[256];
5 int i;
6
7 sprintf(sfile, "data_%06d.dat", n);
8 fp = fopen(sfile, "w");
9 fprintf(fp, "#time!=!%lf\n", (double)n ∗ dt);
10 fprintf(fp, "#x!u\n");
11 for(i = 0; i <= MX; i++)
12 {
13 fprintf(fp, "%lf!%12lf\n", (double)i∗dx, u[i]);
14 }
15 fclose(fp);
16 return;
17 }
Now we can use gnuplot to display data in a single file. First, create a script files shown in Programs
2.3 in the same folder with your program. Then run the following command from the terminal.
$gnuplot sample0.gp
Program 2.3: Scipt of creating a gif still image (sample0.gp)
1 set terminal gif
2 set output "test_still.gif" #Adjust as needed
3 set xrange [0:1] #Adjust as needed
4 set yrange [0:100] #Adjust as needed
5 plot "***.dat" using 1:2 with line
6 unset output
If it works, a gif animation image file named test still.gif will be found in the same folder.
Since the temperature may change over time, it will be a good idea to create a gif animation to visualize
the state of time change. Here’s how it works. First, create two script files shown in Programs 2.4 and
2.5 in one folder. Then run the following command from the terminal.
$gnuplot sample1.gp
If it works, a gif animation image file named test.gif will be found in the same folder.
Program 2.4: Script of creating a gif animation (sample1.gp)
1 reset
2 set nokey
3 set xrange [0:1] #Adjust as needed
4 set yrange [0:100] #Adjust as needed
5 set size square
6
3Lines beginning with # are ignored as comments when plotting with gnuplot.
12
7 set terminal gif animate delay 50 #Display interval can be adjusted by changing 50
8 set output "test.gif" #Specifying the name of the output file
9
10 n0 = 0
11 nmax = 100 #Number of divisions in the time direction
12 dn = 20 #Step size (program T SPAN value)
13
14 load "sample1.plt"
15
16 unset output
Program 2.5: Script of creating a gif animation (sample1.plt)
1 if(exist("n")==0 || n<0) n = n0
2
3 title n = sprintf("t!=!%f", 1.0∗n/nmax)
4 unset label
5 set label 1 at graph 0.05,0.95 title n #The position of the legend can be adjusted by changing (0.05, 0.95)
6
7 filename = sprintf("data_%06d.dat", n)
8 plot filename using 1:2 with line
9
10 n = n + dn
11 if(n <= nmax) reread
12
13 undefine n
Task 2.1
Set L = 1.0, T = 1.0, α = 1.0.
• Initial conditions: u(0, 0) = 20.0, u(x, 0) = 0.0 (0