10. Example Workloads

In this step, we explore two parallel workloads: One harnessing MPI for the parallel computation of Pi using the Monte Carlo method. The other one uses MPI for the parallelization of the stencil-based finite difference method applied to the heat equation.

10.1. Monte Carlo Method

In this section we’ll use the Monte Carlo method to estimate \(\pi\). For this we randomly generate points inside the unit square \([0,1)^2\). Now, a point \((x,y)\) lies within the unit circle if \(x^2 + y^2 < 1\).

Fig. 10.1.1 Illustration of the Monte Carlo methods for the computation of \(\pi\). Shown is the unit square and the respective quadrant of the unit circle. Fifteen randomly sampled points (crosses) lie within the circle. Six samples (circles) lie outside of the unit circle.

This situation is illustrated in Fig. 10.1.1. If we randomly generate some points inside the unit square, we can simply derive how many of those lie inside of the unit circle. If \(N\) is the total number of generated points and \(C\) those inside of the circle, we obtain:

\[\pi \approx \frac{4C}{N}\]

In Fig. 10.1.1’s example we have \(N=21\) and \(C=15\), and thus would obtain \(\pi \approx 2.86\). This is still a pretty bad approximation, but we’ll boost accuracy by simply generating many more random samples.

Task

1. Use the std::mt19937_64 random number generator and std::uniform_real_distribution to generate a sequence of random numbers in \([0,1)\).

2. Parallelize the random number generation using MPI. Ensure that each process has its own independent instance of the random number generator.

3. Implement the described Monte Carlo method for the computation of \(\pi\). Calculate the error and required time-to-solution for 50,000,000 random samples.

4. Parallelize your implementation using your desired MPI communication.

10.2. 2D Heat Equation

A stencil is a computational pattern or algorithm that involves updating elements in a grid or array based on the values of neighboring elements. This pattern is commonly used in numerical simulations, such as those found in scientific and engineering applications like fluid dynamics, heat conduction, and image processing.

The heat diffusion equation is a partial differential equation which describes how heat or temperature spreads through a material over time. The Finite Difference Method (FDM) is a numerical technique used for solving partial differential equations, including heat transfer problems. This section only covers the bare minimum of the theory. For more details you might, for example, read chapter Wärmeleitung of the book Modellbildung und Simulation.

The basic idea is to discretize the spatial and temporal domains, and approximate the partial derivatives with finite differences. The two dimensional heat conduction equation is given by:

\[\frac{\partial T}{\partial t} = \alpha \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right),\]

where:

• \(T\) is the temperature.

• \(t\) is time.

• \(x\) and \(y\) are spatial coordinates.

• \(\alpha\) is the thermal diffusivity.

The partial derivatives (\(\frac{\partial}{\partial t}\), \(\frac{\partial^2}{\partial x^2}\), \(\frac{\partial^2}{\partial y^2}\)) are mathematical expressions representing the rate of \(T\) change with respect to time and spatial coordinates.

We discretize the partial differential equation for heat conduction using the Forward Time, Centered Space (FTCS) scheme. The method uses a five-point stencil as shown in Fig. 10.2.1.

Fig. 10.2.1 Illustration of the two-dimensional five-point stencil where the central point is surrounded by four neighboring grid points.

Let’s break down the discretization:

1. Discretize the Spatial and Temporal Domains:

   Discretize time: \(t^n=n \cdot \Delta t\) (where \(n\) is the ID of the discrete time step and \(\Delta t\) is the length of the time step).

   Discretize space: \(x_i=i \cdot \Delta x\) and \(y_j=j \cdot \Delta y\) (where \(i\) and \(j\) are grid IDs, and \(\Delta x\) and \(\Delta y\) are the distance of two grid points in \(x\) and \(y\) directions).

   Degrees of freedom: \(T_{i,j}^n\) is the discretized temperature at time step \(t^n\) for the grid point \((i,j)\).

2. Approximate the Partial Derivatives:

   Approximate the first derivative in time using a forward difference:

   \[\frac{\partial T}{\partial t} \approx \frac{T_{i,j}^{n+1} - T_{i,j}^n}{\Delta t}\]

   Approximate the second derivatives in space using centered differences:

   \[\frac{\partial^2 T}{\partial x^2} \approx \frac{T_{i+1,j}^n - 2T_{i,j}^n + T_{i-1,j}^n}{\Delta x^2}\]
   \[\frac{\partial^2 T}{\partial y^2} \approx \frac{T_{i,j+1}^n - 2T_{i,j}^n + T_{i,j-1}^n}{\Delta y^2}\]

3. Substitute into the original PDE:

\[\frac{T_{i,j}^{n+1} - T_{i,j}^n}{\Delta t} = \alpha \left( \frac{T_{i+1,j}^n - 2T_{i,j}^n + T_{i-1,j}^n}{\Delta x^2} + \frac{T_{i,j+1}^n - 2T_{i,j}^n + T_{i,j-1}^n}{\Delta y^2} \right)\]

4. Rearrange the equation to solve for the temperature at the next time step \(T_{i,j}^{n+1}\):

\[T_{i,j}^{n+1} = T_{i,j}^n + \frac{\alpha \Delta t}{\Delta x^2} \left( T_{i+1,j}^n - 2T_{i,j}^n + T_{i-1,j}^n \right) + \frac{\alpha \Delta t}{\Delta y^2} \left( T_{i,j+1}^n - 2T_{i,j}^n + T_{i,j-1}^n \right)\]

Here,

• \(T_{i,j}^n\) represents the temperature at the current grid point \((i, j)\) and current time level \(n\).

• \(T_{i,j}^{n+1}\) represents the temperature at the current grid point \((i, j)\) at the next time level \(n+1\).

• \(T_{i+1,j}^n\), \(T_{i-1,j}^n\), \(T_{i,j+1}^n\), and \(T_{i,j-1}^n\) are the temperatures at neighboring grid points of \((i, j)\) at the current time level \(n\) using 5-point-stencil:

  • \(T_{i+1,j}^n\) represents the temperature at the grid point to the right.

  • \(T_{i-1,j}^n\) represents the temperature at the grid point to the left.

  • \(T_{i,j+1}^n\) represents the temperature at the grid point above.

  • \(T_{i,j-1}^n\) represents the temperature at the grid point below.

5. Choosing \(\Delta t\):

In heat transfer simulations, the choice of the time step (\(\Delta t\)) is crucial for stability and accuracy of the simulation. We use the following for the length of \(\Delta t\) of our time steps:

\[\Delta t = \frac{\Delta x^2 \cdot \Delta y^2}{2 \cdot \alpha \cdot (\Delta x^2 + \Delta y^2)}\]

You are give a partly implemented heat transfer simulation in C++ which you can download from FSU-Cloud. The implemented parts of the parallel heat transfer simulation handle the creation of the transfer field and its distribution among the processors. You don’t need to worry about file read and write operations. The code is made to regularly save the simulation’s progress as picture files (PNG). You can adjust how often it does this and how often it writes checkpoints.

This simulation runs for 1000 time steps on a 2000x4000 grid with 0.01 spacing units in both the x and y directions. The material’s temperature evolves, guided by a diffusion constant of 0.5. The temperature inside and outside the disc is respectively set to 65.0°, and 5.0°. Snapshots are captured every 100 steps.

Task

1. Familiarize yourself with the functions and structs outlined in the header file heat.h of the downloaded archive.

2. Complete the implementation heat_update.cpp by coding the following four routines:

• start_halo_exchange

• complete_halo_exchange: If needed, this depends on the MPI communication you are using.

• update_interior_temperature

• update_boundary_temperature

3. Compute the value of \(\Delta t\) represented by the timeStep variable in the main.cpp file.

4. Utilize the functions you’ve implemented in the main.cpp file in the // ToDo section by arranging them in the correct order. In this section, the communication procedure you developed will be repeated for a specific number of iterations.

5. The parameters such as Simulation Steps, Grid Size, \(\Delta x\), \(\Delta y\), \(\alpha\) and Initialized temperatures can be found in the constants.h file. Play around with the parameters to see different results. Feel free to customize other files as well to suit your needs.

6. Submit all files you’ve worked on, changed or generated during execution. Write a one page report explaining exactly what you did and what parallel performance you observed. Share the results of your parameter study.

Hint

If you want to see your output as an infinite gif file you can use the following command:

convert -delay <delay_time> -loop 0 <file_names>.png <output_file_name>.gif

for example:

convert -delay 10 -loop 0 heat_*.png output.gif