[unkn0wn]

Hello all,
We were trying to solve unsteady code for Lid Driven Cavity. Using Alternative Direction Implicit method with TDMA Algorithm. We are getting solution upto few time steps but after that omega is blowing up.

Here the code we are trying.
Any ideas whats going wrong in our code.

#include<iostream>
#include<cstdlib>
#include<cmath>
#include<vector>
#include<cmath>
#include<fstream>
#include<string>

#define C_ERROR 0.000000001
#define grid 101

using namespace std;

vector< vector<double> > psi;
vector< vector<double> > omega;
vector< vector<double> > u;
vector< vector<double> > v;

double dx=0, dy=0, beta=0, U=0;
double dt=0, h=0;

void Initialise(){
    psi.resize(grid, vector<double>(grid,0));
    omega.resize(grid, vector<double>(grid,0));
    u.resize(grid, vector<double>(grid,0));
    v.resize(grid, vector<double>(grid,0));

    for(uint i=0;i<u[0].size();i++)
        u[0][i]=1;

    dx=dy=(double)1/grid;
    beta=dx/dy;
    U=1;
    h=(double)1/grid;
    dt=0.001;
}

void Print(const vector< vector<double> > &v){
    for(uint i=0;i<v.size();i++){
        for(uint j=0;j<v[i].size();j++)
            cout<<"["<<i<<"]"<<"["<<j<<"]"<<v[i][j]<<"\t";
        cout<<endl;
    }
    cout<<endl;
}

void PrintSingle(const vector< double > &v){
    for(uint i=0;i<v.size();i++){
            cout<<"["<<i<<"]"<<v[i]<<"\t";
    }
    cout<<endl<<endl;
}

void SetBoundaryConditions(){
    uint i=0;
    for(i=1;i<grid-1;i++){
        omega[0][i]=-(2*psi[1][i])/(dx*dx);
        omega[grid-1][i]=-(2*psi[grid-2][i])/(dx*dx);
        omega[i][0]=-(2*psi[i][1])/(dy*dy);
        omega[i][grid-1]=-(2*psi[i][grid-2]+2*dy*U)/(dy*dy);
    }
}

double Error(vector< vector<double> > a, vector< vector<double> > b){
    double diffSum=0, nSum=0;
    for(uint i=0;i<a.size();i++)
        for(uint j=0;j<b.size();j++) {
        nSum+=fabs(a[i][j]);
        diffSum+=fabs(a[i][j]-b[i][j]);
    }
    return diffSum/nSum;
}


void thomas(vector<double>& a,
            vector<double>& b,
            vector<double>& c,
            vector<double>& d,
            vector<double>& f) {

     size_t N = d.size();

     for (int i=1; i<N; i++) {
        b[i]=b[i]-c[i-1]*a[i]/b[i-1];

        d[i]=d[i]-d[i-1]*a[i]/b[i-1];
      }

     f[N-1]=d[N-1]/b[N-1];

     for (int i=N-1; i-- > 0; ) {
        f[i] = (d[i]-c[i]*f[i+1])/b[i];
      }

}

void WriteFile(double a,uint Re){
    ofstream File;
    File.open("unsteady.dat", ios::out | ios::app);
    for(uint i=0;i<grid;i++){
        for(uint j=0;j<grid;j++)
            File<<a<<" "<<(i+1)*dx<<" "<<(j+1)*dy<<" "<<psi[i][j]<<" "<<omega[i][j]<<" "<<u[i][j]<<" "<<v[i][j]<<endl;
    }
}

void Solve(double Re){
    uint i=0, j=0, iter=0, mainiter=0;
    double error=0;
    vector< vector<double> > tempomega;
    vector< vector<double> > temppsi;


    do{
    SetBoundaryConditions();
    tempomega=omega;
    mainiter++;

        do{
            //Computing PSI
            temppsi=psi;
            for(i=1;i<psi.size()-1;i++)
                for(j=1;j<psi[i].size()-1;j++)
                    psi[i][j]=((psi[i+1][j]+psi[i-1][j]+beta*beta*(psi[i][j+1]+psi[i][j-1])+(dx*dx)*omega[i][j]))/(2*(1+beta*beta));
        }while(Error(psi,temppsi)>C_ERROR);

    //Print(psi);

    //Solving u and v
    for(i=1;i<grid-1;i++){
        for(j=1;j<grid-1;j++){
            u[i][j]=-(psi[i][j+1]-psi[i][j-1])/(2*dy);
            v[i][j]=(psi[i+1][j]-psi[i-1][j])/(2*dy);
        }
    }

    vector<double> a(grid-2,0);
    vector<double> b(grid-2,0);
    vector<double> c(grid-2,0);
    vector<double> d(grid-2,0);
    vector<double> x(grid-2,0);

    for(j=1;j<grid-1;j++){

        for(i=1;i<grid-1;i++){
            a[i-1]=(-u[i][j]/h - 1/(Re*h*h));
            b[i-1]=((1/dt)+(2/(Re*h*h)));
            c[i-1]=((u[i][j]/h)-(1/(Re*h*h)));
            d[i-1]=(-v[i][j]*((omega[i][j+1]-omega[i][j-1])/h) + (omega[i][j+1]+omega[i][j-1]-2*omega[i][j])/(Re*h*h) + (omega[i][j]*2)/dt);
        }
        omega[0][j]=(-v[0][j]*((omega[0][j+1]-omega[0][j-1])/h) + (omega[0][j+1]+omega[0][j-1]-2*omega[0][j])/(Re*h*h) + (omega[0][j]*2)/dt);;
        omega[grid-1][j]=(-v[grid-1][j]*((omega[grid-1][j+1]-omega[grid-1][j-1])/h) + (omega[grid-1][j+1]+omega[grid-1][j-1]-2*omega[grid-1][j])/(Re*h*h) + (omega[grid-1][j]*2)/dt);

        d[0]=d[0]-a[0]*(-v[0][j]*((omega[0][j+1]-omega[0][j-1])/h) + (omega[0][j+1]+omega[0][j-1]-2*omega[0][j])/(Re*h*h) + (omega[0][j]*2)/dt);
        d[grid-3]=d[grid-3]-c[grid-3]*(-v[grid-1][j]*((omega[grid-1][j+1]-omega[grid-1][j-1])/h) + (omega[grid-1][j+1]+omega[grid-1][j-1]-2*omega[grid-1][j])/(Re*h*h) + (omega[grid-1][j]*2)/dt);

        a[0]=0;
        c[grid-3]=0;

        thomas(a,b,c,d,x);

        for(i=1;i<grid-1;i++){
            omega[i][j]=x[i-1];
        }
    }

    fill(a.begin(), a.end(), 0);
    fill(b.begin(), b.end(), 0);
    fill(c.begin(), c.end(), 0);
    fill(d.begin(), d.end(), 0);
    fill(x.begin(), x.end(), 0);

    for(i=1;i<grid-1;i++){
        for(j=1;j<grid-1;j++){
            a[j-1]=(-v[i][j]/h)-(1/(Re*h*h));
            b[j-1]=((1/dt)+(2/(Re*h*h)));
            c[j-1]=(v[i][j]/h) -(1/(Re*h*h));
            d[j-1]=(-u[i][j]*(omega[i+1][j]-omega[i-1][j]))/h + (omega[i+1][j]+omega[i-1][j]-2*omega[i][j])/(Re*h*h) + (2*omega[i][j])/dt;
        }

        omega[i][0]=(-u[i][0]*(omega[i+1][0]-omega[i-1][0]))/h + (omega[i+1][0]+omega[i-1][0]-2*omega[i][0])/(Re*h*h) + (2*omega[i][0])/dt;
        omega[i][grid-1]=(-u[i][grid-1]*(omega[i+1][grid-1]-omega[i-1][grid-1]))/h + (omega[i+1][grid-1]+omega[i-1][grid-1]-2*omega[i][grid-1])/(Re*h*h) + (2*omega[i][grid-1])/dt;

        d[0]=d[0]-a[0]*(-u[i][0]*(omega[i+1][0]-omega[i-1][0]))/h + (omega[i+1][0]+omega[i-1][0]-2*omega[i][0])/(Re*h*h) + (2*omega[i][0])/dt;
        d[grid-3]=d[grid-3]-c[grid-3]*(-u[i][grid-1]*(omega[i+1][grid-1]-omega[i-1][grid-1]))/h + (omega[i+1][grid-1]+omega[i-1][grid-1]-2*omega[i][grid-1])/(Re*h*h) + (2*omega[i][grid-1])/dt;

        a[0]=0;
        c[grid-3]=0;

        thomas(a,b,c,d,x);

        for(j=1;j<grid-1;j++){
            omega[i][j]=x[j-1];
        }
    }
    Print(omega);
    WriteFile(mainiter*dt,Re);
    cout<<"Time Step : "<<mainiter*dt<<endl;

    }while(Error(omega,tempomega)>C_ERROR);
}


int main(){
    Initialise();

    cout<<"Switching to Re=400\n";
    Solve(100);
    
/*
//Test Case to cCheck Thomas Algorithm
    vector<double> a(4,0);
    vector<double> b(4,0);
    vector<double> c(4,0);
    vector<double> d(4,0);
    vector<double> x(4,0);

    d[0]=8;
    d[1]=116;
    d[2]=47;
    d[3]=40;



    a[0]=0;
    a[1]=4;
    a[2]=2;
    a[3]=4;

    b[0]=2;
    b[1]=5;
    b[2]=5;
    b[3]=7;

    c[0]=3;
    c[1]=34;
    c[2]=7;
    c[3]=0;


    thomas(a,b,c,d,x);

    cout<<"A: "<<a[0]<<" "<<a[1]<<" "<<a[2]<<" "<<a[3]<<"\n";
    cout<<"B: "<<b[0]<<" "<<b[1]<<" "<<b[2]<<" "<<b[3]<<"\n";
    cout<<"C: "<<c[0]<<" "<<c[1]<<" "<<c[2]<<" "<<c[3]<<"\n";
    cout<<"D: "<<d[0]<<" "<<d[1]<<" "<<d[2]<<" "<<d[3]<<"\n";
    cout<<"X: "<<x[0]<<" "<<x[1]<<" "<<x[2]<<" "<<x[3]<<"\n";
*/


}