R.Anush R.Anush - 2 months ago 18
C Question

Difficulty in implementing the neutron transport model

I'm working on simulation of neutron transport model using monte carlo method. I'm first implementing the sequential algorithm that is given in the book parallel programming with openmp and MPI by M. J. Quinn

The below (sequential) code given is a simulation of neutron transport model. when i run the executable file for the first time, it is running in infinite loop. When i stopped and ran again, the output i'm getting is 0, 0, 0 and is not the expected output for all possible values of input parameters, given in the question as mentioned in solution manual of M.J.Quinn parallel programming with openmp and MPI (Q.No 10.8).

The soft copy for the book is given below:
https://docs.google.com/document/d/19aKs6XF3Gt6BRGSGVysZqam3Cb97rgRvTCYP2Nw-C5M/edit

The expected output is: the amount of reflected, absorbed and transmitted rates are: 25.31, 46.13, 28.56. No error was given when the code was compiled. please help me out. here is the following code that i wrote:

struct neutron
{
long double d ; // Direction of the neutron (measured in radians between 0 and pi)
long double x ; // Position of particle in plate(0<=x<=H)where H is the thickness of the plate
};
typedef struct neutron neutron;


int main()
{
srand(time(NULL));
double C ; // mean distance between neutron/atom interactions is 1/C
double C_s; // Scattering component of C
double C_c; // Absorbing component of C
int r=0, b=0, t=0; //Counts of reflected, absorbed, transmitted neutrons
int i=0;
int n; // Number of Samples

//double d; //Direction of the neutron (measured in radians between 0 and pi)
int a; //a==0 implies false 1 otherwise
double u; // Uniform random number
double L; // Distance neutron travels before collision
neutron nt;
double H;
double p,q,o;


printf("\n Enter the value of C_s:\n");
scanf("%lf",&C_s);
printf("\nEnter the value of C_c:\n");
scanf("%lf",&C_c);
//printf("\nEnter the value of L:\n");
//scanf("%lf",&L);
printf("\nEnter the value of H (Thickness of the plate):\n");
scanf("%lf",&H);
printf("\n Enter the number of samples you want to simulate for....\n");
scanf("%d",&n);
for(i=1;i<=n;i++)
{
u=rand()%n;
u=1/u;
nt.d=0;
nt.x=0;
a=1;
while(a)
{
L=-(1/C)*log(u);
nt.x=nt.x+L*cos(nt.d);
if (nt.x<0)
{
// printf("\nparticle %d got reflected\n",i);
r++;
a=0;
}
else
{
if(nt.x>=H)
{
// printf("\n particle %d got transmitted\n",i);
t++;
a=0;
}
else
{
if(u<C_c/C)
{
// printf("\n the particle %d got absorbed\n",i);
b++;
a=0;
}
else
nt.d=u*(22/7);
}
}
}
}
o=(r/n)*100;
p=(a/n)*100;
q=(t/n)*100;
printf("\nthe amount of reflected, absorbed and transmitted rates are: %lf, %lf, %lf\n", o, p, q);
return 0;
}

Answer

You have

o=(r/n)*100;

which performs an integer division before the multiplication. Since your expected output is 25.31 that would mean the result of r / n should be 0.2531, but the result of such an integer division is 0 which supports the 0 you are getting.

I suggest this

o = 100.0 * r / n
p = 100.0 * a / n;
q = 100.0 * t / n;    

This will perform the correct arithmetic.

In addition you are using uninitialised variables.