Calculating Neutron Flux Using Parallel Programming in C with OpenMP and Monte Carlo Neutron Transport

The Monte Carlo Neutron Transport problem involves simulating the behavior of neutrons as they interact and move through a material or a system. Neutrons are subatomic particles that play a crucial role in nuclear reactions and various applications, such as nuclear energy, nuclear weapons, and medical imaging. In this problem, the objective is to model the movement of neutrons in a specific environment and study their interactions with the materials present. Neutrons can scatter, be absorbed, or undergo fission as they travel through the medium, making their behavior complex and challenging to predict analytically.

Overall, the Monte Carlo Neutron Transport problem allows us to understand and optimize the behavior of neutrons in diverse scenarios, aiding in the design and safety assessment of nuclear systems and technologies.

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

#define NUM_PARTICLES 1000000
#define NUM_THREADS_LIST {1, 2, 4, 8, 16}

// Function to perform neutron transport simulation
double monteCarloNeutronTransport(int num_particles) {
    int i, num_collisions = 0;
    double neutron_flux = 0.0;

    // Set seed for random number generator
    unsigned int seed = (unsigned)time(NULL) + omp_get_thread_num();

    // Monte Carlo simulation loop
    for (i = 0; i < num_particles; i++) {
        double position = 0.0; // Neutron's initial position

        // Particle movement loop
        while (position >= 0.0 && position <= 1.0) {
            double random_num = (double)rand_r(&seed) / RAND_MAX; // Generate a random number between 0 and 1

            // Probability of collision = 0.5
            if (random_num < 0.5) {
                num_collisions++;
                break;
            }

            // Move neutron randomly to the left or right
            if (random_num < 0.75) {
                position -= 0.1; // Move left
            } else {
                position += 0.1; // Move right
            }
        }
    }

    // Calculate neutron flux
    neutron_flux = (double)num_collisions / num_particles;
    return neutron_flux;
}

int main() {
    double neutron_flux_total = 0.0;
    int i;

    // Define an array of thread counts to test
    int num_threads_list[] = NUM_THREADS_LIST;
    int num_threads_list_size = sizeof(num_threads_list) / sizeof(num_threads_list[0]);

    // Perform the simulation with a single thread for serial execution
    omp_set_num_threads(1);
    double start_time_serial = omp_get_wtime();
    double neutron_flux_serial = monteCarloNeutronTransport(NUM_PARTICLES);
    double end_time_serial = omp_get_wtime();
    double execution_time_serial = end_time_serial - start_time_serial;

    printf("Serial Execution Time: %lf seconds\n", execution_time_serial);

    double average_neutron_flux = 0.00;
    // Perform the simulation for different thread counts (excluding the first entry since it corresponds to serial execution)
    for (i = 1; i < num_threads_list_size; i++) {
        int num_threads = num_threads_list[i];

        // Set the number of threads for OpenMP
        omp_set_num_threads(num_threads);

        // Record the starting time
        double start_time = omp_get_wtime();

        // Perform the simulation using OpenMP parallel for loop
        #pragma omp parallel for reduction(+:neutron_flux_total)
        for (int j = 0; j < num_threads; j++) {
            double neutron_flux_thread = monteCarloNeutronTransport(NUM_PARTICLES / num_threads);
            neutron_flux_total += neutron_flux_thread;
        }

        // Calculate the average neutron flux over all threads
        average_neutron_flux = neutron_flux_total / num_threads;

        // Record the ending time
        double end_time = omp_get_wtime();

        // Calculate the execution time
        double execution_time = end_time - start_time;

        // Calculate the speedup relative to the serial execution
        double speedup = execution_time_serial / execution_time;
        
        printf("Threads: %d, Execution Time: %lf seconds, Speedup: %lf, Average Neutron Flux: %lf\n",
               num_threads, execution_time, speedup, average_neutron_flux);

    }

    return 0;
}

The purpose of the code is to run the neutron transport simulation in parallel and measure the processing time and acceleration for different thread counts. This illustrates an example of evaluating the performance of the simulation and understanding the benefits of parallel computing.

Resources:

OpenMP official website: https://www.openmp.org

Books

Monte Carlo Neutron Transport: Concepts, Codes, and Applications

This book provides a comprehensive overview of Monte Carlo methods for neutron transport simulations and their applications. It covers the theoretical foundations as well as practical implementation details.

Parallel Programming in C with OpenMP by Michael J. Quinn

This book focuses on parallel programming techniques using OpenMP in C. It covers various parallelization concepts and strategies to optimize code for multi-core architectures.