[原创]Ausm格式的简单源代码

iceArcher

    自己编的很简单的CFD二维求解程序，原先是学习Ausm格式用的。说它简单，因为除了最必须的东东，其余一概没有，刚好拿来玩。况且这个结果的精度不高。至于计算结果的正确性如何，偶也8晓得啦。*^_^*
    大家有空看看，最好对程序提些意见。可别用鸡蛋砸偶呵。
//---------------------------------------------------------------------------
//filename: TEuler2d-mod.cpp
//written by : Dr.Tinghe Lee
//this is additional cpp file to assume the right of my code
//a very neat code for ausm scheme test usage!
//---------------------------------------------------------------------------
#include <fstream>
#include <iostream>
#include <cmath>
#define gama 1.4
using namespace std;
//---------------------------------------------------------------------------
// predeclare variables
int im,jm,imax,jmax,it,itmax;
double rm,cfl;
double *x,*y,*rr,*ru,*rv,*re,*pp,*flxrj,*flxuj,*flxvj,*flxej,*dtime;
double *vol,*rrn,*run,*rvn,*ren,*flxri,*flxui,*flxvi,*flxei;
double *rrt,*rut,*rvt,*ret,p0,r0,u0,v0,e0,c0,q0;
//---------------------------------------------------------------------------
// predeclare functions
void input(),initflow(),bound(),rnk(),resid(),output(),farbnd(),walbnd(),
     timestep(),tecgrid();
//---------------------------------------------------------------------------
template<class Type> Type mmin(Type a,Type b)
{
    return (a<=b)?a:b;
}
//---------------------------------------------------------------------------
template<class Type> Type mmax(Type a,Type b)
{
    return (a>=b)?a:b;
}
//---------------------------------------------------------------------------
int main()
{
    rm  = 3.0;
    cfl = 1.1;
    itmax = 2000;
   
    input();
    tecgrid();
    initflow();
    bound();
    for(it=0;it<itmax;it++)
    {
        rnk();
        resid();
    }
    output();
    delete [] vol; delete [] pp; delete [] rr; delete [] ru; delete [] rv; delete [] re;
    delete [] rrn; delete [] run; delete [] rvn; delete [] ren; delete [] dtime;
    delete [] rrt; delete [] rut; delete [] rvt; delete [] ret;
    delete [] flxri; delete [] flxui; delete [] flxvi; delete [] flxei;
    delete [] flxrj; delete [] flxuj; delete [] flxvj; delete [] flxej;
    delete [] x; delete [] y;
    return 0;
}
//---------------------------------------------------------------------------
void input()
{
    int i,j,ii,i1,i2,i3,i4;
    double rx1,rx2,ry1,ry2,vot;
    fstream mfile;
    mfile.open("gri2d.dat",ios::in);
    mfile>> im >> jm;
    imax = im+2;
    jmax = jm+2;
    x = new double [imax*jmax];
    y = new double [imax*jmax];
    for(i=1;i<=im;i++)
    for(j=1;j<=jm;j++)
    {
        ii = i*jmax+j;
        mfile>> x[ii]>>y[ii];
    }
    mfile.close();
    vol = new double [imax*jmax];
    rr  = new double [imax*jmax];
    ru  = new double [imax*jmax];
    rv  = new double [imax*jmax];
    re  = new double [imax*jmax];
    pp  = new double [imax*jmax];
    flxri = new double [imax*jmax];
    flxui = new double [imax*jmax];
    flxvi = new double [imax*jmax];
    flxei = new double [imax*jmax];
    flxrj = new double [imax*jmax];
    flxuj = new double [imax*jmax];
    flxvj = new double [imax*jmax];
    flxej = new double [imax*jmax];
    rrn  = new double [imax*jmax];
    run  = new double [imax*jmax];
    rvn  = new double [imax*jmax];
    ren  = new double [imax*jmax];
    rrt  = new double [imax*jmax];
    rut  = new double [imax*jmax];
    rvt  = new double [imax*jmax];
    ret  = new double [imax*jmax];
    dtime = new double [imax*jmax];
   
    for(i=1;i<im+1;i++)
    for(j=1;j<jm+1;j++)
    {
        i1 = i*jmax+j;
        i2 = (i+1)*jmax+j;
        i3 = (i+1)*jmax+j+1;
        i4 = i*jmax+j+1;
        rx1 = x[i3]-x[i1];
        ry1 = y[i3]-y[i1];
        rx2 = x[i4]-x[i1];
        ry2 = y[i4]-y[i1];
        vot = rx1*ry2-rx2*ry1;
        rx2 = rx1;
        ry2 = ry1;
        rx1 = x[i2]-x[i1];
        ry1 = y[i2]-y[i1];
        vol[i1] = 0.5*(rx1*ry2-rx2*ry1+vot);
    }
}
//---------------------------------------------------------------------------
void initflow()
{
    int i,j,ii;
    p0 = 1;
    r0 = 1;
    c0 = sqrt(gama*p0/r0);
    q0 = rm*c0;
    u0 = q0;
    v0 = 0.0;
    e0 = p0/((gama-1.0)*r0)+q0*q0*0.5;
    for(i=0;i<imax;i++)
    for(j=0;j<jmax;j++)
    {
        ii = i*jmax+j;
        rr[ii] = r0;
        ru[ii] = r0*u0;
        rv[ii] = r0*v0;
        re[ii] = r0*e0;
        pp[ii] = p0;
    }
}
//---------------------------------------------------------------------------
void farbnd()
{
    int i,j,ii,i1;
    for(i=0;i<=0;i++)
    for(j=1;j<jm;j++)
    {
        ii = i*jmax+j;
        rr[ii] = r0;
        ru[ii] = r0*u0;
        rv[ii] = r0*v0;
        re[ii] = r0*e0;
        pp[ii] = p0;
    }
    for(i=im;i<=im;i++)
    for(j=1;j<jm;j++)
    {
        ii = i*jmax+j;
        i1 = (i-1)*jmax+j;
        rr[ii] = rr[i1];
        ru[ii] = ru[i1];
        rv[ii] = rv[i1];
        re[ii] = re[i1];
        pp[ii] = pp[i1];
    }
}
//---------------------------------------------------------------------------
void walbnd()
{
    int i,j,ii,i1;
    for(i=1;i<im;i++)
    for(j=0;j<=0;j++)
    {
        ii = i*jmax+j;
        i1 = i*jmax+j+1;
        rr[ii] = 0;
        ru[ii] = -ru[i1];
        rv[ii] = -rv[i1];
        re[ii] = 0;
        pp[ii] = pp[i1];
    }
    for(i=1;i<im;i++)
    for(j=jm;j<=jm;j++)
    {
        ii = i*jmax+j;
        i1 = i*jmax+j-1;
        rr[ii] = 0;
        ru[ii] = -ru[i1];
        rv[ii] = -rv[i1];
        re[ii] = 0;
        pp[ii] = pp[i1];
    }
}
//---------------------------------------------------------------------------
void flux()
{
    int i,j,i1,i2,i3,i4,il,ir;
    double rrl,rul,rvl,rel,rhl,rrr,rur,rvr,rer,rhr,ppl,ppr;
    double xn,yn,al,ar,area,ql,qr,ml,mr,ms,pl,pr,alf,arf,ps,as;
   
    for(i=1;i<=im;i++)
    for(j=1;j<jm;j++)
    {
        i1 = i*jmax+j;
        i2 = (i+1)*jmax+j;
        i3 = (i+1)*jmax+j+1;
        i4 = i*jmax+j+1;
        il = (i-1)*jmax+j;
        ir = i1;
        xn = y[i4]-y[i1];
        yn = x[i1]-x[i4];
        area = sqrt(xn*xn+yn*yn);
        rrl = rr[il];
        rul = ru[il];
        rvl = rv[il];
        rel = re[il];
        ppl = pp[il];
        rhl = rel+ppl;
        rrr = rr[ir];
        rur = ru[ir];
        rvr = rv[ir];
        rer = re[ir];
        ppr = pp[ir];
        rhr = rer+ppr;
        al = sqrt(2.0*(gama-1.0)*rhl/rrl/(gama+1.0));
        ar = sqrt(2.0*(gama-1.0)*rhr/rrr/(gama+1.0));
        ql = (rul*xn+rvl*yn)/rrl/area;
        qr = (rur*xn+rvr*yn)/rrr/area;
        al = al*al/mmax(fabs(ql),al);
        ar = ar*ar/mmax(fabs(qr),ar);
        as = mmin(al,ar);
        ml = ql/as;
        mr = qr/as;
        if(fabs(ml)<1.0) ml = 0.25*(ml+1.0)*(ml+1.0);
        if(ml<=-1.0) ml = 0.0;
        if(fabs(mr)<1.0) mr = -0.25*(mr-1.0)*(mr-1.0);
        if(mr>=1.0) mr = 0.0;
        pl = ppl;
        pr = ppr;
        if(fabs(ml)<1.0) pl = 0.25*ppl*(ml+1.0)*(ml+1.0)*(2.0-ml);
        if(ml<=-1.0) pl = 0.0;
        if(fabs(mr)<1.0) pr = 0.25*ppr*(mr-1.0)*(mr-1.0)*(2.0+mr);
        if(mr>=1.0) pr = 0.0;
        ps  = pl+pr;
        ms  = ml+mr;
        alf = 0.5*ms*area;
        arf = 0.5*(fabs(ms)+0.5*(mr-1.0)*(mr-1.0))*area;
        flxri[i1] = alf*(rrl*al+rrr*ar)-arf*(rrr*ar-rrl*al);
        flxui[i1] = alf*(rul*al+rur*ar)-arf*(rur*ar-rul*al)+ps*xn;
        flxvi[i1] = alf*(rvl*al+rvr*ar)-arf*(rvr*ar-rvl*al)+ps*yn;
        flxei[i1] = alf*(rhl*al+rhr*ar)-arf*(rhr*ar-rhl*al);
    }
    for(j=2;j<jm;j++)
    for(i=1;i<im;i++)
    {
        i1 = i*jmax+j;
        i2 = (i+1)*jmax+j;
        i3 = (i+1)*jmax+j+1;
        i4 = i*jmax+j+1;
        il = i*jmax+j-1;
        ir = i1;
        xn = y[i1]-y[i2];
        yn = x[i2]-x[i1];
        area = sqrt(xn*xn+yn*yn);
        rrl = rr[il];
        rul = ru[il];
        rvl = rv[il];
        rel = re[il];
        ppl = pp[il];
        rhl = rel+ppl;
        rrr = rr[ir];
        rur = ru[ir];
        rvr = rv[ir];
        rer = re[ir];
        ppr = pp[ir];
        rhr = rer+ppr;
        al = sqrt(2.0*(gama-1.0)*rhl/rrl/(gama+1.0));
        ar = sqrt(2.0*(gama-1.0)*rhr/rrr/(gama+1.0));
        ql = (rul*xn+rvl*yn)/rrl/area;
        qr = (rur*xn+rvr*yn)/rrr/area;
        al = al*al/mmax(fabs(ql),al);
        ar = ar*ar/mmax(fabs(qr),ar);
        as = mmin(al,ar);
        ml = ql/as;
        mr = qr/as;
        if(fabs(ml)<1.0) ml = 0.25*(ml+1.0)*(ml+1.0);
        if(ml<=-1.0) ml = 0.0;
        if(fabs(mr)<1.0) mr = -0.25*(mr-1.0)*(mr-1.0);
        if(mr>=1.0) mr = 0.0;
        pl = ppl;
        pr = ppr;
        if(fabs(ml)<1.0) pl = 0.25*ppl*(ml+1.0)*(ml+1.0)*(2.0-ml);
        if(ml<=-1.0) pl = 0.0;
        if(fabs(mr)<1.0) pr = 0.25*ppr*(mr-1.0)*(mr-1.0)*(2.0+mr);
        if(mr>=1.0) pr = 0.0;
        ps  = pl+pr;
        ms  = ml+mr;
        alf = 0.5*ms*area;
        arf = 0.5*(fabs(ms)+0.5*(mr-1.0)*(mr-1.0))*area;
   
        flxrj[i1] = alf*(rrl*al+rrr*ar)-arf*(rrr*ar-rrl*al);
        flxuj[i1] = alf*(rul*al+rur*ar)-arf*(rur*ar-rul*al)+ps*xn;
        flxvj[i1] = alf*(rvl*al+rvr*ar)-arf*(rvr*ar-rvl*al)+ps*yn;
        flxej[i1] = alf*(rhl*al+rhr*ar)-arf*(rhr*ar-rhl*al);
    }
    for(j=1;j<=1;j++)
    for(i=1;i<im;i++)
    {
        i1 = i*jmax+j;
        i2 = (i+1)*jmax+j;
        xn = y[i1]-y[i2];
        yn = x[i2]-x[i1];
        ps  = pp[i1];
        flxrj[i1] =0;
        flxuj[i1] =ps*xn;
        flxvj[i1] =ps*yn;
        flxej[i1] =0;
    }
    for(j=jm;j<=jm;j++)
    for(i=1;i<im;i++)
    {
        i1 = i*jmax+j;
        i2 = (i+1)*jmax+j;
        i3 = i*jmax+j-1;
        xn = y[i1]-y[i2];
        yn = x[i2]-x[i1];
        ps  = pp[i3];
        flxrj[i1] =0;
        flxuj[i1] =ps*xn;
        flxvj[i1] =ps*yn;
        flxej[i1] =0;
    }
    for(i=1;i<im;i++)
    for(j=1;j<jm;j++)
    {
        i1 = i*jmax+j;
        i2 = (i+1)*jmax+j;
        i4 = i*jmax+j+1;
        rrn[i1] = -flxri[i1]+flxri[i2]-flxrj[i1]+flxrj[i4];
        run[i1] = -flxui[i1]+flxui[i2]-flxuj[i1]+flxuj[i4];
        rvn[i1] = -flxvi[i1]+flxvi[i2]-flxvj[i1]+flxvj[i4];
        ren[i1] = -flxei[i1]+flxei[i2]-flxej[i1]+flxej[i4];
    }
}
//---------------------------------------------------------------------------
void rnk()
{
    int i,j,n,ii;
    double rq;
    double alf[5]={0.25,1.0/6.0,3.0/8.0,0.5,1.0};
    for(i=1;i<im;i++)
    for(j=1;j<jm;j++)
    {
        ii = i*jmax+j;
        rrt[ii] = rr[ii];
        rut[ii] = ru[ii];
        rvt[ii] = rv[ii];
        ret[ii] = re[ii];
    }
    for(n=0;n<5;n++)
    {
        timestep();
        flux();
        for(i=1;i<im;i++)
        for(j=1;j<jm;j++)
        {
            ii = i*jmax+j;
            rr[ii] = rrt[ii] - alf[n]*dtime[ii]*rrn[ii]/vol[ii];
            ru[ii] = rut[ii] - alf[n]*dtime[ii]*run[ii]/vol[ii];
            rv[ii] = rvt[ii] - alf[n]*dtime[ii]*rvn[ii]/vol[ii];
            re[ii] = ret[ii] - alf[n]*dtime[ii]*ren[ii]/vol[ii];
            rq = ru[ii]*ru[ii]+rv[ii]*rv[ii];
            pp[ii] = (gama-1.0)*(re[ii]-rq/rr[ii]*0.5);
            if(pp[ii]<0||rr[ii]<0)
            {
                rr[ii] = r0;
                ru[ii] = r0*u0;
                rv[ii] = r0*v0;
                re[ii] = r0*e0;
                pp[ii] = p0;
            }
        }
        bound();
    }
}
//---------------------------------------------------------------------------
void bound()
{
    farbnd();
    walbnd();
}
//---------------------------------------------------------------------------
void timestep()
{
    int i,j,i1,i2,i3,i4;
    double xt,yt,xn,yn,a1,a2,c,uu,vv;
    for(i=1;i<im;i++)
    for(j=1;j<jm;j++)
    {
        i1 = i*jmax+j;
        i2 = (i+1)*jmax+j;
        i3 = (i+1)*jmax+j+1;
        i4 = i*jmax+j+1;
        c  = sqrt(gama*pp[i1]/rr[i1]);
        xt = 0.5*(x[i4]-x[i1]+x[i3]-x[i2]);
        yt = 0.5*(y[i4]-y[i1]+y[i3]-y[i2]);
        xn = yt;
        yn = -xt;
        a1 = sqrt(xn*xn+yn*yn);
        uu = (ru[i1]*xn+rv[i1]*yn)/rr[i1]/a1;
        xt = 0.5*(x[i2]-x[i1]+x[i3]-x[i4]);
        yt = 0.5*(y[i2]-y[i1]+y[i3]-y[i4]);
        xn = -yt;
        yn =  xt;
        a2 = sqrt(xn*xn+yn*yn);
        vv = (ru[i1]*xn+rv[i1]*yn)/rr[i1]/a2;
        dtime[i1] = cfl*vol[i1]/(fabs(uu)+fabs(vv)+c*(a1+a2));
    }
}
//---------------------------------------------------------------------------
void resid()
{
    int i,j,ii,ierrm,jerrm;
    double err,errmax,dskip;
    err = 0;
    errmax = 0;
    for(i=1;i<im;i++)
    for(j=1;j<jm;j++)
    {
        ii = i*jmax+j;
        dskip = rr[ii]-rrt[ii];
        err = dskip*dskip+err;
        if(fabs(dskip)>errmax)
        {
            errmax = fabs(dskip);
            ierrm = i;
            jerrm = j;
        }
    }
    err = sqrt(err/(im-1)/(jm-1));
    cout<<"it="<<it<<" err="<<err<<" irm="<<ierrm<<" jrm="<<jerrm
        <<" errm="<<errmax<<endl;
}
//---------------------------------------------------------------------------
void output()
{
    int i,j,i1,i2,i3,i4;
    double uu,vv,qq,cc,ma;
    double *xc,*yc;
    xc = new double[imax*jmax];
    yc = new double[imax*jmax];
    for(i=1;i<im;i++)
    for(j=1;j<jm;j++)
    {
        i1 = i*jmax+j;
        i2 = (i+1)*jmax+j;
        i3 = (i+1)*jmax+j+1;
        i4 = i*jmax+j+1;
        xc[i1] = 0.25*(x[i1]+x[i2]+x[i3]+x[i4]);
        yc[i1] = 0.25*(y[i1]+y[i2]+y[i3]+y[i4]);
    }
    fstream mfile;
    mfile.open("tecress.dat",ios:ut);
    mfile<<"variables=\"x\",\"y\",\"ma\",\"r\",\"p\",\"u\",\"v\"";
    mfile<<endl;
    mfile<<"zone i="<<im-1<<" j="<<jm-1<<" f=point"<<endl;
    for(j=1;j<jm;j++)
    for(i=1;i<im;i++)
    {
        i1 = i*jmax+j;
        uu = ru[i1]/rr[i1];
        vv = rv[i1]/rr[i1];
        qq = sqrt(uu*uu+vv*vv);
        cc = sqrt(gama*pp[i1]/rr[i1]);
        ma = qq/cc;
        mfile<<xc[i1]<<" "<<yc[i1]<<" "<<ma<<" "<<rr[i1]<<" ";
        mfile<<pp[i1]<<" "<<ru[i1]/rr[i1]<<" "<<rv[i1]/rr[i1]<<endl;
    }
    mfile.close();
    delete [] xc;
    delete [] yc;
}
//---------------------------------------------------------------------------
void tecgrid()
{
    int i,j,ii;
    fstream mfile;
    mfile.open("tecMe.dat",ios:ut);
    mfile<<"variables=\"x\",\"y\""<<endl;
    mfile<<"zone i="<<im<<" j="<<jm<<" f=point"<<endl;
    for(j=1;j<=jm;j++)
    for(i=1;i<=im;i++)
    {
        ii = i*jmax+j;
        mfile<<x[ii]<<" "<<y[ii]<<endl;
    }
    mfile.close();
}
//---------------------------------------------------------------------------

boling2000

很不错呀，编程风格也很好，看来功底很深厚！！！
佩服中。。。

tigerwoods

不错，不错。
虽然激波形状很粗，但网格也稀。

solomonzj

大侠能不能对 Ausm 格式，进行一点介绍亚！
渴望了解！有什么优点吗？
谢谢了

songke21

楼主，我也在搞AUSM格式，不少地方弄不明白。请问你的格式里用的是什么限制器？我想用Van Albada限制器，但是其具体形式看不太明白，请指教。

wangdingxi

最近我也了解了一下AUSM格式，AUSM是Advection Upstream Splitting Method,现在AUSM已经发展成一系列格式，如AUSM+,AUSMD,AUSMV,AUSMDV等等。在原理上AUSM格式和Van Leer flux splitting差不多，都是属于矢通量分裂。在AUSM中，将压力进行了分裂，但是Van Leer flux splitting没有。详细可以看看这个链接http://www.chimeracfd.com/java/gryphon/fluxausm.html

laser2k

[这个贴子最后由laser2k在 2006/01/06 05:04pm 第 1 次编辑]

   

zhvickie

有没有例子数据？

oldcat

楼主能把程序拿出来共享，应该鼓励。
我也觉得例子没有算好，所以看了一下程序，发现这个程序只能用来算正交网格，lz给出的例子不恰当，因为网格不是正交的。

gotowuhan

鼓励原创

applytina

这个程序是不是没贴完，还是要从别的地方看。正在学习ausm格式。对其中一些点很是不理解。

kaiqun98

感谢楼主分享