高斯-赛得尔迭代式 c++_(C实现)有限差分法(FDM)求解雷诺方程

1前言

今日所推文章得头条已经对所求解的模型和算法进行了简单的说明，此处不再赘述。相信工科的学生都学习过C，但是真正用C进行数值计算的并不多见。但是，相信大家都听说过C是比较接近底层的语言，拥有更高的执行效率。本文将利用C进行本期文章中的雷诺方程的求解，文中将大量使用C的内存分配特性(malloc 、calloc、free)和指针，C就像一把雕刻刀，没有太多的酷炫特性，灵活与否全看程序员自己的水平。本文程序采用C语言，GCC6.2.0编译器编译。本文程序将会指针漫天，在微信中不便于阅读，如果有读者需要源码，转发本文，截图回复在评论区，并留下邮箱，可以获得本文源码。

2前处理

前处理主要包括求解域设置、网格划分、边界条件设置，以及为了后续求解进行的必要数据准备。在前处理阶段，需要对网格划分参数，边界条件参数进行设置，通过 set_boundary()调用_set_neumann()、_set_boundary()、_set_source()三个子函数对不同类型求解域的网格节点进行标记。

//-------前处理-------

const float P0=100000;

const float Pa=0;

const float Xarea=40;

const float Yarea=40;

const int Nx=40;

const int Ny=40;

const int dim=1681;

const int center_x[2]={0,19};

const int center_y[2]={0,19};

float delta_x=0;

float delta_y=0;

int * classify_field=(int*)calloc(sizeof(int),dim);

//前处理函数声明

void set_boundary();

void _set_neumann();

void _set_boundary();

void _set_source();

void set_boundary()

{

    _set_neumann();

    _set_boundary();

    _set_source();

}

void _set_neumann()

{

    int neumann=3;

    for(int ix=0;ix<=Nx;ix++)

    {

        *(classify_field+ix*(Ny+1))=neumann;

    }

    for(int iy=0;iy<=Ny;iy++)

    {

        *(classify_field+iy)=neumann;

    }

}

void _set_boundary()

{

    int boundary=1;

    for(int iy=0;iy<=Ny;++iy)

    {

        *(classify_field+iy*(Nx+1)+Nx)=boundary;

    }

    for(int ix=0;ix<=Nx;++ix)

    {

        *(classify_field+Nx*(Ny+1)+ix)=boundary;

    }

}

void _set_source()

{

    int source=2;

    (*(classify_field+center_x[0]*(Ny+1)+center_y[0]))=source;

    (*(classify_field+center_x[1]*(Ny+1)+center_y[1]))=source;

}

void get_delta()

{   

    delta_x=Xarea/(float)(Nx);

    delta_y=Yarea/(float)(Ny);

}

3求解器

求解器主要根据偏微分方程构建相邻节点之间的关系，存在多种差分格式可以选择，比如五点差分格式、九点差分格式等。本文选择了五点差分格式进行求解。不同类型的节点其差分格式的构建不同，具体数学过程本文不在赘述，详细过程可以查阅关于偏微分方程的教科书。通过solver()函数，根据不同的节点标记，调用相应的差分方程构建函数_inner_solver()、_dirichlet_solver()和_neumann_solver()，够建差分方程组。

//------求解器-----

//函数声明

void solver(float *,float *);

void _inner_solver(int,int,float*,float*);

void _dirichlet_solver(int,int,float,float*,float*);

void _neumann_solver(int,int,float*,float*);

void gauss(float * ,float *);

int _max_row(float *,int);

void _swap(float *,int,int);

void _row_multiply(float*,int,int,float);

void _swap_1D(float *,int,int);

void _row_multiply_1D(float*,int,int,float);

void solver(float * left,float * right)

{

    for(int ix=0;ix<=Nx;++ix)

    {

        for(int iy=0;iy<=Ny;++iy)

        {

            switch((*(classify_field+ix*(Ny+1)+iy)))

            {

                case 0:

                {

                   _inner_solver(ix,iy,left,right);

                    break;

                }

                case 1:

                {

                    _dirichlet_solver(ix,iy,Pa,left,right);

                    break;

                }

                case 2:

                {

                    _dirichlet_solver(ix,iy,P0,left,right);

                    break;

                }

                case 3:

                {

                    _neumann_solver(ix,iy,left,right);

                    break;

                }

            }

        }

    }

}

void _inner_solver(int ix,int iy,float* left,float* right)

{

    float Alpha_one=1.0/(delta_x*delta_x);

    float Alpha_two=1.0/(delta_y*delta_y);

    float Alpha_three=1.0/(delta_x*delta_x);

    float Alpha_four=1.0/(delta_y*delta_y);

    float Alpha_zero=Alpha_four+Alpha_three+Alpha_two+Alpha_one;

    *(left+(ix + iy*(Nx+1))*dim+ix+1+iy*(Nx+1))=-Alpha_one;

    *(left+(ix + iy*(Nx+1))*dim+ix+(iy+1)*(Nx+1))=-Alpha_two;

    *(left+(ix + iy*(Nx+1))*dim+ix-1+iy*(Nx+1))=-Alpha_three;

    *(left+(ix + iy*(Nx+1))*dim+ix+(iy-1)*(Nx+1))=-Alpha_four;

    *(left+(ix + iy*(Nx+1))*dim+ix+iy*(Nx+1))=Alpha_zero;

    *(right+ix + iy*(Nx+1))=0;

}

void _neumann_solver(int ix,int iy,float* left,float* right)

{

    float Alpha_one=1.0/(delta_x*delta_x);

    float Alpha_two=1.0/(delta_y*delta_y);

    float Alpha_three=1.0/(delta_x*delta_x);

    float Alpha_four=1.0/(delta_y*delta_y);

    float Alpha_zero=Alpha_four+Alpha_three+Alpha_two+Alpha_one;

    if(ix==0)

    {

        *(left+(ix + iy*(Nx+1))*dim+ix+1+iy*(Nx+1))=-Alpha_one;

        *(left+(ix + iy*(Nx+1))*dim+ix+(iy+1)*(Nx+1))=-Alpha_two;

        *(left+(ix + iy*(Nx+1))*dim+ix+(iy-1)*(Nx+1))=-Alpha_four;

        *(left+(ix + iy*(Nx+1))*dim+ix+iy*(Nx+1))=Alpha_zero-Alpha_three;      

    } 

    else if(iy==0)

    {

        *(left+(ix + iy*(Nx+1))*dim+ix+1+iy*(Nx+1))=-Alpha_one;

        *(left+(ix + iy*(Nx+1))*dim+ix+(iy+1)*(Nx+1))=-Alpha_two;

        *(left+(ix + iy*(Nx+1))*dim+ix-1+iy*(Nx+1))=-Alpha_three;

        *(left+(ix + iy*(Nx+1))*dim+ix+iy*(Nx+1))=Alpha_zero-Alpha_four;      

    }

    *(right+ix + iy*(Nx+1))=0;

}

void _dirichlet_solver(int ix,int iy,float pd,float * left,float * right)

{

    *(left+(ix + iy*(Ny+1))*dim+ix+iy*(Nx+1))=1.0;

    *(right+ix + iy*(Nx+1))=pd;

}

4方程组求解

通过各类节点构建起来的差分方程联立，得到差分方程组，求解差分方程组的方法多如牛毛，类似雅克比迭代法、高斯赛德尔迭代法、高斯消元法、多重网格法等。本文采用gauss消元法进行方程组求解，求解过程通过gauss()函数完成，并调用了一系列用于初等变换、选主元等功能的子函数(_max_row()、_swap()、_row_multiply()、_swap_1D()、_row_multiply_1D ())。

void gauss(float * left,float * right)

{

    for(int ix=1;ix<dim;ix++)

    {       

        int max_row=_max_row(left,ix-1);

        if(max_row!=(ix-1))

        {

            _swap(left,ix-1,max_row);

            _swap_1D(right,ix-1,max_row);

        }       

        for(int jx=ix;jx<dim;jx++)

        {

            float temp=(*(left+jx*dim+ix-1))/(*(left+(ix-1)*dim+ix-1));

            _row_multiply(left,ix-1,jx,temp); 

            _row_multiply_1D(right,ix-1,jx,temp);

        }         

    }

    float * result=(float*)calloc(sizeof(float),dim);

    for(int ix=dim-1;ix>=0;ix--)

    {

        float sum=0;

        for(int iy=dim-1;iy>ix;iy--)

        {

            float temp=(*(result+iy));

            sum=(*(left+ix*dim+iy))*temp+sum;

        }

        (*(result+ix))=((*(right+ix))-sum)/(*(left+ix*dim+ix));

    }

    for(int iy=0;iy<dim;iy++)

    {

        (*(right+iy))=(*(result+iy));

    }

    free(result);

}

int _max_row(float * matrix,int col)

{

    float max_value=fabs(*(matrix+col*dim+col));

    int max_row=col;

    for(int ix=col+1;ix<dim;++ix)

    {

        if(fabs(*(matrix+ix*dim+col))>max_value)

        {

            max_row=ix;

            max_value=fabs((*(matrix+ix*dim+col)));

        }     

    }

    return max_row;

}

void _swap(float * matrix,int rowup,int rowdown)

{

    for(int iy=0;iy<dim;iy++)

    {

        double temp=(*(matrix+rowup*dim+iy));

        (*(matrix+rowup*dim+iy))=(*(matrix+rowdown*dim+iy));

        (*(matrix+rowdown*dim+iy))=temp;

    }

}

void _swap_1D(float * matrix,int rowup,int rowdown)

{

        double temp=(*(matrix+rowup));

        (*(matrix+rowup))=(*(matrix+rowdown));

        (*(matrix+rowdown))=temp;

}

void _row_multiply(float* matrix,int firstrow,int secondrow,float a)

{

    for(int iy=0;iy<dim;iy++)

    {

        float temp=(*(matrix+firstrow*dim+iy));

        (*(matrix+secondrow*dim+iy))=(*(matrix+secondrow*dim+iy))-a*temp;      

    }

}

void _row_multiply_1D(float* matrix,int firstrow,int secondrow,float a)

{

        float temp=(*(matrix+firstrow));

        (*(matrix+secondrow))=(*(matrix+secondrow))-a*temp;  

}

5后处理

后处理实际上就是为了将计算结果进行适当的处理，以便直观的看到计算结果。通过求解器得到的计算结果是一个包含(Nx+1)X(Ny+1)个数值的一维向量，为了便于展示，需要将其变成一个(Nx+1)X(Ny+1)矩阵。由于采用了求解域的1/4对称部分进行计算，需要将计算结果进行对称拼接，得到完整的求解域中的结果分布。由于C++中的数值绘图库不太容易上手，所以将计算结果输出为result.txt文件，将result.txt文件导入到Origin或者Matlab软件，就可以轻易的得到计算结果了。在本期头条文章中，采用Matlab实现，网格数为80X80，得到的云图比较细腻，在本文程序中80X80的网格计算机反应太慢，所以采用40X40网格，得到的计算云图略显粗糙。

void grid(float* result)

{

    float * xmesh=(float*)calloc(sizeof(float),(2*Nx+1)*(2*Ny+1));

    float * ymesh=(float*)calloc(sizeof(float),(2*Nx+1)*(2*Ny+1));

    for(int ix=0;ix<=2*Nx;++ix)

    {

        for(int iy=0;iy<=2*Ny;++iy)

        {

            (*(xmesh+ix*(2*Ny+1)+iy))=ix*delta_x;

            (*(ymesh+iy*(2*Nx+1)+ix))=iy*delta_y;

        }

    }

    FILE* xp = fopen("xmesh.txt","w+");

    for(int ix=0;ix<=2*Nx;++ix)

    {

        for(int iy=0;iy<=2*Ny;++iy)

        {

            fprintf(xp,"%8.4f",(*(xmesh+ix*(2*Ny+1)+iy)));

        }

        fprintf(xp,"\n");

    }   

    fclose(xp);

    FILE* yp = fopen("ymesh.txt","w+");

    for(int iy=0;iy<=2*Ny;++iy)

    {

        for(int ix=0;ix<=2*Nx;++ix)

        {

            fprintf(yp,"%8.4f",(*(ymesh+ix*(2*Ny+1)+iy)));

        }

        fprintf(xp,"\n");

    }   

    fclose(yp);

    free(xmesh);

    free(ymesh);

    float * allresult=(float*)calloc(sizeof(float),(2*Nx+1)*(2*Ny+1));

    for(int ix=0;ix<=2*Nx;++ix)

    {

        for(int iy=0;iy<=2*Ny;++iy)

        {

            if (ix<Nx && iy<Ny) (*(allresult+ix*(2*Ny+1)+iy))=(*(result+(Nx-ix)*(Ny+1)+Ny-iy));

            else if(ix<Nx && iy>=Ny) (*(allresult+ix*(2*Ny+1)+iy))=(*(result+(Nx-ix)*(Ny+1)-Ny+iy));

            else if(ix>=Nx && iy<Ny) (*(allresult+ix*(2*Ny+1)+iy))=(*(result+(ix-Nx)*(Ny+1)+Ny-iy));

            else (*(allresult+ix*(2*Ny+1)+iy))=(*(result+(ix-Nx)*(Ny+1)-Ny+iy));

        }

    }

    FILE* re = fopen("result.txt","w+");

    for(int iy=0;iy<=2*Ny;++iy)

    {

        for(int ix=0;ix<=2*Nx;++ix)

        {

            fprintf(re,"%8.4f\t",(*(allresult+iy*(2*Nx+1)+ix)));

        }

        fprintf(re,"\n");

    }   

    fclose(re);

    free(result);

    free(allresult);   

}

计算结果云图如下所示：

6总结

通过编写程序求解简化的雷诺方程，能够更加清晰的体会到各类CAE商用软件中的仿真逻辑，以及其命令流的语言逻辑，甚至可以对其进行二次开发。当遇到CAE软件求解困难的模型的时候，不妨对模型进行适当的简化，通过编程也能够很好的实现仿真的目的。本文仅仅是通过C编程实现，都说C是高效率的语言，Matlab和Python的数值运算库是经过特性优化，且采用的方程组求解算法内置，不知是何种算法，就不参与比对，在同样使用高斯消元法的条件下，C与C++、Fortran的运行时间和可执行程序大小如下所示：

语言	运行时间(s)	可执行文件大小(K)
C	10.1	73
C++(Eigen库)	59.982	7249
C++(allocator)	13.51	166
Fortran	24	77

可以看出C的确是高效精简的语言，采用allocator类编写的C++效率也不错，文件略大一些，Fortran也很精简，但是貌似有点点慢，没有发挥其作为数值计算高效语言的优势，估计是程序没有写好。但是，C++采用Eigen库就很尴尬了，Eigen库对于小程序来说过于庞大，导致生成的可执行文件大了几十倍，执行效率很低。本期的其他文章将会对同样的问题采用其他编程语言(Matlab，C++，Python，Fortran)进行实现，欢迎查阅。

读者需要源码，转发本文，截图回复在微信后台，并留下邮箱，可以获得本文源码。

作者/编辑：剑花潇湘

打造全国最自由

最开放的技术交流平台

同时欢迎技术能手加入讲师团队、工程师团队

优酷频道“ftc正青春”发布有各软件原创视频

感谢，一路有你