2016-09-07 3 views
3

モンテカルロ法を使った中性子輸送モデルのシミュレーションに取り組んでいます。私は最初に書籍で与えられたシーケンシャルアルゴリズムを実装しています。OpenMPとMPIによる並列プログラミングMによってM. Quinn中性子輸送モデルの実装が難しい

与えられた以下の(連続した)コードは、中性子輸送モデルのシミュレーションです。最初に実行ファイルを実行すると、無限ループで実行されています。私が停止して再度実行したとき、出力は0,0、0であり、入力パラメータのすべての可能な値に対して期待される出力ではなく、MJQuinnの解説書でopenmp MPI(Q.No 10.8)。

本のためのソフトコピーは以下の通りである: https://docs.google.com/document/d/19aKs6XF3Gt6BRGSGVysZqam3Cb97rgRvTCYP2Nw-C5M/edit

予想される出力は次のとおりです。反映の量、吸収され、送信料金は以下のとおりです。25.31、46.13、28.56。コードのコンパイル時にエラーは発生しませんでした。私を助けてください。あなたのコードは、UB(未定義の動作)であるので、

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; 
} 
+0

1)_value_は、 'double C 'の' C'とは何ですか? ... L = - (1/C)* log(u); '? 2) 'cos(nt.d)'の代わりに 'cosl(nt.d)'を使うことを提案する。 3)コンパイラの警告を有効にする - あなたの時間を節約します。 – chux

+0

「C」を初期化することはできません。 – 4386427

+1

トピックの:常に 'scanfによって返された値を確認してください – 4386427

答えて

2

あなたは乗算の前に整数の除算を行い

o=(r/n)*100; 

を持っています。あなたの予想される出力は25.31なので、の結果は0.2531でなければなりませんが、そのような整数除算の結果は0をサポートする0です。

私はこの

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

これは正しい算術演算を実行します示唆しています。

さらに、未初期化変数を使用しています。

+0

良いキャッチ..... – 4386427

2

あなたは少なくとも一つの初期化されていない変数を持っている:ここで私が書いた次のコードがあります。

L=-(1/C)*log(u); 
    ^
     uninitialized 

さらに、これは奇妙に見える:

while(a) 
{ 
    .... 
} 

whileを出るときにaは常にゼロになります。しかし、まだあなたは:

p=(a/n)*100; 

pは常にゼロになります。それはおそらくあなたが望むものではありません。

はしかし、私はこの行を参照してください。

b++; 

を私はあなたがbを使用している任意のコードを見つけることができません。

本当にこのようにしたいですか?

p=(b/n)*100; // b instead of a 
0

私は私のクラスメートと議論することにより、自分のコードでミスを発見し、それを整流し、ここにも、コードを並列化次のコードです:

/** 
    To compile the program:   gcc sequential.c -lm 
    To run the program:   ./a.out 
**/ 
#include<stdio.h> 
#include<math.h> 
#include<stdlib.h> 
#include<time.h> 

#define TOTAL_NUMBER_OF_NEUTRONS 10000000        //Total number of neutrons 
#define SCATTERING_COMPONENT 0.7         // Scattering component of C 
#define ABSORBING_COMPONENT 0.3         // Absorbing component of C 
#define THICKNESS_OF_THE_PLATE 8         //Thickness of the Plate 
#define NUMBER_OF_ITERATIONS 1000 
int no_of_particles_reflected=0, no_of_particles_absorbed=0, no_of_particles_transmitted=0; 
double total_amount_absorbed,total_amount_reflected,total_amount_transmitted; 
double distance_travelled_before_collision;          // Distance travelled by neutron before colliding from the plate 
double C ;            // mean distance between neutron/atom interactions is 1/C 
int flag= 1;                       int counter= 0; 
struct neutron 
{ 
    double d;           // Direction of the neutron (measured in radians between 0 and pi) 
    double x;           // Position of particle in plate(0<=x<=H),H is the thickness of the plate 
}; 
struct neutron nt; 
void neutron_model() 
{ 
    double random_number_generator;           // Uniform random number generator 
    total_amount_absorbed=0,total_amount_reflected=0,total_amount_transmitted=0; 
    for(int i=1;i<=TOTAL_NUMBER_OF_NEUTRONS;i++) 
    { 
     random_number_generator=rand()%TOTAL_NUMBER_OF_NEUTRONS; 
     random_number_generator=1/random_number_generator; 
     nt.d=0; 
     nt.x=0; 
     flag=1; 
     while(flag) 
       { 
      distance_travelled_before_collision=-(1/C)*log(random_number_generator); 
      nt.x=distance_travelled_before_collision*cos(nt.d); 

      if (nt.x<0) 
      { 
       no_of_particles_reflected++; 
       flag=0; 
      } 
      else if(nt.x>=THICKNESS_OF_THE_PLATE) 
      { 
       no_of_particles_transmitted++; 
       flag=0; 
      } 
      else 
      { 
       if(random_number_generator<(ABSORBING_COMPONENT/C)) 
       { 
        no_of_particles_absorbed++; 
        flag=0; 
       } 
       else 
       { 
        nt.d=random_number_generator*(22.0/7.0); 
        counter++; 
        if(counter==NUMBER_OF_ITERATIONS) 
        { 
         no_of_particles_absorbed++; 
         flag=0; 
         counter=0; 
        } 
       } 
      } 
     } 
    } 
} 
int main(int argc, char *argv[]) 
{ 
    C=SCATTERING_COMPONENT + ABSORBING_COMPONENT; 


    clock_t start, end; 
    double diff_t; 
    diff_t = 0.0; 
    //start the clock 
    start= clock(); 
    neutron_model(); 
    //stop the clock 
    end = clock(); 
    diff_t = ((double) (end - start))/CLOCKS_PER_SEC; 
    total_amount_reflected=((double)no_of_particles_reflected/(double)TOTAL_NUMBER_OF_NEUTRONS)*100.0; 
    total_amount_absorbed=((double)no_of_particles_absorbed/(double)TOTAL_NUMBER_OF_NEUTRONS)*100.0; 
    total_amount_transmitted=((double)no_of_particles_transmitted/(double)TOTAL_NUMBER_OF_NEUTRONS)*100.0; 

    printf("\nDATA VALUE ASSUMED:-\n"); 
    printf("***********************************************************************\n"); 
    printf("TOTAL_NUMBER_OF_NEUTRONS: %d\n",TOTAL_NUMBER_OF_NEUTRONS); 
    printf("SCATTERING_COMPONENT: %f\n",SCATTERING_COMPONENT); 
    printf("ABSORBING_COMPONENT: %f\n",ABSORBING_COMPONENT); 
    printf("THICKNESS_OF_THE_PLATE: %d\n",THICKNESS_OF_THE_PLATE); 
    printf("***********************************************************************\n\n"); 
    printf("Total particles reflected,absorbed and transmitted respectively are: %d, %d, %d  \n",no_of_particles_reflected,no_of_particles_absorbed,no_of_particles_transmitted); 
     printf("Percentage of particles reflected,absorbed and transmitted respectively are: %lf, %lf, %lf\n", total_amount_reflected, total_amount_absorbed, total_amount_transmitted); 
    printf("Total time taken in sec %10.6lf\n",diff_t); 
    return 0; 
} 

MPIバージョン:

/** 
    To compile the program:   mpicc -o neutron_transport with_MPI.c -lm 
    To run the program:   mpirun -np <specify no. of processes> neutron_transport 
**/ 
#include<stdio.h> 
#include<math.h> 
#include<stdlib.h> 
#include<time.h> 
#include<mpi.h> 

#define TOTAL_NUMBER_OF_NEUTRONS 1000      //Total number of neutrons 
#define SCATTERING_COMPONENT 0.7       // Scattering component of C 
#define ABSORBING_COMPONENT 0.3       // Absorbing component of C 
#define THICKNESS_OF_THE_PLATE 4       //Thickness of the Plate 
#define NUMBER_OF_ITERATIONS 1000 
#define BLOCK_LOW(id,p,n) ((id)*(n)/(p)) 
#define BLOCK_HIGH(id,p,n) (BLOCK_LOW((id)+1,p,n)-1) 
#define BLOCK_SIZE(id,p,n) (BLOCK_LOW((id)+1,p,n)-BLOCK_LOW(id,p,n))  //Block Size for each process 

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 
}; 

int main(int argc, char *argv[]) 
{ 
    int process_id;         //Process Identifier 
    int number_of_processes;       //Number of Processes in the communicator 
    MPI_Init(&argc,&argv);        //Starting of MPI program 
    MPI_Comm_rank(MPI_COMM_WORLD,&process_id);     //Determine the process ID number in the communicator 
    MPI_Comm_size(MPI_COMM_WORLD,&number_of_processes);    //Determine the number of processes in the communicator 
    int i=0;         //Loop iterator 
    int counter=0; 
    int buffer_share_for_each_process=0; 
    double C ;          // mean distance between neutron/atom interactions is 1/C 
    C=SCATTERING_COMPONENT + ABSORBING_COMPONENT; 
    int flag= 1;              
    double random_number_generator;         // Uniform random number 
    double distance_travelled_before_collision;       // Distance travelled by neutron before collision 
    struct neutron nt; 
    double total_amount_absorbed,total_amount_reflected,total_amount_transmitted; 
    int count_of_reflected_absorbed_transmitted[3]; 
    count_of_reflected_absorbed_transmitted[0]=0; 
    count_of_reflected_absorbed_transmitted[1]=0; 
    count_of_reflected_absorbed_transmitted[2]=0; 
    int global_count_of_reflected_absorbed_transmitted[3]; 
    buffer_share_for_each_process= BLOCK_SIZE(process_id,number_of_processes,TOTAL_NUMBER_OF_NEUTRONS); 
    double elapsed_time1;     //Timer variables 
    elapsed_time1= -MPI_Wtime(); 
    if(number_of_processes > TOTAL_NUMBER_OF_NEUTRONS) 
    { 
     if(!process_id) 
      printf("Too many processes. Kindly provide lesser number of processes.\n"); 
     MPI_Finalize(); 
     exit(1); 
    } 

    for(int i=1;i<=buffer_share_for_each_process;i++) 
    { 
     random_number_generator=rand()%TOTAL_NUMBER_OF_NEUTRONS; 
     random_number_generator=1/random_number_generator; 
     nt.d=0; 
     nt.x=0; 
     flag=1; 
     distance_travelled_before_collision=-(1/C)*log(random_number_generator); 
     nt.x=distance_travelled_before_collision*cos(nt.d); 

     while(flag) 
       { 
      distance_travelled_before_collision=-(1/C)*log(random_number_generator); 
      nt.x=distance_travelled_before_collision*cos(nt.d); 

      if (nt.x<0) 
      { 
       count_of_reflected_absorbed_transmitted[0]++; 
       flag=0; 
      } 
      else if(nt.x>=THICKNESS_OF_THE_PLATE) 
      { 
       count_of_reflected_absorbed_transmitted[2]++; 
       flag=0; 
      } 
      else 
      { 
       if(random_number_generator<(ABSORBING_COMPONENT/C)) 
       { 
        count_of_reflected_absorbed_transmitted[1]++; 
        flag=0; 
       } 
       else 
       { 
        nt.d=random_number_generator*(22.0/7.0); 
        counter++; 
        if(counter==NUMBER_OF_ITERATIONS) 
        { 
         count_of_reflected_absorbed_transmitted[1]++; 
         flag=0; 
         counter=0; 
        } 
       } 
      } 
     } 
    } 
    for(int i= 0 ; i < 3 ; i ++) 
     MPI_Reduce(&count_of_reflected_absorbed_transmitted[i], &global_count_of_reflected_absorbed_transmitted[i], 1, MPI_INT,MPI_SUM,0, MPI_COMM_WORLD) ; 
    elapsed_time1+= MPI_Wtime(); 
    MPI_Barrier(MPI_COMM_WORLD); 
    if(!process_id) 
    { 
     total_amount_reflected=((double)global_count_of_reflected_absorbed_transmitted[0]/(double)TOTAL_NUMBER_OF_NEUTRONS)*100.0; 
     total_amount_absorbed=((double)global_count_of_reflected_absorbed_transmitted[1]/(double)TOTAL_NUMBER_OF_NEUTRONS)*100.0; 
     total_amount_transmitted=((double)global_count_of_reflected_absorbed_transmitted[2]/(double)TOTAL_NUMBER_OF_NEUTRONS)*100.0; 

     printf("\nDATA VALUE ASSUMED:-\n"); 
     printf("*****************************************************************************************\n"); 
     printf("TOTAL_NUMBER_OF_NEUTRONS: %d\n",TOTAL_NUMBER_OF_NEUTRONS); 
     printf("SCATTERING_COMPONENT: %f\n",SCATTERING_COMPONENT); 
     printf("ABSORBING_COMPONENT: %f\n",ABSORBING_COMPONENT); 
     printf("THICKNESS_OF_THE_PLATE: %d\n",THICKNESS_OF_THE_PLATE); 
     printf("TOTAL NUMBER OF PROCESSES: %d\n",number_of_processes); 
     printf("*****************************************************************************************\n\n"); 
     printf("Total particles reflected,absorbed and transmitted respectively are: %d, %d, %d \n",global_count_of_reflected_absorbed_transmitted[0],global_count_of_reflected_absorbed_transmitted[1],global_count_of_reflected_absorbed_transmitted[2]); 
      printf("Percentage of particles reflected,absorbed and transmitted respectively are: %lf, %lf, %lf\n", total_amount_reflected, total_amount_absorbed, total_amount_transmitted); 
     printf("Total time taken in sec %10.6f\n",elapsed_time1); 
    } 
    MPI_Finalize(); 
    return 0; 
} 

OpenMPのバージョン:

/** 
    To compile the program:   gcc -fopenmp with_open_mp.c -lm 
    To run the program:   ./a.out 
**/ 
#include<stdio.h> 
#include<math.h> 
#include<stdlib.h> 
#include<time.h> 
#include<omp.h> 

#define TOTAL_NUMBER_OF_NEUTRONS 10000000        //Total number of neutrons 
#define SCATTERING_COMPONENT 0.7         // Scattering component of C 
#define ABSORBING_COMPONENT 0.3         // Absorbing component of C 
#define THICKNESS_OF_THE_PLATE 5         //Thickness of the Plate 
#define NUMBER_OF_ITERATIONS 1000 
#define TOTAL_NUMBER_OF_THREADS 500 
struct neutron 
{ 
    double d;           // Direction of the neutron (measured in radians between 0 and pi) 
    double x;           // Position of particle in plate(0<=x<=H) where H is the thickness of the plate 
}; 

int main(int argc, char *argv[]) 
{ 
    double C ;            // mean distance between neutron/atom interactions is 1/C 
    C=SCATTERING_COMPONENT + ABSORBING_COMPONENT; 
    int no_of_particles_reflected=0, no_of_particles_absorbed=0, no_of_particles_transmitted=0; 
    int flag= 1;              
    double random_number_generator;           // Uniform random number 
    double distance_travelled_before_collision;         // Distance travelled by neutron before collision 
    struct neutron nt; 
    double total_amount_absorbed,total_amount_reflected,total_amount_transmitted; 
    int counter= 0; 
     srand(time(NULL)); 
    clock_t start, end; 
    double diff_t; 
    diff_t = 0.0; 
    //start the clock 
    start= clock(); 
    //int t=omp_get_num_procs(); 
    omp_set_num_threads(TOTAL_NUMBER_OF_THREADS); 
    #pragma omp parallel for reduction(+:no_of_particles_reflected,no_of_particles_absorbed,no_of_particles_transmitted) 
    for(int i=1;i<=TOTAL_NUMBER_OF_NEUTRONS;i++) 
    { 
     random_number_generator=rand()%TOTAL_NUMBER_OF_NEUTRONS; 
     random_number_generator=1/random_number_generator; 
     nt.d=0; 
     nt.x=0; 
     flag=1; 
     while(flag) 
       { 
      distance_travelled_before_collision=-(1/C)*log(random_number_generator); 
      nt.x=distance_travelled_before_collision*cos(nt.d); 

      if (nt.x<0) 
      { 
       no_of_particles_reflected++; 
       flag=0; 
      } 
      else if(nt.x>=THICKNESS_OF_THE_PLATE) 
      { 
       no_of_particles_transmitted++; 
       flag=0; 
      } 
      else 
      { 
       if(random_number_generator<(ABSORBING_COMPONENT/C)) 
       { 
        no_of_particles_absorbed++; 
        flag=0; 
       } 
       else 
       { 
        nt.d=random_number_generator*(22.0/7.0); 
        counter++; 
        if(counter==NUMBER_OF_ITERATIONS) 
        { 

         no_of_particles_absorbed++; 
         flag=0; 
         counter=0; 
        } 
       } 
      } 
     } 
    } 
    //stop the clock 
    end = clock(); 
    diff_t = ((double) (end - start))/CLOCKS_PER_SEC; 
    printf("\nDATA VALUE ASSUMED:-\n"); 
    printf("*****************************************************************************************\n"); 
    printf("TOTAL_NUMBER_OF_NEUTRONS: %d\n",TOTAL_NUMBER_OF_NEUTRONS); 
    printf("SCATTERING_COMPONENT: %f\n",SCATTERING_COMPONENT); 
    printf("ABSORBING_COMPONENT: %f\n",ABSORBING_COMPONENT); 
    printf("THICKNESS_OF_THE_PLATE: %d\n",THICKNESS_OF_THE_PLATE); 
    printf("TOTAL NUMBER OF THREADS: %d\n",TOTAL_NUMBER_OF_THREADS); 
    printf("*****************************************************************************************\n\n"); 
    total_amount_reflected=((double)no_of_particles_reflected/(double)TOTAL_NUMBER_OF_NEUTRONS)*100.0; 
    total_amount_absorbed=((double)no_of_particles_absorbed/(double)TOTAL_NUMBER_OF_NEUTRONS)*100.0; 
    total_amount_transmitted=((double)no_of_particles_transmitted/(double)TOTAL_NUMBER_OF_NEUTRONS)*100.0; 
printf("Total particles reflected,absorbed and transmitted respectively are: %d, %d, %d \n",no_of_particles_reflected,no_of_particles_absorbed,no_of_particles_transmitted); 
     printf("The amount reflected,absorbed and transmitted respectively are: %lf, %lf, %lf\n", total_amount_reflected, total_amount_absorbed, total_amount_transmitted); 
    printf("Total time taken in sec %10.6lf\n",diff_t); 
    return 0; 
} 

ハイブリッドバージョンのMPIとOpenMP:

/** 
    To compile the program:   mpicc -fopenmp mpi_openmp.c -o hybrid -lm 
    To run the program:   mpirun -np <specify number of processes> ./hybrid 
**/ 
#include<stdio.h> 
#include<math.h> 
#include<stdlib.h> 
#include<time.h> 
#include<omp.h> 
#include<mpi.h> 

#define TOTAL_NUMBER_OF_NEUTRONS 10000000        //Total number of neutrons 
#define SCATTERING_COMPONENT 0.7         // Scattering component of C 
#define ABSORBING_COMPONENT 0.3         // Absorbing component of C 
#define THICKNESS_OF_THE_PLATE 3         //Thickness of the Plate 
#define BLOCK_LOW(id,p,n) ((id)*(n)/(p)) 
#define BLOCK_HIGH(id,p,n) (BLOCK_LOW((id)+1,p,n)-1) 
#define BLOCK_SIZE(id,p,n) (BLOCK_LOW((id)+1,p,n)-BLOCK_LOW(id,p,n))    //Block Size for each process 
#define NUMBER_OF_ITERATIONS 1000 
#define TOTAL_NUMBER_OF_THREADS 500         //Total number of threads 
struct neutron 
{ 
    double d;           // Direction of the neutron (measured in radians between 0 and pi) 
    double x;           // Position of particle in plate(0<=x<=H) where H is the thickness of the plate 
}; 

int main(int argc, char *argv[]) 
{ 
    struct neutron nt; 
    int process_id;           //Process Identifier 
    int number_of_processes;         //Number of Processes in the communicator 
    MPI_Init(&argc,&argv);          //Starting of MPI program 
    MPI_Comm_rank(MPI_COMM_WORLD,&process_id);       //Determine the process ID number in the communicator 
    MPI_Comm_size(MPI_COMM_WORLD,&number_of_processes);      //Determine the number of processes in the communicator 

    double C ;            // mean distance between neutron/atom interactions is 1/C 
    C=SCATTERING_COMPONENT + ABSORBING_COMPONENT; 
    int flag= 1;              
    double random_number_generator;           // Uniform random number 
    double distance_travelled_before_collision;         // Distance travelled by neutron before collision 
    int no_of_particles_reflected=0, no_of_particles_absorbed=0, no_of_particles_transmitted=0; 
    double total_amount_absorbed,total_amount_reflected,total_amount_transmitted; 

    int counter= 0; 
    int buffer_share_for_each_process=0; 
     srand(time(NULL)); 

    int count_of_reflected_absorbed_transmitted[3]; 
    count_of_reflected_absorbed_transmitted[0]=0; 
    count_of_reflected_absorbed_transmitted[1]=0; 
    count_of_reflected_absorbed_transmitted[2]=0; 
    int global_count_of_reflected_absorbed_transmitted[3]; 
    buffer_share_for_each_process= BLOCK_SIZE(process_id,number_of_processes,TOTAL_NUMBER_OF_NEUTRONS); 
    double elapsed_time1;     //Timer variables 
    elapsed_time1= -MPI_Wtime(); 

    //int t=omp_get_num_procs(); 
    omp_set_num_threads(TOTAL_NUMBER_OF_THREADS); 
    if((number_of_processes) > TOTAL_NUMBER_OF_NEUTRONS) 
    { 
     if(!process_id) 
      printf("Too many processes. Kindly provide lesser number of processes.\n"); 
     MPI_Finalize(); 
     exit(1); 
    } 

    #pragma omp parallel for 
    for(int i=1;i<=buffer_share_for_each_process;i++) 
    { 
     //printf("Buffer_Share= %d by process %d (thread %d out of %d)\n",buffer_share_for_each_process,process_id,omp_get_thread_num(),omp_get_num_threads()); 
     random_number_generator=rand()%TOTAL_NUMBER_OF_NEUTRONS; 
     random_number_generator=1/random_number_generator; 
     nt.d=0; 
     nt.x=0; 
     flag=1; 
     while(flag) 
       { 
      distance_travelled_before_collision=-(1/C)*log(random_number_generator); 
      nt.x=distance_travelled_before_collision*cos(nt.d); 

      if (nt.x<0) 
      { 
       count_of_reflected_absorbed_transmitted[0]++; 
       flag=0; 
      } 
      else if(nt.x>=THICKNESS_OF_THE_PLATE) 
      { 
       count_of_reflected_absorbed_transmitted[2]++; 
       flag=0; 
      } 
      else 
      { 
       if(random_number_generator<(ABSORBING_COMPONENT/C)) 
       { 
        count_of_reflected_absorbed_transmitted[1]++; 
        flag=0; 
       } 
       else 
       { 
        nt.d=random_number_generator*(22.0/7.0); 
        counter++; 
        if(counter==NUMBER_OF_ITERATIONS) 
        { 
         count_of_reflected_absorbed_transmitted[1]++; 
         flag=0; 
         counter=0; 
        } 
       } 
      } 
     } 
    } 
    for(int i= 0 ; i < 3 ; i ++) 
     MPI_Reduce(&count_of_reflected_absorbed_transmitted[i], &global_count_of_reflected_absorbed_transmitted[i], 1, MPI_INT,MPI_SUM,0, MPI_COMM_WORLD) ; 
    elapsed_time1+= MPI_Wtime(); 
    MPI_Barrier(MPI_COMM_WORLD); 
    if(!process_id) 
    { 
     total_amount_reflected=((double)global_count_of_reflected_absorbed_transmitted[0]/(double)TOTAL_NUMBER_OF_NEUTRONS)*100.0; 
     total_amount_absorbed=((double)global_count_of_reflected_absorbed_transmitted[1]/(double)TOTAL_NUMBER_OF_NEUTRONS)*100.0; 
     total_amount_transmitted=((double)global_count_of_reflected_absorbed_transmitted[2]/(double)TOTAL_NUMBER_OF_NEUTRONS)*100.0; 

     printf("\nDATA VALUE ASSUMED:-\n"); 
     printf("*****************************************************************************************\n"); 
     printf("TOTAL_NUMBER_OF_NEUTRONS: %d\n",TOTAL_NUMBER_OF_NEUTRONS); 
     printf("SCATTERING_COMPONENT: %f\n",SCATTERING_COMPONENT); 
     printf("ABSORBING_COMPONENT: %f\n",ABSORBING_COMPONENT); 
     printf("THICKNESS_OF_THE_PLATE: %d\n",THICKNESS_OF_THE_PLATE); 
     printf("TOTAL NUMBER OF PROCESSES: %d\n",number_of_processes); 
     printf("TOTAL NUMBER OF THREADS: %d\n",TOTAL_NUMBER_OF_THREADS); 
     printf("*****************************************************************************************\n\n"); 
     printf("Total particles reflected,absorbed and transmitted respectively are: %d, %d, %d \n",global_count_of_reflected_absorbed_transmitted[0],global_count_of_reflected_absorbed_transmitted[1],global_count_of_reflected_absorbed_transmitted[2]); 
      printf("Percentage of particles reflected,absorbed and transmitted respectively are: %lf, %lf, %lf\n", total_amount_reflected, total_amount_absorbed, total_amount_transmitted); 
     printf("Total time taken in sec %10.6f\n",elapsed_time1); 
    } 
    MPI_Finalize(); 
    return 0; 
} 
関連する問題