[MPI / SPRNG / OpenMP] A neutron transport Monte Carlo class project

neutrony.cpp

#include <omp.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <mpi.h>
#include <string.h>

#define SPRNG_SIMPLE 1
#define SEED 985456376
#define SPRNG_GENERATOR_TYPE DEFAULT_RNG_TYPE
// Define random number function

#include "sprng.h"
#include "sprng_cpp.h"

#define M_PI 3.14159265358979323846

//computation parameters
int H;
int ITERS;
double CC;
double CS;

// computed after
double C;
double NINV_C, P_ABSOR, P_BOUNC;

int base_stream_id = 0;	// identifier of zeroth-threads SPRNG stream
int total_streams = 0;	// total amount of streams
Sprng ** random_streams;	// array of streams
/**
* Simulates running a single neutron. 
*
* Returns below constants on given event
*/

#define TRANSMITTED 0
#define ABSORBED 1
#define REFLECTED 2
int simulate_single_neutron(int h, Sprng * local_stream) {
    double x = 0;       // Position in the plate
    double dir = 0;     // Current direction of the neutron
    while (1) {
        x += NINV_C * log(local_stream->sprng()) * cos(dir);
        if (x < 0) return REFLECTED;              // Reflected, out of plate!
        if (x >= h) return TRANSMITTED;           // Transmitted, out of the plate!
        if (local_stream->sprng() < P_ABSOR) return ABSORBED;     // Absorbed!
        else dir = local_stream->sprng() * M_PI;                  // Else is used here because P_ABSOR+P_BOUNC=1
    }
    abort();            // Should never get here!
}

/**
* Simulates a particular number of neutrons
*
* Stores result in passed 3-element array. Caller has to take care
* of cleaning the array.
*/
#pragma omp private(local_stream_id) reduce(+: locresults)
void simulate_neutrons(long long rounds, long long results[][3], int h) {
	Sprng * local_stream = random_streams[omp_get_thread_num()];

	int locresults[] = {0, 0, 0};

    while (rounds--)
		locresults[simulate_single_neutron(h, local_stream)]++;

#pragma omp critical
	{
		for (int i=0; i<3; i++) results[h-1][i] = locresults[i];
	}
}

int main(int argc, char * argv[]) {
    int n, myid, h;    // MPI bookkeeping
    
    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &n);
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
	if (myid==0) printf("Started. Got %d nodes.\n", n);

	// Get amount of threads for OpenMP - every thread
	int * threads_requested_c = new int[n];	// input for reduce
	int * threads_requested = new int[n];	// output for reduce
    
    for (int i=0; i<n; i++)
	   threads_requested_c[i] = (i == myid) ? omp_get_max_threads() : 0;
    
	MPI_Allreduce(threads_requested_c, threads_requested, n, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

	for (int i=0; i<n; i++) total_streams += threads_requested[i];		// calculate total stream
	
	if (myid == 0)	printf("Got %d threads total.\n", total_streams);

	for (int i=0; i<myid; i++) base_stream_id += threads_requested[i];	// which stream do I start from?

	// Allocate streams
	random_streams = new Sprng*[threads_requested[myid]];
	for (int i=0; i<threads_requested[myid]; i++) {
        random_streams[i] = SelectType(SPRNG_GENERATOR_TYPE);
        random_streams[i]->init_rng(base_stream_id+i, total_streams, SEED, 0);
    }

    // Read parameters from stdin. Broadcast them. Compute dependencies
    if (myid == 0) {
        printf("Enter H, total iterations, CC and CS: ");
        scanf("%i %i %lf %lf", &H, &ITERS, &CC, &CS);
        printf("H=%i, iterations=%i, CC=%lf, CS=%lf\nComputation started...\n", H, ITERS, CC, CS);
        H++;
    }
    
    MPI_Bcast(&H, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast(&ITERS, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast(&CC, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(&CS, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    
    C = CC + CS;
    NINV_C = -1/C, P_ABSOR = CC / C, P_BOUNC = CS / C;
    
    // Run the simulation
	long long results[H][3];
    for (int h=1; h<H; h++) memset(results[h-1], 0, sizeof(long long)*3);
        
#pragma omp parallel for
    for (int h=1; h<H; h++) simulate_neutrons((long long)ceil(ITERS/n), results, h);
    
    // Reduce the calculations
    long long final_results[H][3];
    MPI_Reduce(&results, &final_results, 3*H, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
    
    // Output the information
    if (!myid) {
        printf("------------\nSimulation results:\n\n");
        
        for (int h=1; h<H; h++) {
            double total_iters_ran = final_results[h-1][ABSORBED] + final_results[h-1][TRANSMITTED] + final_results[h-1][REFLECTED];
            printf("---\tH = %d\n", h);
            printf("   Absorption ratio:  %lf\t", final_results[h-1][ABSORBED]/total_iters_ran);
            printf("   Transmission ratio:  %lf\t", final_results[h-1][TRANSMITTED]/total_iters_ran);
            printf("   Reflection ratio:  %lf\n", final_results[h-1][REFLECTED]/total_iters_ran);
        }
    }
    
    MPI_Finalize();

    for (int i=0; i<threads_requested[myid]; i++)
        delete random_streams[i];
    delete random_streams;
    delete threads_requested_c;
    delete threads_requested;

    return 0;
}