[selig5576]

I have recently starting making my numerical solver 3D from 2D and I have ran into a problem when solving my PPE. For my 2D code, I have a convergent solution and pressure field never exceeds values of 1.2. For my 3D solver after 1 iteration I get pressure values of 12000.

I am almost certain it has to do with my pressure solver as my velocity components are not reaching values even near those of the pressure. The method in which I am solving the PPE is SOR. My relaxation factor I am using is 1.93 and my tolerance is 1e-6. Also my max iterations is 30,000.

for (int j=2; j < ny-2; j++)
        {
            for (int k=2; k < nz-2; k++)
            {
                R[ID(1,j,k)] = (rho/dt)*(idx*(ufs[ID(1,j,k)] - u[ID(0,j,k)]) + idy*(vfs[ID(1,j,k)] - vfs[ID(1,j-1,k)]) + idz*(wfs[ID(1,j,k)] - wfs[ID(1,j,k-1)]));
                P[ID(1,j,k)] = (1.0-omega)*P[ID(1,j,k)] + Aw*(idx2*P[ID(2,j,k)] + idy2*(P[ID(1,j+1,k)] + P[ID(1,j-1,k)]) + idz2*(P[ID(1,j,k+1)] + P[ID(1,j,k-1)]) - R[ID(1,j,k)]);

                R[ID(nx-2,j,k)] = (rho/dt)*(idx*(u[ID(nx-1,j,k)] - ufs[ID(nx-3,j,k)]) + idy*(vfs[ID(nx-2,j,k)] - vfs[ID(nx-2,j-1,k)]) + idz*(wfs[ID(nx-2,j,k)] - wfs[ID(nx-2,j,k-1)]));
                P[ID(nx-2,j,k)] = (1.0-omega)*P[ID(nx-2,j,k)] + Ae*(idx2*P[ID(nx-3,j,k)] + idy2*(P[ID(nx-2,j+1,k)] + P[ID(nx-2,j,j-1)]) + idz2*(P[ID(nx-2,j,k+1)] + P[ID(nx-2,j,k-1)]) - R[ID(nx-2,j,k)]);
            }
        }

        for (int i=2; i < nx-2; i++)
        {
            for (int k=2; k < nz-2; k++)
            {
                R[ID(i,1,k)] = (rho/dt)*(idx*(ufs[ID(i,1,k)] - ufs[ID(i-1,1,k)]) + idy*(vfs[ID(i,1,k)] - v[ID(i,0,k)]) + idz*(wfs[ID(i,1,k)] - wfs[ID(i,1,k-1)]));
                P[ID(i,1,k)] = (1.0-omega)*P[ID(i,1,k)] + An*(idx2*(P[ID(i+1,1,k)] + P[ID(i-1,1,k)]) + idy2*P[ID(i,2,k)] + idz2*(P[ID(i,1,k+1)] + P[ID(i,1,k-1)]) - R[ID(i,1,k)]);

                R[ID(i,ny-2,k)] = (rho/dt)*(idx*(ufs[ID(i,ny-2,k)] - ufs[ID(i-1,ny-2,k)]) + idy*(v[ID(i,ny-1,k)] - vfs[ID(i,ny-3,k)]) + idz*(wfs[ID(i,ny-2,k)] - wfs[ID(i,ny-2,k-1)]));
                P[ID(i,ny-2,k)] = (1.0-omega)*P[ID(i,ny-2,k)] + As*(idx2*(P[ID(i+1,ny-2,k)] + P[ID(i-1,ny-2,k)]) + idy2*P[ID(i,ny-3,k)] + idz2*(P[ID(i,ny-2,k+1)] + P[ID(i,ny-2,k-1)]) - R[ID(i,ny-2,k)]);
            }
        }

        for (int i=2; i < nx-2; i++)
        {
            for (int j=2; j < ny-2; j++)
            {
                R[ID(i,j,1)] = (rho/dt)*(idx*(ufs[ID(i,j,1)] - ufs[ID(i-1,j,1)]) + idy*(vfs[ID(i,j,1)] - vfs[ID(i,j-1,1)]) + idz*(wfs[ID(i,j,1)] - w[ID(i,j,0)]));
                P[ID(i,j,1)] = (1.0-omega)*P[ID(i,j,1)] + Ab*(idx2*(P[ID(i+1,j,1)] + P[ID(i-1,j,1)]) + idy2*(P[ID(i,j+1,1)] + P[ID(i,j-1,1)]) + idz2*P[ID(i,j,2)] - R[ID(i,j,1)]);

                R[ID(i,j,nz-2)] = (rho/dt)*(idx*(ufs[ID(i,j,nz-2)] - ufs[ID(i-1,j,nz-2)]) + idy*(vfs[ID(i,j,nz-2)] - vfs[ID(i,j-1,nz-2)]) + idz*(w[ID(i,j,nz-1)] - wfs[ID(i,j,nz-3)]));
                P[ID(i,j,nz-2)] = (1.0-omega)*P[ID(i,j,nz-2)] + At*(idx2*(P[ID(i+1,j,nz-2)] + P[ID(i-1,j,nz-2)]) + idy2*(P[ID(i,j+1,nz-2)] + P[ID(i,j-1,nz-2)]) + idz2*P[ID(i,j,nz-3)] - R[ID(i,j,nz-2)]);
            }
        }