Finite Volume LBM Cell Vertex Scheme

The code below is posted as a support for the following thread on the forum.

/*
FVLBM code
 
Problem:
 Poiseuille flow
 
Assumptions:
1- this model is cell vertex finite volume metohd
2- each vertex  has 8 neighbors
3- dx is the same  between all vertices
 
        4  3  2
         \ | /
          \|/
        5--0--1
          /|\
         / | \
        6  7  8
 
this is the numbering of neighbors and also the distributio functions
 
Pnt is a help sruct to hold a pier of data
*/
#include <iostream>
#include <fstream>
#include <math.h>
using namespace std;
 
template <typename T>struct Pnt
{
    Pnt(){v[0]=0,v[1]=0;}
    Pnt(T v1, T v2) { v[0]=v1,v[1]=v2;}
    T v[2];
    T& operator[](int index){ return v[index];}
};
 
void Q9(int& q,int**& c, double*& t, int*& ops, double& csq )
{
    q=9;
    t=new double[q];
    ops=new int[q];
    c=new int*[2];
    for (int iD=0; iD<2; ++iD)
        c[iD]=new int[q];
    t[0]=4./9.; t[1]=1./9.; t[2]=1./9.; t[3]=1./36.; t[4]=1./36;
                t[5]=1./9.; t[6]=1./9.; t[7]=1./36.; t[8]=1./36;
    c[0][0]=0;      c[0][1]=-1;     c[0][2]= 0; c[0][3]=-1;     c[0][4]=-1;
    c[1][0]=0;      c[1][1]= 0;     c[1][2]=-1; c[1][3]= 1;     c[1][4]=-1;
 
                    c[0][5]=+1;     c[0][6]= 0; c[0][7]=+1;     c[0][8]=+1;
                    c[1][5]= 0;     c[1][6]=+1; c[1][7]=-1;     c[1][8]=+1;
    ops[0]=0;
    ops[1]=5;
    ops[2]=6;
    ops[3]=7;
    ops[4]=8;
    ops[5]=1;
    ops[6]=2;
    ops[7]=3;
    ops[8]=4;
    csq=1./3.;
}
 
void Q9v(int& q,int**& c, double*& t, int*& ops, double& csq )
{
    q=9;
    t=new double[q];
    ops=new int[q];
    c=new int*[2];
    for (int iD=0; iD<2; ++iD)
        c[iD]=new int[q];
    t[0]=4./9.; t[1]=1./9.; t[2]=1./36.; t[3]=1./9.; t[4]=1./36;
                t[5]=1./9.; t[6]=1./36.; t[7]=1./9.; t[8]=1./36;
    c[0][0]=0;      c[0][1]=+1;     c[0][2]=+1; c[0][3]= 0;     c[0][4]=-1;
    c[1][0]=0;      c[1][1]= 0;     c[1][2]=+1; c[1][3]= 1;     c[1][4]=+1;
 
                    c[0][5]=-1;     c[0][6]=-1; c[0][7]= 0;     c[0][8]=+1;
                    c[1][5]= 0;     c[1][6]=-1; c[1][7]=-1;     c[1][8]=-1;
    ops[0]=0;
    ops[1]=5;
    ops[2]=6;
    ops[3]=7;
    ops[4]=8;
    ops[5]=1;
    ops[6]=2;
    ops[7]=3;
    ops[8]=4;
    csq=1./3.;
}
 
class Core_2D
{
    private:
    int half ;
    int nextX ,nextY;
    double*** f;
    double*** ftemp;
    //int nX,nY;
    int q;
    int d;
    int** c;
    double* t;
    int* ops;
    double csq;
    double temp;
 
    ofstream out;
 
    double** Rho, **U [2];
 
 
    double ( Core_2D::*Feq ) ( int iF, double rho, Pnt<double> u, double uSqr );
 
    Pnt<double> ** P;
 
    public:
    double Omega;
    double Visc;
    int nX,nY,nT;
 
    Core_2D(int I, int J, int** ck,double* wk, int* opsit ,double& cs, int Q , int type)
    {
        out.open("out.txt",ios::trunc);
        d=2;
        q=Q;
        half = (Q-1)/2;
        nX=I;
        nY=J;
        csq=cs;
 
        f  =new double** [q];
        ftemp  =new double** [q];
        ops=new int      [q];
        for(int iF=0;iF<q;iF++)
        {
            f[iF]=new double* [nX];
            ftemp[iF]=new double* [nX];
            for(int iX=0;iX<nX;iX++)
            {
                f[iF][iX]=new double [nY];
                ftemp[iF][iX]=new double [nY];
            }
        }
 
        c=new int* [d];
 
        for(int iD=0;iD<d;iD++)
            c[iD]=ck[iD];
 
        t=new double [q];
        for(int iF=0;iF<q;iF++)
        {
            t[iF]=wk[iF];
            ops[iF]=opsit[iF];
        }
 
        if ( type==1 )
			Feq= &Core_2D::fEq1;
		else if ( type==2 )
			Feq= &Core_2D::fEq2;
		else if ( type==4 )
			Feq= &Core_2D::fEq4;
 
        Rho=new double* [nX];
        U[0]=new double* [nX];
        U[1]=new double* [nX];
        P=new Pnt<double>* [nX];
        for(int iX=0;iX<nX;iX++)
        {
            Rho[iX]=new double [nY];
            U[0][iX]=new double [nY];
            U[1][iX]=new double [nY];
            P[iX]=new Pnt<double> [nY];
        }
    }
 
    void Init_Geo(double L, double dx)
    {
        for(int iX=0;iX<nX;iX++)
        {
            for(int iY=0;iY<nY; iY++)
            {
                P[iX][iY][0]=iX*dx;
                P[iX][iY][1]=iY*dx;
            }
        }
    }
 
// initual all the field with rho and u
    void Init (double rho, Pnt<double> u)
    {
        for (int iX=0; iX<nX; ++iX)
            for (int iY=0; iY<nY; ++iY)
            {
                for (int iF=0; iF<q; ++iF)
                    f[iF][iX][iY]= (this->*Feq)(iF, rho, u,u[0]*u[0]+u[1]*u[1]);
                Rho[iX][iY]=rho;
                U[0][iX][iY]=u[0];
                U[1][iX][iY]=u[1];
            }
 
    }
 
// initual all the field with rho
    void Init (double rho)
    {
        for (int iX=0; iX<nX; ++iX)
            for (int iY=0; iY<nY; ++iY)
            {
                for (int iF=0; iF<q; ++iF)
                    f[iF][iX][iY]= rho*t[iF];
                Rho[iX][iY]=rho;
            }
 
    }
 
    void Init (Pnt<int> ix, Pnt<int> iy, double rho, Pnt<double> u)
    {
        for (int iX=ix[0]; iX<=ix[1]; ++iX)
            for (int iY=iy[0]; iY<=iy[1]; ++iY)
            {
                for (int iF=0; iF<q; ++iF)
                    f[iF][iX][iY]= (this->*Feq)(iF, rho, u,u[0]*u[0]+u[1]*u[1]);
                Rho[iX][iY]=rho;
                U[0][iX][iY]=u[0];
                U[1][iX][iY]=u[1];
            }
    }
 
    void Init (Pnt<int> ix, Pnt<int> iy, double rho)
    {
        for (int iX=ix[0]; iX<=ix[1]; ++iX)
            for (int iY=iy[0]; iY<=iy[1]; ++iY)
                Init (rho);
    }
 
    void computeRho(int iX, int iY, double& rho)
    {
        rho = 0.;
        for (int iF=0; iF<q; ++iF)
        {
            rho =rho +f[iF][iX][iY];
        }
    }
 
    void computeU(int iX,  int iY,  double& rho, Pnt<double>& u, double& uSqr)
    {
        uSqr = 0.;
        for(int iD=0; iD<d; ++iD)
        {
            u[iD] = 0.;
            for (int iF=0; iF<q; ++iF)
            {
                u[iD] += f[iF][iX][iY] * c[iD][iF];
            }
            //if(rho==0)
            //    u[iD]=0;
            //else
                u[iD] /= rho;
 
            uSqr += u[iD]*u[iD];
        }
    }
 
    double fEq1 (int iF, double rho, Pnt<double> u, double uSqr)
    {
        return rho*t[iF];
    }
 
    double fEq2(int iF, double rho, Pnt<double> u, double uSqr)
    {
        double c_u = 0.; // scalar product between c_{iF} and u
        for (int iD=0; iD<d; iD++)
        {
            c_u += c[iD][iF] * u[iD];
        }
        return rho*t[iF]*(1+c_u/csq);
    }
 
    double fEq4(int iF, double rho, Pnt<double> u, double uSqr)
    {
        double c_u = 0.; // scalar product between c_{iF} and u
        for (int iD=0; iD<d; iD++)
        {
            c_u += c[iD][iF] * u[iD];
        }
        return rho*t[iF]*(1. + c_u/csq + c_u*c_u/csq/csq/2.0 - uSqr/csq/2.0);
    }
 
    void FVon(int Ix , int Iy, double dx, double dt, double tau, int nb_min, int nb_max)
    {
        const int node=8;
        bool isBoundary = (nb_max-nb_min==node)?false:true ;
        double Edge_Mid[node+1][2];
        double SubVol_Cntr[node+1][2];
 
        double f_Edge_Mid[node+1];
        double f_SubVol_Cntr[node];
 
        double feq_Edge_Mid[node+1];
        double feq_SubVol_Cntr[node+1];
 
        double Ynorm[9][2];
        double Xnorm[9][2];
 
        int N;
        int Nb;
        int Nb_min=(Nb_min+node-1)%node+1;
        int Nb_max=(Nb_max+node-1)%node+1;;
 
        Pnt<double> u;
        double uSqr;
 
        Edge_Mid   [0 ][0]=P[Ix] [Iy][0];
        Edge_Mid   [0 ][1]=P[Ix] [Iy][1];
        SubVol_Cntr[0 ][0]=P[Ix] [Iy][0];
        SubVol_Cntr[0 ][1]=P[Ix] [Iy][1];
 
        for(int iF=0; iF<q ;iF++)
        {
            f_Edge_Mid   [0 ]=f[iF][Ix][Iy];
            f_SubVol_Cntr[0 ]=f[iF][Ix][Iy];
 
            u[0]=U[0][Ix][Iy];//+0.000012/omega;
            u[1]=U[1][Ix][Iy];
            uSqr=u[0]*u[0]+u[1]*u[1];
            feq_Edge_Mid   [0 ]=(this->*Feq)(iF, Rho[Ix][Iy] ,u , uSqr);
            feq_SubVol_Cntr[0 ]=(this->*Feq)(iF, Rho[Ix][Iy] ,u , uSqr);
 
            Ynorm[0][0]= Edge_Mid   [0     ][0] - Edge_Mid   [nb_min][0];
            Ynorm[0][1]= Edge_Mid   [nb_max][0] - Edge_Mid   [0     ][0];
 
            Xnorm[0][0]= Edge_Mid   [nb_min][1] - Edge_Mid   [0     ][1];
            Xnorm[0][1]= Edge_Mid   [0     ][1] - Edge_Mid   [nb_max][1];
 
            for(int nb=nb_min; nb<nb_max ;nb++)
            {
                N =(nb+node-1)%node+1;
                Nb=(nb+node)%node+1;
                Edge_Mid   [N ][0]=(    P[Ix            ] [Iy           ][0]+
                                        P[Ix+c[0][N ]   ] [Iy+c[1][N ]  ][0] ) /2.0;
                SubVol_Cntr[N ][0]=(    P[Ix            ] [Iy           ][0]+
                                        P[Ix+c[0][N ]   ] [Iy+c[1][N ]  ][0]+
                                        P[Ix+c[0][Nb]   ] [Iy+c[1][Nb]  ][0]
                                    ) /3.0;
 
                Edge_Mid   [N ][1]=(    P[Ix            ] [Iy           ][1]+
                                        P[Ix+c[0][N ]   ] [Iy+c[1][N ]  ][1] ) /2.0;
                SubVol_Cntr[N ][1]=(    P[Ix            ] [Iy           ][1]+
                                        P[Ix+c[0][N ]   ] [Iy+c[1][N ]  ][1]+
                                        P[Ix+c[0][Nb]   ] [Iy+c[1][Nb]  ][1]
                                    ) /3.0;
            }
 
            for(int nb=nb_min; nb<nb_max ;nb++)
            {
                N =(nb+node-1)%node+1;
                Nb=(nb+node)%node+1;
                f_Edge_Mid[N ]    =(f[iF][Ix][Iy]+f[iF][Ix+c[0][N ]][Iy+c[1][N ]])/2.0;
                f_SubVol_Cntr[N ]=(    f[iF][Ix          ][Iy         ] +
                                        f[iF][Ix+c[0][N ] ][Iy+c[1][N ]] +
                                        f[iF][Ix+c[0][Nb] ][Iy+c[1][Nb]]
                                    )/3.0;
 
                u[0]=U[0][Ix][Iy];//+0.000012/omega;
                u[1]=U[1][Ix][Iy];
                uSqr=u[0]*u[0]+u[1]*u[1];
                double f0=(this->*Feq)(iF, Rho[Ix][Iy] ,u , uSqr);
 
                u[0]=U[0][Ix+c[0][N ]][Iy+c[1][N ]];//+0.000012/omega;
                u[1]=U[1][Ix+c[0][N ]][Iy+c[1][N ]];
                uSqr=u[0]*u[0]+u[1]*u[1];
                double f1=(this->*Feq)(iF, Rho[Ix+c[0][N ]][Iy+c[1][N ]] ,u , uSqr);
 
                u[0]=U[0][Ix+c[0][Nb]][Iy+c[1][Nb]];//+0.000012/omega;
                u[1]=U[1][Ix+c[0][Nb]][Iy+c[1][Nb]];
                uSqr=u[0]*u[0]+u[1]*u[1];
                double f2=(this->*Feq)(iF, Rho[Ix+c[0][Nb]][Iy+c[1][Nb]] ,u , uSqr);
 
                feq_Edge_Mid[N ]    =(f0+f1)/2.0;
                feq_SubVol_Cntr[N ]=( f0 + f1 + f2 )/3.0;
            }
 
            //double Ynorm=xi - xi+1
            //double Xnorm=yi+1 - yi
            for(int nb=nb_min; nb<nb_max ;nb++)
            {
                N =(nb+node-1)%node+1;
                Nb=(nb+node)%node+1;
                Ynorm[N ][0]= Edge_Mid   [N ][0] - SubVol_Cntr[N ][0];
                Ynorm[N ][1]=-Edge_Mid   [Nb][0] + SubVol_Cntr[N ][0];
 
                Xnorm[N ][0]= SubVol_Cntr[N ][1] - Edge_Mid   [N ][1] ;
                Xnorm[N ][1]= Edge_Mid   [Nb][1] - SubVol_Cntr[N ][1];
            }
 
            double flux=0;
            for(int nb=nb_min; nb<nb_max ;nb++)
            {
                N =(nb+node-1)%node+1;
                Nb=(nb+node)%node+1;
 
                flux+=c[0][iF] * Xnorm[N ][0] * (f_Edge_Mid[N ]+ f_SubVol_Cntr[N ])/2.0;  //x
                flux+=c[1][iF] * Ynorm[N ][0] * (f_Edge_Mid[N ]+ f_SubVol_Cntr[N ])/2.0;  //y
 
                flux+=c[0][iF] * Xnorm[N ][1] * (f_Edge_Mid[Nb]+ f_SubVol_Cntr[N ])/2.0;  //x
                flux+=c[1][iF] * Ynorm[N ][1] * (f_Edge_Mid[Nb]+ f_SubVol_Cntr[N ])/2.0;  //y
 
            }
 
            if(isBoundary)
            {
                Nb=(nb_min+node)%node+1;
                flux+=c[0][iF] * Xnorm[0][0] * (f_Edge_Mid[Nb]+ f_SubVol_Cntr[0])/2.0;  //x
                flux+=c[1][iF] * Ynorm[0][0] * (f_Edge_Mid[Nb]+ f_SubVol_Cntr[0])/2.0;  //y
 
                Nb=(nb_max+node-1)%node+1;
                flux+=c[0][iF] * Xnorm[0][1] * (f_Edge_Mid[Nb]+ f_SubVol_Cntr[0])/2.0;  //x
                flux+=c[1][iF] * Ynorm[0][1] * (f_Edge_Mid[Nb]+ f_SubVol_Cntr[0])/2.0;  //y
            }
 
 
            double collision=0;
            for(int nb=nb_min; nb<nb_max ;nb++)
            {
                N =(nb+node-1)%node+1;
                Nb=(nb+node)%node+1;
 
                collision+= -dx/12.0/tau*(f_Edge_Mid   [0 ]+ f_Edge_Mid  [N ]+ f_SubVol_Cntr  [N ]
                                        -feq_Edge_Mid [0 ]- feq_Edge_Mid[N ]- feq_SubVol_Cntr[0 ])/3.0;
 
                collision+= -dx*dx/12.0/tau*(f_Edge_Mid   [0 ]+ f_Edge_Mid  [Nb]+ f_SubVol_Cntr  [N ]
                                        -feq_Edge_Mid [0 ]- feq_Edge_Mid[Nb]- feq_SubVol_Cntr[0 ])/3.0;
            }
 
            ftemp[iF][Ix][Iy]=f[iF][Ix][Iy]+dt/dx/dx*(collision-flux)+c[0][iF]/6.0*0.000012;
        }
 
    }
 
    void FVon(Pnt<int> ix, Pnt<int> iy, double dx, double dt, double tau)
    {
        for (int iX=ix[0]; iX<=ix[1]; ++iX)
            for (int iY=iy[0]; iY<=iy[1]; ++iY)
                FVon(iX , iY, dx, dt, tau,1,9);
 
        //for (int iX=ix[0]; iX<=ix[1]; ++iX)
            for (int iY=iy[0]; iY<=iy[1]; ++iY)  //w
                FVon(0 , iY, dx, dt, tau,7,8+3);
 
        //for (int iX=ix[0]; iX<=ix[1]; ++iX)
            for (int iY=iy[0]; iY<=iy[1]; ++iY)  //e
                FVon(nX-1 , iY, dx, dt, tau,3,7);
 
        for (int iX=ix[0]; iX<=ix[1]; ++iX)  //s
            //for (int iY=iy[0]; iY<=iy[1]; ++iY)
                FVon(iX , 0, dx, dt, tau,1,5);
 
        for (int iX=ix[0]; iX<=ix[1]; ++iX)  //n
            //for (int iY=iy[0]; iY<=iy[1]; ++iY)
                FVon(iX , nY-1, dx, dt, tau,5,8);
 
        for (int iX=ix[0]; iX<=ix[1]; ++iX)
            for (int iY=iy[0]; iY<=iy[1]; ++iY)
                for(int iF=0; iF<9 ;iF++)
                    f[iF][iX][iY]=ftemp[iF][iX][iY];
 
    }
 
    void collosion(double omega)
    {
        double uSqr;
        Pnt<double> u;
 
        for (int iX=0; iX<nX; ++iX)
        {
            for (int iY=0; iY<nY; ++iY)
            {
                u[0]=U[0][iX][iY]+0.000012/omega; //2.604e-5/omega;
                u[1]=U[1][iX][iY];
 
                uSqr=U[0][iX][iY]*U[0][iX][iY]+U[1][iX][iY]*U[1][iX][iY];
                uSqr=u[0]*u[0]+u[1]*u[1];
 
                for (int iF=0; iF<q; ++iF)
                {
                    f[iF][iX][iY] *= (1.-omega);
                    f[iF][iX][iY] += omega * (this->*Feq)(iF, Rho[iX][iY], u, uSqr);
                }
 
                for (int iF=1; iF<q; ++iF)
                {
                    if( c[0][iF]==-1 || (c[0][iF]==0 && c[1][iF] !=1) )
                    swap(f[iF][iX][iY], f[ops[iF]][iX][iY]);
                }
            }
        }
    }
 
    void streaming()
    {
        for (int iX=0; iX<nX; ++iX)
        {
            for (int iY=0; iY<nY; ++iY)
            {
                for (int iF=1; iF<q; ++iF)
                {
                    nextX = iX + c[0][iF];
                    nextY = iY + c[1][iF];
 
                    if( c[0][iF]==-1 || (c[0][iF]==0 && c[1][iF] !=1) )
                    if (nextX>=0 && nextX<nX && nextY>=0 && nextY<nY)
                    {
                        swap(f[ops[iF]][iX][iY], f[iF][nextX][nextY]);
                    }
                }
            }
        }
    }
 
    void update(int iX, int iY)
    {
        double ssum=0,usum=0,vsum=0;
        for (int iF=0; iF<q; ++iF)
        {
            ssum+=f[iF][iX][iY];
            usum+=f[iF][iX][iY]*c[0][iF];
            vsum+=f[iF][iX][iY]*c[1][iF];
        }
        Rho[iX][iY]=ssum;
        if (ssum==0)
        {
            U[0][iX][iY]=0;
            U[1][iX][iY]=0;
        }
        else
        {
            U[0][iX][iY]=usum/ssum;
            U[1][iX][iY]=vsum/ssum;
        }
    }
 
    void update(Pnt<int> ix, Pnt<int> iy)
    {
        for (int iX=ix[0]; iX<=ix[1]; ++iX)
            for (int iY=iy[0]; iY<=iy[1]; ++iY)
               update(iX,iY);
    }
 
    void update()
    {
        for (int iX=0; iX<nX; ++iX)
            for (int iY=0; iY<nY; ++iY)
               update(iX,iY);
    }
 
    void BounceBackBoundOn(int iX , int iY,  Pnt<int> IF)
    {
        for (int iF=1; iF<q; ++iF)
            //if( ( c[0][iF]==IF[0] && c[0][iF]!=0 ) || ( c[1][iF]==IF[1]  && c[1][iF]!=0 ) )
            if( Dir(iF, IF) )
                f[iF][iX][iY]=f[ops[iF]][iX][iY];
        update(iX,iY);
    }
 
    void BounceBackBoundOn(Pnt<int> ix, Pnt<int> iy,  Pnt<int> IF)
    {
        for (int iX=ix[0]; iX<=ix[1]; ++iX)
            for (int iY=iy[0]; iY<=iy[1]; ++iY)
               BounceBackBoundOn(iX ,iY, IF);
    }
 
    void OpenBoundOn(int iX , int iY,  Pnt<int> IF)
    {
        //for (int iF=0 ;iF<q; ++iF)
           // if( (c[0][iF]==-IF[0] && c[0][iF]!=0  ) || ( c[1][iF]==-IF[1] && c[1][iF]!=0 )   )
                //f[iF][iX][iY]=4.0/3.0*f [iF] [ iX+IF[0] ] [ iY+IF[1] ] - 1.0/3.0*f[iF] [ iX+2*IF[0] ] [ iY+2*IF[1]];
                //f[iF][iX][iY]=f[iF][iX+IF[0]][iY+IF[1]];
 
 
        for (int iF=0 ;iF<q; ++iF)
            if( (c[0][iF]==IF[0] && c[0][iF]!=0  ) || ( c[1][iF]==IF[1] && c[1][iF]!=0 )   )
                //f[iF][iX][iY]=4.0/3.0*f [iF] [ iX+IF[0] ] [ iY+IF[1] ] - 1.0/3.0*f[iF] [ iX+2*IF[0] ] [ iY+2*IF[1]];
                f[iF][iX][iY]=f[iF][iX+IF[0]][iY+IF[1]];
 
    update(iX,iY);
    for (int iD=0 ;iD<d; ++iD)
    if(IF[iD]==0)
        U[iD][iX][iY]=0;
    }
 
    void OpenBoundOn(Pnt<int> ix, Pnt<int> iy,  Pnt<int> IF)
    {
        for (int iX=ix[0]; iX<=ix[1]; ++iX)
            for (int iY=iy[0]; iY<=iy[1]; ++iY)
               OpenBoundOn(iX ,iY, IF);
    }
 
    void TecPrint(int t)
    {
        double rho;
        Pnt<double> u;
        //ofstream out ( "out.txt",ios::ate );
 
        double** str;
        str=new double*[nX];
        for(int i=0; i<nX ; i++)
            str[i]=new double [nY];
 
        for(int i=0; i<nX ; i++)
            for(int j=0; j<nY ; j++)
                str[i][j]=0;
 
        for(int i=0; i<nX ; i++)
            for(int j=0; j<nY ; j++)
                str[i][j]=IntgU( 0 , i , j) - IntgV(0 , 0 , i);
 
 
        out<<" VARIABLES = \"x\" \"y\" \"r\" \"u\" \"v\" \"s\" "<<endl;
        out<<" ZONE T= \"Z"<< t<<"\" "<<endl;
        out<<"I="<<nY<<", J="<<nX<<", K="<<1<<endl;
 
        for ( int iX=0;iX<nX;iX++ )
        {
            for ( int iY=0;iY<nY;iY++ )
            {
                computeRho(iX,iY,rho);
                computeU(iX, iY,rho, u, temp);
                out<<iX<<"\t"<<iY<<"\t"<<rho <<"\t"<<u[0]<<"\t"<<u[1]<<"\t"<<str[iX][iY] <<endl;
            }
            out<<endl;
        }
        //system ( "wgnuplot 3plot.plt -" );
    }
 
    bool Dir(int iF, Pnt<int> IF)
    {
        if(IF[0]*IF[1])
            return ( c[0][iF]!=-IF[0]  )  && ( c[1][iF]!=-IF[1])  ;
        else
            return (c[0][iF]==IF[0] && c[0][iF]!=0      ) || ( c[1][iF]==IF[1] && c[1][iF]!= 0 );
            //return (c2[0][iF]==IF[0] && c2[0][iF]!=-IF[0] ) || ( c2[1][iF]==IF[1] && c2[1][iF]!=-IF[1]);
    }
 
    double IntgV(int X0 , int Y0 , int X)
    {
        double I=0;
        for(int i=X0 ; i<X ; i++)
            I+= (U[1][i][Y0]+U[1][i+1][Y0])/2.0;
        return I;
    }
 
    double IntgU(int Y0 , int X , int Y)
    {
        double I=0;
        for(int j=Y0 ; j<Y ; j++)
            I+= (U[0][X][j]+U[0][X][j+1])/2.0;
        return I;
    }
};
 
 
int main()
{
    int q2;
    double* t2;
    int* ops2;
    int** c2;
    double csq2;
    int nX=32;
    int nY=nX;
 
    Q9v(q2,c2,t2,ops2,csq2 );
 
    double dx_lu=0.5;
    double dt_lu=0.001;
    double tau_lu=0.1;
    double visc=tau_lu/3.0;
 
    //double omega2=1./(visc/csq2+0.5); //for FDLBM
    double omega2=1/tau_lu;             //for FVLBM
 
    cout<< 1/omega2<<"    "<<endl;
    char syspause;
    cin>>syspause;
 
    Core_2D f2(nX,nY ,c2,t2,ops2,csq2,q2,4);
    f2.Init_Geo(nX,dx_lu);
 
    //init the field with rho=1  and u=0
    f2.Init (1, Pnt<double>(0.0,0.00) );
 
    for(int i=0;i<6000;i++)
    {
        cout<<i<<endl;
 
        //f2.collosion(omega2);  //for FDLBM
        //f2.streaming();        //for FDLBM
 
        f2.FVon(Pnt<int>(1,nX-2), Pnt<int>(1,nY-2), dx_lu, dt_lu, tau_lu);   //for FVLBM
 
        f2.update();
 
        f2.OpenBoundOn      (Pnt<int>(0,0)       , Pnt<int>(0,nY-1)    ,  Pnt<int>(1,0));  //w
        f2.BounceBackBoundOn(Pnt<int>(0,nX-1)    , Pnt<int>(nY-1,nY-1) ,  Pnt<int>(0,-1));  //n
        f2.BounceBackBoundOn(Pnt<int>(0,nX-1)    , Pnt<int>(0,0)       ,  Pnt<int>(0,+1));  //s
        f2.OpenBoundOn      (Pnt<int>(nX-1,nX-1) , Pnt<int>(0,nY-1)    ,  Pnt<int>(-1,0));  //e
 
 
        if(i%500==0) f2.TecPrint(i);
    }
    f2.TecPrint(400);
}
 
community/finite_volume_lbm_cell_vertex_scheme.txt · Last modified: 2011/12/17 17:33 by yann_at_flowkit
 
RSS - 2011 © FlowKit.com
Except where otherwise noted, text, pictures and animations on this site are licensed under a Creative Commons Attribution-Share Alike 3.0 License.