c - 中子输运模型的难点

Question

我正在使用蒙特卡洛方法模拟中子传输模型。我首先要实现 M. J. Quinn 的《parallel programming with openmp and MPI》一书中给出的顺序算法

下面给出的(顺序)代码是中子输运模型的模拟。当我第一次运行可执行文件时，它在无限循环中运行。当我停止并再次运行时，我得到的输出是 0、0、0，并且不是输入参数所有可能值的预期输出，如 M.J.Quinn parallel programming with openmp 的解决方案手册中提到的问题中给出的那样MPI(Q.No 10.8)。

本书的软拷贝如下: https://docs.google.com/document/d/19aKs6XF3Gt6BRGSGVysZqam3Cb97rgRvTCYP2Nw-C5M/edit

预期输出为:反射量、吸收量和透射率分别为:25.31、46.13、28.56。编译代码时没有报错。请帮帮我。这是我写的以下代码:

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 1

你有

o=(r/n)*100;

在乘法之前执行整数除法。由于您的预期输出是 25.31 这意味着 r/n 的结果应该是 0.2531，但是这样的整数除法的结果是 0 支持您得到的 0。

我建议这样做

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

这将执行正确的算法。

此外，您还使用了未初始化的变量。